simulate_burninthining_single: simulate_burninthining_single

In isci1102/ERPM:

Description

Une observation de

Usage

simulate_burninthining_single(partition, theta, nodes, effects, objects, num.steps, neighborhood, sizes.allowed, sizes.simulated, max.thining)

Arguments

partition	Observed partition.
theta	Initial model parameters
nodes	Data frame containing the nodes.
effects	Effects or sufficient statistics. A list with a vector "names" and a vector "objects".
objects	Objects used for statistics calculation. A list with a vector "name" and a vector "object".
num.steps	Number of samples wanted.
neighborhood	Way of choosing partitions: probability vector (probability actors swap, probability merge/division, probability single actor move).
sizes.allowed	Vector of group sizes allowed in sampling.
sizes.simulated	Vector of group sizes allowed in the Markov chain but not necessraily sampled.
max.thining	Where to stop adding thining.

Examples

##---- Should be DIRECTLY executable !! ----
##-- ==>  Define data, use random,
##--	or do  help(data=index)  for the standard data sets.

## The function is currently defined as
function (partition, theta, nodes, effects, objects, num.steps,
    neighborhood, sizes.allowed, sizes.simulated, max.thining)
{
    num.effects <- length(effects$names)
    chain <- draw_Metropolis_single(theta, partition, nodes,
        effects, objects, 1, 1, num.steps, neighborhood, sizes.allowed,
        sizes.simulated, return.all.partitions = F)
    ok <- F
    current.draws <- chain$draws
    current.last <- chain$last.partition
    current.thining <- 1
    allautocors <- c()
    while (!ok) {
        draws <- current.draws[seq(1, nrow(current.draws), current.thining),
            ]
        autocors <- rep(0, num.effects)
        for (e in 1:num.effects) {
            autocors[e] <- cor(draws[1:(num.steps - 1), e], draws[2:num.steps,
                e])
        }
        if (current.thining%%50 == 0) {
            print("thining")
            print(current.thining)
            print("autocorrelations")
            print(autocors)
        }
        allautocors <- rbind(allautocors, autocors)
        if (current.thining >= max.thining) {
            ok <- T
        }
        else {
            new.chain <- draw_Metropolis_single(theta, current.last,
                nodes, effects, objects, 1, 1, num.steps, neighborhood,
                sizes.allowed, sizes.simulated, return.all.partitions = F)
            current.draws <- rbind(current.draws, new.chain$draws)
            current.last <- new.chain$last.partition
            current.thining <- current.thining + 1
        }
    }
    allmeans <- matrix(0, nrow(draws), num.effects)
    allmeans[1, ] <- draws[1, ]
    for (d in 2:nrow(draws)) {
        allmeans[d, ] <- colMeans(draws[1:d, ])
    }
    smoothedautocors <- allautocors
    smoothedmeans <- allmeans
    for (eff in 1:num.effects) {
        lo <- loess(y ~ x, data.frame(x = 1:max.thining, y = allautocors[,
            eff]))
        smoothedautocors[!is.na(smoothedautocors[, eff]), eff] <- lo$fitted
        lo <- loess(y ~ x, data.frame(x = 1:nrow(draws), y = allmeans[,
            eff]))
        smoothedmeans[!is.na(smoothedmeans[, eff]), eff] <- lo$fitted
    }
    return(list(draws = current.draws, autocorrelations = allautocors,
        moving.means = allmeans, autocorrelations.smoothed = smoothedautocors,
        moving.means.smoothed = smoothedmeans))
  }

isci1102/ERPM documentation built on Jan. 18, 2022, 12:25 a.m.