estimate_multipleERPM: estimate_multipleERPM
In isci1102/ERPM:

estimate_multipleERPM(partitions, presence.tables, nodes, objects, effects, startingestimates, multiplicationfactor = 30, gainfactor = 0.1, a.scaling = 0.2, r.truncation.p1 = 2, r.truncation.p2 = 5, burnin = 30, thining = 10, length.p1 = 100, min.iter.p2 = NULL, max.iter.p2 = NULL, multiplication.iter.p2 = 200, num.steps.p2 = 6, length.p3 = 1000, neighborhood = c(0.7, 0.3, 0), fixed.estimates = NULL, sizes.allowed = NULL, sizes.simulated = NULL, double.averaging = F, inv.zcov = NULL, inv.scaling = NULL, parallel = F, parallel2 = F, cpus = 1)

`partitions`	Observed partitions.
`presence.tables`	Matrix indicating which actors were present for each observations (mandatory).
`nodes`	Data frame containing the nodes.
`objects`	Objects used for statistics calculation. A list with a vector "name" and a vector "object".
`effects`	Effects or sufficient statistics. A list with a vector "names" and a vector "objects".
`startingestimates`	First guess for the model parameters.
`multiplicationfactor`	USELESS. Default=30.
`gainfactor`	Numeric used to decrease the size of steps made in the Newton optimization. Default = 0.1
`a.scaling`	Numeric used to reduce the influence of non-diagonal elements in the scaling matrix for stability. Default = 0.2
`r.truncation.p1`	Numeric used to limit extreme values in the covariance matrix for stability. Default = 2
`r.truncation.p2`	Numeric used to limit extreme values in the covariance matrix for stability. Default = 5
`burnin`	Integer for the number of burn-in steps before sampling. Default = 30
`thining`	Integer for the number of thining steps between sampling.Default = 10
`length.p1`	Number of samples in phase 1. Default = 100.
`min.iter.p2`	Minimum number of sub-steps in phase 2. Default = NULL
`max.iter.p2`	Maximum number of substeps in phase 2. Default = NULL
`multiplication.iter.p2`	value for the lengths of sub-steps in phase 2 (multiplied by 2.52^k). Default = 100
`num.steps.p2`	Number of optimisation steps in phase 2. Dedault = 6
`length.p3`	Number of samples in phase 3. Default = 1000
`neighborhood`	Way of choosing partitions: probability vector (actors swap, merge/division, single actor move). Default = (0.7,0.3,0)
`fixed.estimates`	If some parameters are fixed, list with as many elements as effects, these elements equal a fixed value if needed, or NULL if they should be estimated. Default = NULL
`sizes.allowed`	Vector of group sizes allowed in sampling. Default = NULL
`sizes.simulated`	Vector of group sizes allowed in the Markov chain but not necessraily sampled (now, it only works for vectors like size_min:size_max). Default = NULL
`double.averaging`	Option to average the statistics sampled in each sub-step of phase 2. Default = False
`inv.zcov`	Initial value of the inverted covariance matrix (if a phase 3 was run before) to bypass the phase 1. Default = NULL
`inv.scaling`	Initial value of the inverted scaling matrix (if a phase 3 was run before) to bypass the phase 1. Default = NULL
`parallel`	Whether the phase 1 and 3 should be parallelized. Default = False
`parallel2`	Whether there should be several phases 2 run in parallel. Default = False
`cpus`	How many cores can be used. Default = 1

Marion Hoffman

