[Rsiena-commits] r186 - in pkg: RSiena RSiena/R RSiena/inst/doc RSiena/man RSiena/src RSiena/src/model RSiena/src/model/ml RSienaTest RSienaTest/R RSienaTest/doc RSienaTest/inst/doc RSienaTest/inst/examples RSienaTest/man RSienaTest/src RSienaTest/src/model RSienaTest/src/model/ml
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Dec 4 18:51:37 CET 2011
Author: ripleyrm
Date: 2011-12-04 18:51:35 +0100 (Sun, 04 Dec 2011)
New Revision: 186
Added:
pkg/RSienaTest/R/algorithms.r
pkg/RSienaTest/man/algorithms.Rd
pkg/RSienaTest/man/profileLikelihoods.Rd
Removed:
pkg/RSienaTest/inst/examples/algorithms.r
Modified:
pkg/RSiena/DESCRIPTION
pkg/RSiena/R/bayes.r
pkg/RSiena/R/initializeFRAN.r
pkg/RSiena/R/maxlikec.r
pkg/RSiena/R/phase1.r
pkg/RSiena/R/phase2.r
pkg/RSiena/R/phase3.r
pkg/RSiena/R/robmon.r
pkg/RSiena/R/simstatsc.r
pkg/RSiena/changeLog
pkg/RSiena/inst/doc/RSiena_Manual.pdf
pkg/RSiena/man/RSiena-package.Rd
pkg/RSiena/man/bayes.Rd
pkg/RSiena/man/maxlikec.Rd
pkg/RSiena/man/simstats0c.Rd
pkg/RSiena/src/model/EpochSimulation.cpp
pkg/RSiena/src/model/EpochSimulation.h
pkg/RSiena/src/model/Model.cpp
pkg/RSiena/src/model/Model.h
pkg/RSiena/src/model/State.cpp
pkg/RSiena/src/model/ml/Chain.cpp
pkg/RSiena/src/model/ml/MLSimulation.cpp
pkg/RSiena/src/model/ml/MLSimulation.h
pkg/RSiena/src/siena07models.cpp
pkg/RSiena/src/siena07models.h
pkg/RSiena/src/siena07setup.cpp
pkg/RSiena/src/siena07utilities.cpp
pkg/RSienaTest/DESCRIPTION
pkg/RSienaTest/NAMESPACE
pkg/RSienaTest/R/bayes.r
pkg/RSienaTest/R/initializeFRAN.r
pkg/RSienaTest/R/maxlikec.r
pkg/RSienaTest/R/phase1.r
pkg/RSienaTest/R/phase2.r
pkg/RSienaTest/R/phase3.r
pkg/RSienaTest/R/robmon.r
pkg/RSienaTest/R/simstatsc.r
pkg/RSienaTest/changeLog
pkg/RSienaTest/doc/RSiena_Manual.tex
pkg/RSienaTest/inst/doc/RSiena_Manual.pdf
pkg/RSienaTest/inst/examples/runalg.r
pkg/RSienaTest/man/RSiena-package.Rd
pkg/RSienaTest/man/bayes.Rd
pkg/RSienaTest/man/maxlikec.Rd
pkg/RSienaTest/man/simstats0c.Rd
pkg/RSienaTest/src/model/EpochSimulation.cpp
pkg/RSienaTest/src/model/EpochSimulation.h
pkg/RSienaTest/src/model/Model.cpp
pkg/RSienaTest/src/model/Model.h
pkg/RSienaTest/src/model/State.cpp
pkg/RSienaTest/src/model/ml/Chain.cpp
pkg/RSienaTest/src/model/ml/MLSimulation.cpp
pkg/RSienaTest/src/model/ml/MLSimulation.h
pkg/RSienaTest/src/siena07models.cpp
pkg/RSienaTest/src/siena07models.h
pkg/RSienaTest/src/siena07setup.cpp
pkg/RSienaTest/src/siena07utilities.cpp
Log:
Alter interface for bays and algorithms
Modified: pkg/RSiena/DESCRIPTION
===================================================================
--- pkg/RSiena/DESCRIPTION 2011-11-27 16:23:01 UTC (rev 185)
+++ pkg/RSiena/DESCRIPTION 2011-12-04 17:51:35 UTC (rev 186)
@@ -1,15 +1,15 @@
Package: RSiena
Type: Package
Title: Siena - Simulation Investigation for Empirical Network Analysis
-Version: 1.0.12.185
-Date: 2011-11-27
+Version: 1.0.12.186
+Date: 2011-12-04
Author: Various
Depends: R (>= 2.10.0)
Imports: Matrix
Suggests: tcltk, snow, rlecuyer, network, codetools, lattice, MASS, parallel,
xtable, tools
SystemRequirements: GNU make, tcl/tk 8.5, Tktable
-Maintainer: Ruth Ripley <ruth at stats.ox.ac.uk>
+Maintainer: RSiena developers <rsiena at stats.ox.ac.uk>
Description: Fits models to longitudinal networks
License: GPL (>=2)
LazyLoad: yes
Modified: pkg/RSiena/R/bayes.r
===================================================================
--- pkg/RSiena/R/bayes.r 2011-11-27 16:23:01 UTC (rev 185)
+++ pkg/RSiena/R/bayes.r 2011-12-04 17:51:35 UTC (rev 186)
@@ -12,7 +12,8 @@
##@bayes Bayesian fit a Bayesian model
bayes <- function(data, effects, model, nwarm=100, nmain=100, nrunMHBatches=20,
plotit=FALSE, nbrNodes=1, dfra=NULL, n=10,
- priorSigma=NULL, prevAns=NULL, getDocumentation=FALSE)
+ priorSigma=NULL, prevAns=NULL, clusterType=c("PSOCK", "FORK"),
+ getDocumentation=FALSE)
{
##@createStores internal bayes Bayesian set up stores
createStores <- function()
@@ -22,7 +23,7 @@
z$posteriorTot <<- matrix(0, nrow=z$nGroup, ncol=npar)
z$posteriorMII <<- array(0, dim=c(z$nGroup, npar, npar))
z$candidates <<- array(NA, dim=c(numberRows, z$nGroup, npar))
- z$acceptances <<- matrix(NA, nrow=z$nGroup, ncol=numberRows)
+ z$acceptances <<- matrix(NA, ncol=z$nGroup, nrow=numberRows)
z$MHacceptances <<- array(NA, dim=c(numberRows, z$nGroup,
z$nDependentVariables, 9))
z$MHrejections <<- array(NA, dim=c(numberRows, z$nGroup,
@@ -37,7 +38,7 @@
nrun <- nrow(z$parameters)
npar <- length(z$theta)
end <- start + nrun - 1
- z$acceptances[, start:end] <<- z$accepts
+ z$acceptances[start:end, ] <<- z$accepts
z$candidates[start:end,, ] <<- z$parameters
z$posteriorTot <<- z$posteriorTot + colSums(z$parameters)
for (group in 1:z$nGroup)
@@ -99,9 +100,8 @@
}
cat('fine tuning took ', iter, ' iterations. Scalefactor:',
z$scaleFactors, '\n')
-
}
- ##@MCMCcycle algorithms do some loops of (MH steps and sample parameters)
+ ##@MCMCcycle internal Bayes do some loops of (MH steps and sample parameters)
MCMCcycle <- function(nrunMH, nrunMHBatches, change=TRUE)
{
z$accepts <<- matrix(NA, nrow=z$nGroup, nrunMHBatches)
@@ -116,15 +116,16 @@
z$nrunMH <<- nrunMH
for (i in 1:nrunMHBatches)
{
- # c1 <- proc.time()[1]
- ans <- z$FRAN(z, returnChains=TRUE, byGroup=TRUE)
+ # cc <- proc.time()[1]
+
+ ans <- z$FRAN(z, byGroup=TRUE, returnLoglik=TRUE, onlyLoglik=TRUE)
# c2 <- proc.time()[1]
- # cat ('fran',c2-c1,'\n')
- z$chain <<- ans$chain
- # c1 <- proc.time()[1]
+ # cat ('fran',c2-cc,'\n')
+ z$loglik <<- ans$loglik
+ # cc <- proc.time()[1]
sampleParameters(change)
- # c2 <- proc.time()[1]
- # cat ('sam',c2-c1,'\n')
+ # cc1 <- proc.time()[1]
+ # cat('samp',cc1-cc, '\n')
z$accepts[, i] <<- z$accept
z$parameters[i, , ] <<- z$thetaMat
z$MHaccepts[i, , , ] <<-
@@ -141,7 +142,6 @@
}
z$BayesAcceptances <<- rowSums(z$accepts)
z$nrunMH <<- storeNrunMH
-
}
##@sampleParameters algorithms propose new parameters and accept them or not
sampleParameters <- function(change=TRUE)
@@ -181,16 +181,13 @@
sigma=z$priorSigma[use, use])
}
)
- logpOld <- getProbs(z, z$chain)
+ logpOld <- z$loglik
- storeNrunMH <- z$nrunMH
- z$nrunMH <<- 0
thetaNew[, z$basicRate] <- exp(thetaNew[, z$basicRate])
z$thetaMat <<- thetaNew
- resp <- z$FRAN(z, returnChains=TRUE, byGroup=TRUE)
- logpNew <- getProbs(z, resp$chain)
-
+ logpNew <- getProbabilitiesFromC(z)[[1]]
proposalProbability <- priorNew - priorOld + logpNew - logpOld
+##cat(proposalProbability, priorNew, priorOld, logpNew, logpOld, '\n')
z$accept <<- log(runif(length(proposalProbability))) < proposalProbability
thetaOld[, z$basicRate] <- exp(thetaOld[, z$basicRate])
@@ -203,11 +200,7 @@
##print(z$thetaMat)
z$thetaMat[!z$accept, ] <<- thetaOld[!z$accept, ]
}
- z$nrunMH <<- storeNrunMH
- ## cat(thetaNew, priorNew, priorOld, logpNew, logpOld,
- ## exp(proposalProbability),
- ## z$accept,'\n')
-
+## print(thetaNew)
}
## ################################
## start of function proper
@@ -234,7 +227,7 @@
ctime <- proc.time()[1]
z <- initializeBayes(data, effects, model, nbrNodes, priorSigma,
- prevAns=prevAns)
+ prevAns=prevAns, clusterType=clusterType)
createStores()
z$sub <- 0
@@ -249,16 +242,15 @@
{
z$dfra <- dfra
}
- lambda <- z$theta[z$basicRate]
- dfra <- z$dfra
- z$dfra[z$basicRate, ] <- z$dfra[z$basicRate,] * lambda
- z$dfra[, z$basicRate] <- z$dfra[, z$basicRate] * lambda
- diag(z$dfra)[z$basicRate] <- lambda *
- diag(dfra)[z$basicRate] + lambda * lambda *
- colMeans(z$sf)[z$basicRate]
- # print(z$dfra)
- z$dfra <- solve(z$dfra)
- # print(z$dfra)
+ else
+ {
+ if (is.null(z$sf))
+ {
+ stop("need some scores to scale dfra")
+ }
+ z$dfra <- scaleDfra(z)
+ }
+
}
ctime1 <- proc.time()[1]
cat(ctime1-ctime,'\n')
@@ -284,7 +276,6 @@
{
MCMCcycle(nrunMH=4, nrunMHBatches=20)
}
-
print('endof warm')
ctime3<- proc.time()[1]
@@ -296,6 +287,7 @@
ctime4<- proc.time()[1]
cat('main', ii, ctime4-ctime3,'\n')
ctime3 <- ctime4
+
if (ii %% 10 == 0 && plotit) ## do some plots
{
cat('main after ii', ii, '\n')
@@ -329,7 +321,16 @@
"DelPerm", "InsMissing", "DelMissing",
"BayesAccepts")
varnames <- paste(names(thetadf)[-1], sep="", collapse= " + ")
- varcall <- paste("~ ", varnames, " | Group", sep="", collapse="")
+ if (z$nGroup > 1)
+ {
+ varcall <- paste("~ ", varnames, " | Group", sep="",
+ collapse="")
+ }
+ else
+ {
+ varcall <- paste("~ ", varnames, sep="", collapse="")
+
+ }
print(histogram(as.formula(varcall), data=thetadf, scales="free",
outer=TRUE, breaks=NULL, type="density",
panel=function(x, ...)
@@ -350,8 +351,18 @@
}
))
varnames <- paste(names(thetadf)[-1], sep="", collapse= " + ")
- varcall <- paste(varnames, "~ 1:", ii * nrunMHBatches * z$nGroup,
+ if (z$nGroup > 1)
+ {
+ varcall <- paste(varnames, "~ 1:", ii *
+ nrunMHBatches * z$nGroup,
" | Group", sep="", collapse="")
+ }
+ else
+ {
+ varcall <- paste(varnames, "~ 1:", ii *
+ nrunMHBatches * z$nGroup,
+ sep="", collapse="")
+ }
dev.set(tseriesplot)
print(xyplot(as.formula(varcall), data=thetadf, scales="free",
outer=TRUE))
@@ -380,7 +391,7 @@
##@initializeBayes algorithms do set up for Bayesian model
initializeBayes <- function(data, effects, model, nbrNodes, priorSigma,
- prevAns)
+ prevAns, clusterType=c("PSOCK", "FORK"))
{
## initialise
Report(openfiles=TRUE, type="n") #initialise with no file
@@ -417,24 +428,33 @@
WriteOutTheta(z)
- if (nbrNodes > 1)
+ if (nbrNodes > 1 && z$observations > 1)
{
require(parallel)
+ clusterType <- match.arg(clusterType)
+ if (clusterType == "PSOCK")
+ {
clusterString <- rep("localhost", nbrNodes)
z$cl <- makeCluster(clusterString, type = "PSOCK",
outfile = "cluster.out")
+ }
+ else
+ {
+ z$cl <- makeCluster(nbrNodes, type = "FORK",
+ outfile = "cluster.out")
+ }
clusterCall(z$cl, library, pkgname, character.only = TRUE)
clusterCall(z$cl, storeinFRANstore, FRANstore())
clusterCall(z$cl, FRANstore)
clusterCall(z$cl, initializeFRAN, z, model,
initC = TRUE, profileData=FALSE, returnDeps=FALSE)
- clusterSetRNGStream(z$cl,
- iseed = as.integer(runif(1, max=.Machine$integer.max)))
+ clusterSetRNGStream(z$cl, iseed = as.integer(runif(1,
+ max=.Machine$integer.max)))
}
z$scaleFactors <- rep(1, z$nGroup)
## z$returnDataFrame <- TRUE # chains come back as data frames not lists
- z$returnChains <- TRUE
+ z$returnChains <- FALSE
if (is.null(priorSigma))
{
z$priorSigma <- diag(z$pp) * 10000
@@ -451,9 +471,6 @@
use[!z$basicRate] <- TRUE
z$thetaMat[i, !use] <- NA
}
- z$callGrid <- cbind(rep(1:z$nGroup, z$groupPeriods),
- as.vector(unlist(sapply(z$groupPeriods,
- function(x) 1:x))))
z
}
##@getDFRA algorithms do a few ML iterations and calculate a derivative matrix
@@ -462,76 +479,32 @@
## do n MLmodelsteps with the initial thetas and get
## derivs
z$sdf <- vector("list", n)
- z$ssc <- matrix(0, nrow=n, ncol=z$pp)
+ z$sf <- matrix(0, nrow=n, ncol=z$pp)
z$Deriv <- TRUE
for (i in 1:n)
{
ans <- z$FRAN(z)
z$sdf[[i]] <- ans$dff
- z$ssc[i, ] <- colSums(ans$fra)
+ z$sf[i, ] <- colSums(ans$fra)
}
- dfra <- t(as.matrix(Reduce("+", z$sdf))) / length(z$sdf)
- ssc <- colMeans(z$ssc)
+ dfra <- t(as.matrix(Reduce("+", z$sdf) / length(z$sdf)))
z$dfra <- dfra
+ z$dfra <- scaleDfra(z)
+ z$Deriv <- FALSE
+ z
+}
+scaleDfra <- function(z)
+{
lambda <- z$theta[z$basicRate]
+ dfra <- z$dfra
z$dfra[z$basicRate, ] <- z$dfra[z$basicRate,] * lambda
z$dfra[, z$basicRate] <- z$dfra[, z$basicRate] * lambda
diag(z$dfra)[z$basicRate] <- lambda *
- diag(dfra)[z$basicRate] + lambda * lambda * ssc[z$basicRate]
- ##print(z$dfra)
- z$dfra <- solve(z$dfra)
- ##print(z$dfra)
- ##$eS <- eigen(z$dfra, symmetric=TRUE, EISPACK=TRUE)
- ##z$ev <- z$eS$values
- ##z$covFactor <- z$eS$vectors %*% diag(sqrt(pmax(z$ev, 0)), z$pp)
- z
+ diag(dfra)[z$basicRate] + lambda * lambda *
+ colMeans(z$sf)[z$basicRate]
+ chol2inv(chol(z$dfra))
}
-##@getLikelihood algorithms calculated likelihood for one chain
-getLikelihood <- function(chain, nactors, lambda, simpleRates)
-{
- loglik <- 0
- ncvals <- sapply(chain, function(x)x[[3]])
- nc <- nactors
- nc[] <- 0
- ncvals <- table(ncvals)
- nc[names(ncvals)] <- ncvals
- logChoiceProb <- sapply(chain, function(x)x[[9]])
- logOptionSetProb <- sapply(chain, function(x)x[[8]])
- loglik <- sum(logChoiceProb) + sum(logOptionSetProb)
- ##print(sum(logOptionSetProb))
- if (simpleRates)
- {
- loglik <- loglik - sum(nactors * lambda) + sum(nc * log(lambda))# -
- ## sum(lfactorial(nc)) don't need factorial in bayes!
- }
- else
- {
- mu <- attr(chain, "mu")
- sigma <- sqrt(attr(chain, "sigma2"))
- finalReciprocalRate <- attr(chain, "finalReciprocalRate")
- loglik <- loglik + dnorm(1, mu, sigma, log=TRUE) +
- log(finalReciprocalRate)
- }
- loglik
-}
-##@getProbs algorithms calculates likelihood sum over nested list of chains
-getProbs <- function(z, chain)
-{
- sapply(1: length(chain), function(i)
- {
- groupChain <- chain[[i]]
- sum(sapply(1:length(groupChain), function(j)
- {
- periodChain <- groupChain[[j]]
- theta1 <- z$thetaMat[i, ]
- k <- z$rateParameterPosition[[i]][[j]]
- getLikelihood(periodChain, z$nactors[[i]], theta1[k],
- z$simpleRates)
- }
- ))
- }
- )
-}
+
##@flattenChains algorithms converts a nested list of chains to a single list
flattenChains <- function(zz)
{
@@ -558,3 +531,56 @@
logdet <- sum(log(eigen(sigma, symmetric=TRUE, only.values=TRUE)$values))
-(ncol(x) * log(2 * pi) + logdet + distval) / 2
}
+
+##@getProbabilitiesFromC bayes gets loglik from chains in C
+getProbabilitiesFromC <- function(z, index=1, getScores=FALSE)
+{
+ ## expects maximum likelihood parallelisations
+ f <- FRANstore()
+
+ callGrid <- z$callGrid
+ ## z$int2 is the number of processors if iterating by period, so 1 means
+ ## we are not. Can only parallelize by period1
+ if (nrow(callGrid) == 1)
+ {
+ theta <- z$thetaMat[1,]
+ ans <- .Call("getChainProbabilities", PACKAGE = pkgname, f$pData,
+ f$pModel, as.integer(1), as.integer(1),
+ as.integer(index), f$myeffects, theta, getScores)
+ anss <- list(ans)
+ }
+ else
+ {
+ if (z$int2 == 1 )
+ {
+ anss <- apply(callGrid, 1,
+ doGetProbabilitiesFromC, z$thetaMat, index, getScores)
+ }
+ else
+ {
+ use <- 1:(min(nrow(callGrid), z$int2))
+ anss <- parRapply(z$cl[use], callGrid,
+ doGetProbabilitiesFromC, z$thetaMat, index,
+ getScores)
+ }
+ }
+ ans <- list()
+ ans[[1]] <- sum(sapply(anss, "[[", 1))
+ if (getScores)
+ {
+ ans[[2]] <- rowSums(sapply(anss, "[[", 2))
+ }
+ ans[[3]] <- sapply(anss, "[[", 3)
+ ans
+}
+
+##@doGetProbabilitiesFromC Maximum likelihood
+doGetProbabilitiesFromC <- function(x, thetaMat, index, getScores)
+{
+ f <- FRANstore()
+ theta <- thetaMat[x[1], ]
+ .Call("getChainProbabilities", PACKAGE = pkgname, f$pData,
+ f$pModel, as.integer(x[1]), as.integer(x[2]),
+ as.integer(index), f$myeffects, theta, getScores)
+
+}
Modified: pkg/RSiena/R/initializeFRAN.r
===================================================================
--- pkg/RSiena/R/initializeFRAN.r 2011-11-27 16:23:01 UTC (rev 185)
+++ pkg/RSiena/R/initializeFRAN.r 2011-12-04 17:51:35 UTC (rev 186)
@@ -9,7 +9,8 @@
# * and for ML making the initial chain.
# *****************************************************************************/
##@initializeFRAN siena07 reformat data and send to C. get targets.
-initializeFRAN <- function(z, x, data, effects, prevAns, initC, profileData,
+initializeFRAN <- function(z, x, data, effects, prevAns=NULL, initC,
+ profileData=FALSE,
returnDeps)
{
## fix up the interface so can call from outside robmon framework
@@ -25,7 +26,6 @@
{
z$int2 <- 1
}
-
if (!initC) ## ie first time round
{
if (!inherits(data,'siena'))
@@ -424,6 +424,12 @@
ans <- .Call("getTargets", PACKAGE=pkgname, pData, pModel, myeffects,
z$parallelTesting)
##stop('done')
+ ## create a grid of periods with group names in case want to
+ ## parallelize using this or to access chains easily
+ groupPeriods <- attr(f, "groupPeriods")
+ z$callGrid <- cbind(rep(1:nGroup, groupPeriods - 1),
+ as.vector(unlist(sapply(groupPeriods - 1,
+ function(x) 1:x))))
if (!x$maxlike)
{
z$targets <- rowSums(ans)
@@ -446,12 +452,6 @@
round(z$nrunMH / 100) * 100, z$nrunMH)
##thetaMat is to allow different thetas for each group in Bayes
z$thetaMat <- matrix(z$theta, nrow=nGroup, ncol=z$pp, byrow=TRUE)
- ## create a grid of periods with group names in case want to
- ## parallelize using this
- groupPeriods <- attr(f, "groupPeriods")
- z$callGrid <- cbind(rep(1:nGroup, groupPeriods - 1),
- as.vector(unlist(sapply(groupPeriods - 1,
- function(x) 1:x))))
}
}
@@ -479,19 +479,20 @@
CONDTARGET <- NULL
}
- ans <- .Call("setupModelOptions", PACKAGE=pkgname,
- pData, pModel, MAXDEGREE, CONDVAR, CONDTARGET,
- profileData, z$parallelTesting, x$modelType)
-
- if (x$maxlike)
- {
simpleRates <- TRUE
if (any(!z$effects$basicRate & z$effects$type =="rate"))
{
- # browser()
+ ## browser()
simpleRates <- FALSE
}
z$simpleRates <- simpleRates
+
+ ans <- .Call("setupModelOptions", PACKAGE=pkgname,
+ pData, pModel, MAXDEGREE, CONDVAR, CONDTARGET,
+ profileData, z$parallelTesting, x$modelType, z$simpleRates)
+
+ if (x$maxlike)
+ {
if (!initC)
{
## set up chains and do initial steps
@@ -523,17 +524,15 @@
x$prdrms)
##cat(z$probs,'\n')
ans <- .Call("mlMakeChains", PACKAGE=pkgname, pData, pModel,
- simpleRates, z$probs, z$prmin, z$prmib,
+ z$probs, z$prmin, z$prmib,
x$minimumPermutationLength,
x$maximumPermutationLength,
x$initialPermutationLength)
f$minimalChain <- ans[[1]]
f$chain <- ans[[2]]
- f$simpleRates <- simpleRates
## print(nrow(ans[[1]][[1]]))
## print(nrow(ans[[2]][[1]]))
- ## browser()
}
else ## set up the initial chains in the sub processes
{
@@ -542,11 +541,11 @@
simpleRates, z$probs, z$prmin, z$prmib,
x$minimumPermutationLength,
x$maximumPermutationLength,
- x$initialPermutationLength, ff$chain, ff$missingChain)
+ x$initialPermutationLength, ff$chain)
f$chain <- ff$chain
- f$missingChain <- ff$missingChain
}
}
+ f$simpleRates <- simpleRates
f$myeffects <- myeffects
f$myCompleteEffects <- myCompleteEffects
if (!initC)
Modified: pkg/RSiena/R/maxlikec.r
===================================================================
--- pkg/RSiena/R/maxlikec.r 2011-11-27 16:23:01 UTC (rev 185)
+++ pkg/RSiena/R/maxlikec.r 2011-12-04 17:51:35 UTC (rev 186)
@@ -12,10 +12,10 @@
##@maxlikec siena07 ML Simulation Module
maxlikec <- function(z, x, INIT=FALSE, TERM=FALSE, initC=FALSE, data=NULL,
effects=NULL, profileData=FALSE, prevAns=NULL,
- returnChains=FALSE, byGroup=FALSE,
- returnDataFrame=FALSE, byWave=FALSE)
+ returnChains=FALSE, byGroup=FALSE, returnDataFrame=FALSE,
+ byWave=FALSE, returnLoglik=FALSE, onlyLoglik=FALSE)
{
- if (INIT || initC) ## initC is to initialise multiple C processes in phase3
+ if (INIT || initC) ## initC is to initialise multiple C processes in
{
z <- initializeFRAN(z, x, data, effects, prevAns, initC,
profileData=profileData, returnDeps=FALSE)
@@ -66,7 +66,9 @@
f$pModel, f$myeffects, theta,
1, 1, z$nrunMH, z$addChainToStore,
z$needChangeContributions, z$returnDataFrame,
- z$returnChains)
+ z$returnChains, returnLoglik, onlyLoglik)
+ if (!onlyLoglik)
+ {
ans[[6]] <- list(ans[[6]])
ans[[7]] <- list(ans[[7]])
if (byGroup)
@@ -78,13 +80,22 @@
}
else
{
+ ans[[2]] <- list(ans[[2]])
+ ans[[3]] <- list(ans[[3]])
+ ans[[4]] <- list(ans[[4]])
+
+ }
+ }
+ else
+ {
if (z$int2 == 1)
{
anss <- apply(cbind(callGrid, 1:nrow(callGrid)),
1, doMLModel, z$Deriv, z$thetaMat,
z$nrunMH, z$addChainToStore,
z$needChangeContributions, z$returnDataFrame,
- z$returnChains, byGroup, z$theta)
+ z$returnChains, byGroup, z$theta, returnLoglik,
+ onlyLoglik)
}
else
{
@@ -94,35 +105,19 @@
z$nrunMH, z$addChainToStore,
z$needChangeContributions,
z$returnDataFrame, z$returnChains, byGroup,
- z$theta)
+ z$theta, returnLoglik, onlyLoglik)
}
## reorganize the anss so it looks like the normal one
- ans <- NULL
+ ans <- list()
+ if (!onlyLoglik)
+ {
ans[[1]] <- sapply(anss, "[[", 1) ## statistics
ans[[2]] <- NULL ## scores
ans[[3]] <- NULL ## seeds
ans[[4]] <- NULL ## ntim
ans[[5]] <- NULL # randomseed
- ## if (returnDeps) ## this is for dependent variables but
- ## ## these are not returned yet
- ## {
- ## fff <- lapply(anss, "[[", 6)
- ## fff <- split(fff, callGrid[1, ]) ## split by group
- ## ans[[6]] <-
- ## lapply(fff, function(x)
- ## {
- ## lapply(1:length(f$depNames), function(x, z)
- ## lapply(z, "[[", x), z=x)
- ## }
- ## )
- ## }
- ## else
- ## {
- ## ans[[6]] <- NULL
- ## }
if (z$returnChains)
{
- ##TODO put names on these?
fff <- lapply(anss, function(x) x[[6]][[1]])
fff <- split(fff, callGrid[, 1 ]) ## split by group
ans[[6]] <- fff
@@ -137,38 +132,20 @@
ans[[9]] <- Reduce("+", ans[[9]]) ## rejects
ans[[10]] <- Reduce("+", ans[[10]]) ## aborts
}
+ ans[[11]] <- sapply(anss, "[[", 11)
+ }
+ else ##onlyLoglik is always byGroup (bayes)
+ {
+ ans[[1]] <- sum(sapply(anss, '[[', 1))## loglik
+ ans[[2]] <- lapply(anss, "[[", 2)
+ ans[[3]] <- lapply(anss, "[[", 3)
+ ans[[4]] <- lapply(anss, "[[", 4)
+ }
}
- fra <- -t(ans[[1]]) ##note sign change
-
FRANstore(f)
- #if (returnDeps || z$returnChains)
- #{
- # sims <- ans[[6]]
- #}
- #else
- #{
- sims <- NULL
- #}
- ##print(length(sims[[1]]))
- ##if (returnDeps) ## this is not finished yet!
- ##{
- ## ## attach the names
- ## names(sims) <- f$groupNames
- ## periodNo <- 1
- ## for (i in 1:length(sims))
- ## {
- ## names(sims[[i]]) <- f$depNames
- ## for (j in 1:length(sims[[i]]))
- ## {
- ## periodNos <- periodNo:(periodNo + length(sims[[i]][[j]]) - 1)
- ## names(sims[[i]][[j]]) <- periodNos
- ## }
- ## periodNo <- periodNos[length(periodNos)] + 2
- ## }
- ##}
- if (z$Deriv) ## need to reformat the derivatives
+ if (z$Deriv && !onlyLoglik) ## need to reformat the derivatives
{
resp <- reformatDerivs(z, f, ans[[7]])
dff <- resp[[1]]
@@ -179,18 +156,26 @@
dff <- NULL
dff2 <- NULL
}
- ## browser()
- ##print(length(ans[[6]][[1]][[1]]))
+ if (!onlyLoglik)
+ {
+ fra <- -t(ans[[1]]) ##note sign change
+
list(fra = fra, ntim0 = NULL, feasible = TRUE, OK = TRUE,
- sims=sims, dff = dff, dff2=dff2,
+ sims=NULL, dff = dff, dff2=dff2,
chain = ans[[6]], accepts=ans[[8]],
- rejects= ans[[9]], aborts=ans[[10]])
+ rejects= ans[[9]], aborts=ans[[10]], loglik=ans[[11]])
+ }
+ else
+ {
+ list(loglik=ans[[1]], accepts=ans[[2]],
+ rejects= ans[[3]], aborts=ans[[4]])
}
+}
##@doMLModel Maximum likelihood
doMLModel <- function(x, Deriv, thetaMat, nrunMH, addChainToStore,
needChangeContributions, returnDataFrame, returnChains,
- byGroup, theta)
+ byGroup, theta, returnLoglik, onlyLoglik)
{
f <- FRANstore()
if (byGroup)
@@ -204,9 +189,11 @@
.Call("mlPeriod", PACKAGE=pkgname, Deriv, f$pData,
f$pModel, f$myeffects, theta,
as.integer(x[1]), as.integer(x[2]), nrunMH[x[3]], addChainToStore,
- needChangeContributions, returnDataFrame, returnChains)
+ needChangeContributions, returnDataFrame, returnChains,
+ returnLoglik, onlyLoglik)
}
+##@reformatDerivs Maximum likelihood move this back to inline or internal function
reformatDerivs <- function(z, f, derivList)
{
## current format is a list of vectors of the lower? triangle
Modified: pkg/RSiena/R/phase1.r
===================================================================
--- pkg/RSiena/R/phase1.r 2011-11-27 16:23:01 UTC (rev 185)
+++ pkg/RSiena/R/phase1.r 2011-12-04 17:51:35 UTC (rev 186)
@@ -351,7 +351,7 @@
}
else
{
- zz <- clusterCall(z$cl, usesim, zdummy, xsmall,
+ zz <- clusterCall(z$cl, simstats0c, zdummy, xsmall,
INIT=FALSE, fromFiniteDiff=TRUE)
}
if (int == 1)
Modified: pkg/RSiena/R/phase2.r
===================================================================
--- pkg/RSiena/R/phase2.r 2011-11-27 16:23:01 UTC (rev 185)
+++ pkg/RSiena/R/phase2.r 2011-12-04 17:51:35 UTC (rev 186)
@@ -12,11 +12,6 @@
## args: z: internal control object
## x: model object (readonly as not returned)
-##@usesim siena07 Used to avoid Namespace problems with multiple processes
-usesim <- function(...)
-{
- simstats0c(...)
-}
##@storeinFRANstore siena07 Used to avoid Namespace problems with multiple processes
storeinFRANstore <- function(...)
{
@@ -258,7 +253,7 @@
}
}
zsmall$nit <- z$nit
- if (z$int == 1)
+ if (z$int == 1) ## not parallel runs at this level
{
zz <- x$FRAN(zsmall, xsmall)
## browser()
@@ -271,7 +266,7 @@
}
else
{
- zz <- clusterCall(z$cl, usesim, zsmall, xsmall)
+ zz <- clusterCall(z$cl, simstats0c, zsmall, xsmall)
fra <- sapply(zz, function(x) colSums(x$fra)- z$targets)
dim(fra) <- c(z$pp, z$int)
fra <- rowMeans(fra)
@@ -316,7 +311,9 @@
z$truncated[z$nit] <- TRUE
}
else
+ {
maxrat <- 1.0
+ }
fchange<- z$gain * fra * maxrat / diag(z$dfra)
}
else
@@ -450,7 +447,7 @@
}
else
{
- zz <- clusterCall(z$cl, usesim, zsmall, xsmall)
+ zz <- clusterCall(z$cl, simstats0c, zsmall, xsmall)
fra <- sapply(zz, function(x) colSums(x$fra) - z$targets)
dim(fra) <- c(z$pp, z$int)
fra <- rowMeans(fra)
@@ -538,9 +535,9 @@
{
nGroup <- z$f$nGroup
z$chains <- lapply(1:nGroup, function(x)
- lapply(1:z$groupPeriods[x], function(y)
+ lapply(1:(z$f$groupPeriods[x] - 1), function(y)
vector("list", niter)))
- z$lik0 <- rep(0, niter * sum(z$groupPeriods - 1))
+ z$lik0 <- rep(0, niter * sum(z$f$groupPeriods - 1))
z$iterFra <- matrix(0, ncol=z$pp, nrow=niter)
z$thetaStore <- matrix(0, ncol=z$pp, nrow=z$n2max)
z$sf <- matrix(0, ncol=z$pp, nrow=z$n2max)
@@ -552,7 +549,7 @@
{
for (ii in 1:z$nGroup)
{
- for (jj in 1:z$groupPeriods[ii])
+ for (jj in 1:(z$f$groupPeriods[ii] - 1))
{
z$chains[[ii]][[jj]][[i]] <- zz$chain[[ii]][[jj]]
}
@@ -569,7 +566,7 @@
storeSub <- 1
for (ii in 1:z$nGroup)
{
- for (jj in 1:z$groupPeriods[ii])
+ for (jj in 1:(z$f$groupPeriods[ii] - 1))
{
chains <- z$chains[[ii]][[jj]]
for (chain in chains)
@@ -747,9 +744,9 @@
return(z)
}
nGroup <- z$f$nGroup
+ groupPeriods <- attr(z$f, "groupPeriods")
z$nDependentVariables <- length(z$f$depNames)
- atts <- attributes(z$f)
- z$groupPeriods <- atts$groupPeriods - 1
+ ## atts <- attributes(z$f)
netnames <- names(z$f[[1]]$depvars)
z$nactors <- lapply(1:nGroup, function(i, periods, data)
{
@@ -757,7 +754,7 @@
dim(x)[1])
tmp <- tmp[match(netnames, names(data[[i]]$depvars))]
tmp
- }, periods=z$groupPeriods, data=z$f[1:nGroup]
+ }, periods=groupPeriods - 1, data=z$f[1:nGroup]
)
z$rateParameterPosition <-
lapply(1:nGroup, function(i, periods, data)
@@ -776,7 +773,7 @@
tmp
}
)
- }, periods=z$groupPeriods, data=z$f[1:nGroup]
+ }, periods=groupPeriods - 1, data=z$f[1:nGroup]
)
z$evalParameterPosition <-
lapply(netnames, function(x)
Modified: pkg/RSiena/R/phase3.r
===================================================================
--- pkg/RSiena/R/phase3.r 2011-11-27 16:23:01 UTC (rev 185)
+++ pkg/RSiena/R/phase3.r 2011-12-04 17:51:35 UTC (rev 186)
@@ -480,9 +480,9 @@
z$phase1Its <- z$phase1Its + 1
}
}
- else
+ else ### only do parallels at this level for forward simulation
{
- zz <- clusterCall(z$cl, usesim, zsmall, xsmall)
+ zz <- clusterCall(z$cl, simstats0c, zsmall, xsmall)
z$n <- z$n + z$int
if (phase == 1)
{
Modified: pkg/RSiena/R/robmon.r
===================================================================
--- pkg/RSiena/R/robmon.r 2011-11-27 16:23:01 UTC (rev 185)
+++ pkg/RSiena/R/robmon.r 2011-12-04 17:51:35 UTC (rev 186)
@@ -68,6 +68,7 @@
outfile = "cluster.out")
}
clusterCall(cl, library, pkgname, character.only = TRUE)
+ ##parLapply(cl, c('f1','f2'), sink)
z$oldRandomNumbers <- .Random.seed
if (R.version$minor < 14.0) ## fake this to recreate old results
## if (TRUE)
@@ -83,7 +84,14 @@
clusterCall(cl, storeinFRANstore, FRANstore())
if (initC)
{
- clusterCall(cl, usesim, z, x, INIT=FALSE, initC = initC)
+ if (!x$maxlike)
+ {
+ clusterCall(cl, simstats0c, z, x, INIT=FALSE, initC = initC)
+ }
+ else
+ {
+ clusterCall(cl, maxlikec, z, x, INIT=FALSE, initC = initC)
+ }
}
z$cl <- cl
}
Modified: pkg/RSiena/R/simstatsc.r
===================================================================
--- pkg/RSiena/R/simstatsc.r 2011-11-27 16:23:01 UTC (rev 185)
+++ pkg/RSiena/R/simstatsc.r 2011-12-04 17:51:35 UTC (rev 186)
@@ -12,13 +12,13 @@
simstats0c <- function(z, x, INIT=FALSE, TERM=FALSE, initC=FALSE, data=NULL,
effects=NULL, fromFiniteDiff=FALSE,
profileData=FALSE, prevAns=NULL, returnDeps=FALSE,
- returnChains=FALSE, byWave=FALSE)
+ returnChains=FALSE, byWave=FALSE, returnLoglik=FALSE)
{
if (INIT || initC) ## initC is to initialise multiple C processes in phase3
{
z <- initializeFRAN(z, x, data, effects, prevAns, initC,
profileData=profileData, returnDeps=returnDeps)
-
+ z$returnChains <- returnChains
z$byWave <- byWave
if (initC)
{
@@ -82,12 +82,7 @@
}
## cat(randomseed2, '\n')
}
- ## create a grid of periods with group names in case want to parallelize
- ## using this
- groupPeriods <- attr(f, "groupPeriods")
- callGrid <- cbind(rep(1:f$nGroup, groupPeriods - 1),
- as.vector(unlist(sapply(groupPeriods - 1,
- function(x) 1:x))))
+ ##callGrid <- z$callGrid
if (R.version$minor < 14.0) ##fake this to repeat old results
## if (TRUE)
{
@@ -98,41 +93,41 @@
useStreams <- FALSE
}
## z$int2 is the number of processors if iterating by period, so 1 means
- ## we are not
- if (z$int2==1 || nrow(callGrid) == 1)
+ ## we are not. Now we have removed option to parallelize by period
+ # if (z$int2==1 || nrow(callGrid) == 1)
{
## cat("theta", z$theta, "\n")
ans <- .Call('model', PACKAGE=pkgname, z$Deriv, f$pData, seeds,
fromFiniteDiff, f$pModel, f$myeffects, z$theta,
randomseed2, returnDeps, z$FinDiff.method,
!is.null(z$cl) && useStreams, z$addChainToStore,
- z$needChangeContributions, returnChains)
+ z$needChangeContributions, returnChains, returnLoglik)
}
- else
- {
- use <- 1:(min(nrow(callGrid), z$int2))
- anss <- parRapply(z$cl[use], callGrid, doModel,
- z$Deriv, seeds, fromFiniteDiff, z$theta,
- randomseed2, returnDeps, z$FinDiff.method, useStreams,
- returnChains)
- ## reorganize the anss so it looks like the normal one
- ## browser()
- ans <- NULL
- ans[[1]] <- sapply(anss, "[[", 1) ## statistics
- ans[[2]] <- sapply(anss, "[[", 2) ## scores
- ans[[3]] <- split(lapply(anss, "[[", 3), callGrid[, 1]) ## seeds
- ans[[4]] <- sapply(anss, "[[", 4) # ntim
- ans[[5]] <- NULL # randomseed not sensible here
- fff <- lapply(anss, "[[", 6)
- fff <- split(fff, callGrid[, 1])
- ans[[6]] <-
- lapply(fff, function(x)
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/rsiena -r 186
More information about the Rsiena-commits
mailing list