library("networksis") bipartite.graph <- c(1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0) bipartite.graph <- matrix(bipartite.graph, nrow = 3, byrow = TRUE) example.net <- network(bipartite.graph) example.net %v% "set" <- c(rep(1, 3), rep(2, 4)) color <- shape <- example.net %v% "set" plot(example.net, vertex.col = color, vertex.sides = 2 + shape, vertex.cex = 2.5) sim <- simulate(example.net, nsim = 1000) exp(sim %n% "log.graphspace.size") exp(sim %n% "log.graphspace.SE") data("finch") plot(finch, vertex.col = c(rep(2, 13), rep(3, 17)), vertex.cex = 2.5) sim <- simulate(finch, nsim = 10000) exp(sim %n% "log.graphspace.size") exp(sim %n% "log.graphspace.SE") importance.weights <- 1 / exp(sim %n% "log.prob") hist(importance.weights, breaks = 75, xlab = "Inverse Graph Probability", main = "") nsim <- 100000 prob.vec <- rep(0, nsim) s.bar.squared.vec <- rep(0, nsim) for(i in 1:nsim) { sim <- simulate(finch, nsim = 1) new.graph <- as.matrix.network(sim) s.bar.squared <- (sum((new.graph %*% t(new.graph)) ^ 2) - sum(diag((new.graph %*% t(new.graph)) ^ 2))) / (13 * 12) s.bar.squared.vec[i] <- s.bar.squared prob.vec[i] <- exp(sim %n% "log.prob") } nbreaks <- 75 w <- (1 / prob.vec) / (sum(1 / prob.vec)) intervals <- cut(s.bar.squared.vec, nbreaks, include.lowest = TRUE) weights <- sapply(split(w, f = intervals), sum) x <- seq(min(s.bar.squared.vec), max(s.bar.squared.vec), length.out = nbreaks) plot(x, weights, type = "h", xlab = "", ylab = "Probability Density", lwd = 2) abline(v = 53.115) sim <- simulate_sis(finch ~ coincidences(0:17), nsim = 10000) observed.stats <- summary(finch ~ coincidences(0:17)) sampled.stats <- sim %n% "samplestatistics" library("Hmisc") p <- exp(sim %n% "log.prob") p <- p / sum(p) maxs <- apply(sampled.stats, 2, wtd.quantile, weights = p, probs = 0.975, normwt = TRUE) mins <- apply(sampled.stats, 2, wtd.quantile, weights = p, probs = 0.025, normwt = TRUE) means <- apply(sampled.stats, 2, wtd.mean, weights = p) plot(0:17, means, type = "b", ylim = c(0, 24.5), lwd = 3, lty = 3, xlab = "Number of Islands", ylab = "Pairs of Finches") for(i in 1:18) { points(rep(i - 1, 2), c(maxs[i], mins[i]), type = "l", lwd = 2) } points(0:17, observed.stats, type = "b", pch = 4, lwd = 3) r0 <- (p %*% sweep(sampled.stats, 2, observed.stats, "<"))[1, ] r1 <- (p %*% sweep(sampled.stats, 2, observed.stats, ">"))[1, ] round(apply(cbind(r0, r1), 1, min), digits = 8)[1] threshold <- 10 ^ (-18) breakloop <- 0 while(breakloop == 0) { sissim <- simulate(finch, nsim = 1) if(sissim %n% "log.prob" < log(threshold)){breakloop <- 1} } sim <- ergm(finch ~ coincidences(0), constraints = ~degrees, theta0 = 0.6108, MCMCsamplesize = 10000, control = control.ergm(initial.network = sissim)) sim <- simulate_sis(finch ~ coincidences(0:17), control = simulate.control(parallel = 4)), nsim = 10000)