##---- 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 (partitions, presence.tables, nodes, objects, effects,
    startingestimates, multiplicationfactor = 30, gainfactor = 0.1,
    a.scaling = 0.2, r.truncation.p1 = 2, r.truncation.p2 = 5,
    burnin = 30, thining = 10, length.p1 = 100, min.iter.p2 = NULL,
    max.iter.p2 = NULL, multiplication.iter.p2 = 200, num.steps.p2 = 6,
    length.p3 = 1000, neighborhood = c(0.7, 0.3, 0), fixed.estimates = NULL,
    sizes.allowed = NULL, sizes.simulated = NULL, double.averaging = F,
    inv.zcov = NULL, inv.scaling = NULL, parallel = F, parallel2 = F,
    cpus = 1)
{
    z.obs <- rowSums(computeStatistics_multiple(partitions, presence.tables,
        nodes, effects, objects))
    gainfactors <- rep(0, num.steps.p2)
    for (i in 1:num.steps.p2) {
        gainfactors[i] <- gainfactor/(2^(i - 1))
    }
    num.effects <- length(effects$names)
    if (!is.null(fixed.estimates)) {
        for (e in 1:num.effects) {
            if (!is.null(fixed.estimates[[e]])) {
                startingestimates[e] <- fixed.estimates[[e]]
            }
        }
    }
    print("Observed statistics")
    print(z.obs)
    print("Burn-in")
    print(burnin)
    print("Thining")
    print(thining)
    if (!is.null(inv.zcov)) {
        estimates.phase1 <- startingestimates
        autocorrelations.phase1 <- NULL
    }
    else {
        results.phase1 <- run_phase1_multiple(partitions, startingestimates,
            z.obs, presence.tables, nodes, effects, objects,
            burnin, thining, gainfactor, a.scaling, r.truncation.p1,
            length.p1, neighborhood, fixed.estimates, sizes.allowed,
            sizes.simulated, parallel, cpus)
        estimates.phase1 <- results.phase1$estimates
        inv.zcov <- results.phase1$inv.zcov
        inv.scaling <- results.phase1$inv.scaling
        autocorrelations.phase1 <- results.phase1$autocorrelations
    }
    if (parallel2) {
        sfExport("partitions", "estimates.phase1", "inv.zcov",
            "inv.scaling", "z.obs", "presence.tables", "nodes",
            "effects", "objects", "burnin", "thining", "num.steps.p2",
            "gainfactors", "r.truncation.p2", "min.iter.p2",
            "max.iter.p2", "neighborhood", "fixed.estimates",
            "sizes.allowed", "sizes.simulated", "double.averaging")
        res <- sfLapply(1:cpus, fun = function(k) {
            set.seed(k)
            subres <- run_phase2_multiple(partitions, estimates.phase1,
                inv.zcov, inv.scaling, z.obs, presence.tables,
                nodes, effects, objects, burnin, thining, num.steps.p2,
                gainfactors, r.truncation.p2, min.iter.p2, max.iter.p2,
                multiplication.iter.p2, neighborhood, fixed.estimates,
                sizes.allowed, sizes.simulated, double.averaging)
            return(subres)
        })
        final.estimates <- c()
        for (k in 1:cpus) final.estimates <- final.estimates +
            res[[k]]$final.estimates
        estimates.phase2 <- final.estimates/cpus
    }
    else {
        results.phase2 <- run_phase2_multiple(partitions, estimates.phase1,
            inv.zcov, inv.scaling, z.obs, presence.tables, nodes,
            effects, objects, burnin, thining, num.steps.p2,
            gainfactors, r.truncation.p2, min.iter.p2, max.iter.p2,
            multiplication.iter.p2, neighborhood, fixed.estimates,
            sizes.allowed, sizes.simulated, double.averaging)
        estimates.phase2 <- results.phase2$final.estimates
    }
    results.phase3 <- run_phase3_multiple(partitions, estimates.phase2,
        z.obs, presence.tables, nodes, effects, objects, burnin,
        thining, a.scaling, length.p3, neighborhood, sizes.allowed,
        sizes.simulated, fixed.estimates, parallel, cpus)
    draws <- results.phase3$draws
    means <- results.phase3$means
    standard.deviations <- results.phase3$standard.deviations
    standard.errors <- results.phase3$standard.errors
    convergence.ratios <- results.phase3$convergence.ratios
    autocorrelations.phase3 <- results.phase3$autocorrelations
    results <- data.frame(effect = effects$names, object = effects$objects,
        est = as.vector(estimates.phase2), std.err = standard.errors,
        conv = convergence.ratios)
    print_results(results)
    objects.phase1 <- list(autocorrelations = autocorrelations.phase1)
    objects.phase2 <- list(estimates = results.phase2$all.estimates,
        lengths.subphases = results.phase2$lengths.subphases)
    objects.phase3 <- list(draws = draws, means = means, standard.deviations = standard.deviations,
        inv.zcov = inv.zcov, inv.scaling = inv.scaling, autocorrelations = autocorrelations.phase3)
    return(list(results = results, objects.phase1 = objects.phase1,
        objects.phase2 = objects.phase2, objects.phase3 = objects.phase3))
  }