[Rsiena-commits] r178 - in pkg: RSiena RSiena/R RSiena/inst RSiena/inst/doc RSiena/inst/examples RSiena/inst/scripts RSiena/man RSiena/src RSiena/src/model/ml RSiena/src/model/variables RSiena/tests RSienaTest RSienaTest/R RSienaTest/doc RSienaTest/inst/doc RSienaTest/man RSienaTest/src/model/ml RSienaTest/tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Oct 27 19:30:41 CEST 2011
Author: ripleyrm
Date: 2011-10-27 19:30:41 +0200 (Thu, 27 Oct 2011)
New Revision: 178
Added:
pkg/RSiena/inst/examples/s50e.dat
pkg/RSiena/inst/scripts/
pkg/RSiena/inst/scripts/RSienaSNADescriptives.R
pkg/RSiena/inst/scripts/Rscript01DataFormat.R
pkg/RSiena/inst/scripts/Rscript02SienaVariableFormat.R
pkg/RSiena/inst/scripts/Rscript03SienaRunModel.R
pkg/RSiena/inst/scripts/Rscript04SienaBehaviour.R
pkg/RSiena/man/bayes.Rd
pkg/RSiena/man/maxlikec.Rd
pkg/RSiena/man/updateTheta.Rd
pkg/RSiena/man/xtable.Rd
pkg/RSiena/tests/scriptfile.Rout.save
pkg/RSiena/tests/scriptfile.Rout.win
pkg/RSiena/tests/scripts.R
pkg/RSiena/tests/scripts.Rout.save
pkg/RSienaTest/doc/_emacs
Modified:
pkg/RSiena/DESCRIPTION
pkg/RSiena/NAMESPACE
pkg/RSiena/R/bayes.r
pkg/RSiena/R/effects.r
pkg/RSiena/R/effectsMethods.r
pkg/RSiena/R/initializeFRAN.r
pkg/RSiena/R/maxlikec.r
pkg/RSiena/R/phase1.r
pkg/RSiena/R/phase3.r
pkg/RSiena/R/robmon.r
pkg/RSiena/R/siena07.r
pkg/RSiena/R/sienaDataCreate.r
pkg/RSiena/R/sienaModelCreate.r
pkg/RSiena/R/sienaTimeTest.r
pkg/RSiena/R/sienaprint.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/siena07.Rd
pkg/RSiena/man/sienaModelCreate.Rd
pkg/RSiena/man/sienaTimeTest.Rd
pkg/RSiena/man/simstats0c.Rd
pkg/RSiena/src/model/ml/Chain.cpp
pkg/RSiena/src/model/ml/MLSimulation.cpp
pkg/RSiena/src/model/ml/MLSimulation.h
pkg/RSiena/src/model/ml/MiniStep.cpp
pkg/RSiena/src/model/variables/DependentVariable.cpp
pkg/RSiena/src/model/variables/DependentVariable.h
pkg/RSiena/src/model/variables/NetworkVariable.cpp
pkg/RSiena/src/siena07models.cpp
pkg/RSiena/src/siena07models.h
pkg/RSiena/src/siena07setup.cpp
pkg/RSiena/tests/parallel.Rout.save
pkg/RSienaTest/DESCRIPTION
pkg/RSienaTest/R/phase3.r
pkg/RSienaTest/R/sienaGOF.r
pkg/RSienaTest/R/sienaTimeTest.r
pkg/RSienaTest/changeLog
pkg/RSienaTest/doc/RSienaDeveloper.tex
pkg/RSienaTest/doc/RSiena_Manual.tex
pkg/RSienaTest/inst/doc/RSiena_Manual.pdf
pkg/RSienaTest/man/RSiena-package.Rd
pkg/RSienaTest/src/model/ml/MLSimulation.cpp
pkg/RSienaTest/tests/scriptfile.Rout.win
pkg/RSienaTest/tests/scripts.Rout.save
Log:
New parallel package for R 2.14.0. RSiena updated with most of recent changes to RSienaTest.
Modified: pkg/RSiena/DESCRIPTION
===================================================================
--- pkg/RSiena/DESCRIPTION 2011-10-27 13:28:48 UTC (rev 177)
+++ pkg/RSiena/DESCRIPTION 2011-10-27 17:30:41 UTC (rev 178)
@@ -1,12 +1,13 @@
Package: RSiena
Type: Package
Title: Siena - Simulation Investigation for Empirical Network Analysis
-Version: 1.0.12.174
-Date: 2011-10-04
+Version: 1.0.12.178
+Date: 2011-10-27
Author: Various
Depends: R (>= 2.10.0)
Imports: Matrix
-Suggests: tcltk, snow, rlecuyer, network, codetools, lattice, MASS, xtable
+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>
Description: Fits models to longitudinal networks
Modified: pkg/RSiena/NAMESPACE
===================================================================
--- pkg/RSiena/NAMESPACE 2011-10-27 13:28:48 UTC (rev 177)
+++ pkg/RSiena/NAMESPACE 2011-10-27 17:30:41 UTC (rev 178)
@@ -1,11 +1,11 @@
useDynLib(RSiena)
export(coCovar, coDyadCovar, getEffects, model.create, print01Report,
-siena01Gui, siena07, sienaCompositionChange,
+siena01Gui, siena07, sienaCompositionChange, bayes, updateTheta,
sienaCompositionChangeFromFile, sienaDataCreate, sienaDataCreateFromSession,
sienaGroupCreate, sienaModelCreate, sienaNet, sienaNodeSet, xtable.sienaFit,
varCovar, varDyadCovar, setEffect, includeEffects, includeInteraction,
effectsDocumentation, sienaDataConstraint, print.xtable.sienaFit,
- installGui, siena08, iwlsm, sienaTimeTest, includeTimeDummy)
+ installGui, siena08, iwlsm, sienaTimeTest, includeTimeDummy, xtable)
import(Matrix)
Modified: pkg/RSiena/R/bayes.r
===================================================================
--- pkg/RSiena/R/bayes.r 2011-10-27 13:28:48 UTC (rev 177)
+++ pkg/RSiena/R/bayes.r 2011-10-27 17:30:41 UTC (rev 178)
@@ -10,7 +10,7 @@
## ****************************************************************************/
##@bayes Bayesian fit a Bayesian model
bayes <- function(data, effects, model, nwarm=100, nmain=100, nrunMHBatches=20,
- nrunMH=100, save=TRUE, nbrNodes=1, seed=1, dfra=NULL,
+ plotit=FALSE, nbrNodes=1, dfra=NULL, n=10,
priorSigma=NULL, prevAns=NULL, getDocumentation=FALSE)
{
##@createStores internal bayes Bayesian set up stores
@@ -22,9 +22,9 @@
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$MHacceptances <- array(NA, dim=c(numberRows, z$nGroup, 7))
- z$MHrejections <- array(NA, dim=c(numberRows, z$nGroup, 7))
- z$MHproportions <- array(NA, dim=c(numberRows, z$nGroup, 7))
+ z$MHacceptances <- array(NA, dim=c(numberRows, z$nGroup, 9))
+ z$MHrejections <- array(NA, dim=c(numberRows, z$nGroup, 9))
+ z$MHproportions <- array(NA, dim=c(numberRows, z$nGroup, 9))
z
}
##@storeData internal bayes Bayesian put data in stores
@@ -121,7 +121,7 @@
}
}
- z <- initializeBayes(data, effects, model, nbrNodes, seed, priorSigma,
+ z <- initializeBayes(data, effects, model, nbrNodes, priorSigma,
prevAns=prevAns)
z <- createStores(z)
@@ -129,7 +129,7 @@
if (is.null(z$dfra) && is.null(dfra))
{
- z <- getDFRA(z, 10)
+ z <- getDFRA(z, n)
}
else
{
@@ -152,7 +152,7 @@
}
z <- improveMH(z)
- if (!save)
+ if (plotit)
{
require(lattice)
dev.new()
@@ -172,10 +172,10 @@
for (ii in 1:nmain)
{
- z <- MCMCcycle(z, nrunMH=nrunMH, nrunMHBatches=nrunMHBatches)
+ z <- MCMCcycle(z, nrunMH=z$nrunMH, nrunMHBatches=nrunMHBatches)
z <- storeData()
- if (ii %% 10 == 0 && !save) ## do some plots
+ if (ii %% 10 == 0 && plotit) ## do some plots
{
cat('main after ii', ii, '\n')
dev.set(thetaplot)
@@ -261,9 +261,9 @@
{
z$accepts <- matrix(NA, nrow=z$nGroup, nrunMHBatches)
z$parameters <- array(NA, dim=c(nrunMHBatches, z$nGroup, z$pp))
- z$MHaccepts <- array(NA, dim=c(nrunMHBatches, z$nGroup, 7))
- z$MHrejects <- array(NA, dim=c(nrunMHBatches, z$nGroup, 7))
- z$MHaborts <- array(NA, dim=c(nrunMHBatches, z$nGroup, 7))
+ z$MHaccepts <- array(NA, dim=c(nrunMHBatches, z$nGroup, 9))
+ z$MHrejects <- array(NA, dim=c(nrunMHBatches, z$nGroup, 9))
+ z$MHaborts <- array(NA, dim=c(nrunMHBatches, z$nGroup, 9))
storeNrunMH <- z$nrunMH
z$nrunMH <- nrunMH
for (i in 1:nrunMHBatches)
@@ -357,11 +357,10 @@
}
##@initializeBayes algorithms do set up for Bayesian model
-initializeBayes <- function(data, effects, model, nbrNodes, seed, priorSigma,
+initializeBayes <- function(data, effects, model, nbrNodes, priorSigma,
prevAns)
{
## initialise
- set.seed(seed)
Report(openfiles=TRUE, type="n") #initialise with no file
z <- NULL
z$Phase <- 1
@@ -377,8 +376,19 @@
if (!is.null(model$randomSeed))
{
set.seed(model$randomSeed)
+ ##seed <- model$randomSeed
}
- z$FRAN <- getFromNamespace(model$FRANname, pos=grep("RSiena",
+ else
+ {
+ if (exists(".Random.seed"))
+ {
+ rm(.Random.seed, pos=1)
+ }
+ newseed <- trunc(runif(1) * 1000000)
+ set.seed(newseed) ## get R to create a random number seed for me.
+ ##seed <- NULL
+ }
+ z$FRAN <- getFromNamespace(model$FRANname, pos=grep("RSiena",
search())[1])
z <- z$FRAN(z, model, INIT=TRUE, data=data, effects=effects,
prevAns=prevAns)
Modified: pkg/RSiena/R/effects.r
===================================================================
--- pkg/RSiena/R/effects.r 2011-10-27 13:28:48 UTC (rev 177)
+++ pkg/RSiena/R/effects.r 2011-10-27 17:30:41 UTC (rev 178)
@@ -209,7 +209,6 @@
}
}
-### not sure we need this: if so then check relevant combinations of nodesets
if (length(xx$cCovars) + length(xx$vCovars) +
length(xx$dycCovars) + length(xx$dyvCovars) +
length(types=='behavior') > 0)
@@ -606,7 +605,6 @@
rateEffects <- rbind(rateEffects, tmp$rateEff)
}
}
-### not sure we need this: if so then check relevant combinations of nodesets
if (length(xx$cCovars) + length(xx$vCovars) +
length(xx$dycCovars) + length(xx$dyvCovars) +
length(types=='behavior') > 0)
@@ -665,9 +663,6 @@
objEffects$effectName <- paste(varname, ': ',
objEffects$effectName, sep = '')
}
- ## now create the real effects, extra rows for endowment effects etc
- ##objEffects <- createObjEffectList(objEffects, varname)
- ##rateEffects <- createRateEffectList(rateEffects, varname)
objEffects$functionName[objEffects$type == 'endow'] <-
paste('Lost ties:',
@@ -695,12 +690,6 @@
objEffects <-
objEffects[!objEffects$shortName == "density", ]
}
- ## if (attr(xx$depvars[[i]],'uponly') || attr(xx$depvars[[i]],
- ## 'downonly'))
- ## {
- ## objEffects <-
- ## objEffects[!objEffects$shortName == "density", ]
- ## }
rateEffects$basicRate[1:observations] <- TRUE
@@ -1418,4 +1407,3 @@
list(startRate=startRate, degree=alphaf1, alpha=alpha, prec=prec, tmp=tmp,
untrimmed = untrimmed)
}
-
Modified: pkg/RSiena/R/effectsMethods.r
===================================================================
--- pkg/RSiena/R/effectsMethods.r 2011-10-27 13:28:48 UTC (rev 177)
+++ pkg/RSiena/R/effectsMethods.r 2011-10-27 17:30:41 UTC (rev 178)
@@ -150,6 +150,10 @@
edit.sienaEffects <- function(name, ...)
{
+ if (!interactive())
+ {
+ return(name)
+ }
## store the original column order
originalNames <- names(name)
## move function name and other things out of the way
Modified: pkg/RSiena/R/initializeFRAN.r
===================================================================
--- pkg/RSiena/R/initializeFRAN.r 2011-10-27 13:28:48 UTC (rev 177)
+++ pkg/RSiena/R/initializeFRAN.r 2011-10-27 17:30:41 UTC (rev 178)
@@ -69,35 +69,6 @@
}
effects$initialValue <- defaultEffects$initialValue
}
- else
- {
- if (!is.null(prevAns) && inherits(prevAns, "sienaFit"))
- {
- prevEffects <- prevAns$requestedEffects
- prevEffects$initialValue <- prevAns$theta
- if (prevAns$cconditional)
- {
- condEffects <- attr(prevAns$f, "condEffects")
- condEffects$initialValue <- prevAns$rate
- prevEffects <- rbind(prevEffects, condEffects)
- }
- oldlist <- apply(prevEffects, 1, function(x)
- paste(x[c("name", "shortName",
- "type", "groupName",
- "interaction1", "interaction2",
- "period")],
- collapse="|"))
- efflist <- apply(effects, 1, function(x)
- paste(x[c("name", "shortName",
- "type", "groupName",
- "interaction1", "interaction2",
- "period")],
- collapse="|"))
- use <- efflist %in% oldlist
- effects$initialValue[use] <-
- prevEffects$initialValue[match(efflist, oldlist)][use]
- }
- }
## get data object into group format to save coping with two
## different formats
if (inherits(data, 'sienaGroup'))
@@ -109,10 +80,17 @@
nGroup <- 1
data <- sienaGroupCreate(list(data), singleOK=TRUE)
}
- ## add any effects needed for time dummies
+ ## add any effects needed for time dummies
tmp <- sienaTimeFix(effects, data)
data <- tmp$data
effects <- tmp$effects
+ if (!x$useStdInits)
+ {
+ if (!is.null(prevAns) && inherits(prevAns, "sienaFit"))
+ {
+ effects <- updateTheta(effects, prevAns)
+ }
+ }
## find any effects not included which are needed for interactions
tmpEffects <- effects[effects$include, ]
interactionNos <- unique(c(tmpEffects$effect1, tmpEffects$effect2,
@@ -241,7 +219,7 @@
attr(data, "numberNonMissingBehavior")
attr(f, "numberMissingBehavior") <- attr(data, "numberMissingBehavior")
- # attr(f, "totalMissings") <- attr(data, "totalMissings")
+ ## attr(f, "totalMissings") <- attr(data, "totalMissings")
if (x$maxlike && x$FinDiff.method)
{
@@ -310,7 +288,7 @@
pData <- .Call('setupData', PACKAGE=pkgname,
lapply(f, function(x)(as.integer(x$observations))),
lapply(f, function(x)(x$nodeSets)))
- ans <- .Call('OneMode', PACKAGE=pkgname,
+ ans <- .Call('OneMode', PACKAGE=pkgname,
pData, lapply(f, function(x)x$nets))
ans <- .Call('Bipartite', PACKAGE=pkgname,
pData, lapply(f, function(x)x$bipartites))
@@ -320,11 +298,11 @@
pData, lapply(f, function(x)x$cCovars))
ans <-.Call('ChangingCovariates', PACKAGE=pkgname,
pData, lapply(f, function(x)x$vCovars))
- ans <-.Call('DyadicCovariates', PACKAGE=pkgname,
+ ans <-.Call('DyadicCovariates', PACKAGE=pkgname,
pData, lapply(f, function(x)x$dycCovars))
- ans <-.Call('ChangingDyadicCovariates', PACKAGE=pkgname,
+ ans <-.Call('ChangingDyadicCovariates', PACKAGE=pkgname,
pData, lapply(f, function(x)x$dyvCovars))
- ans <-.Call('ExogEvent', PACKAGE=pkgname,
+ ans <-.Call('ExogEvent', PACKAGE=pkgname,
pData, lapply(f, function(x)x$exog))
## split the names of the constraints
higher <- attr(f, "allHigher")
@@ -396,7 +374,7 @@
interactionEffectsl <- ff$interactionEffectsl
types <- ff$types
}
- ans <- .Call('effects', PACKAGE=pkgname, pData, basicEffects)
+ ans <- .Call('effects', PACKAGE=pkgname, pData, basicEffects)
pModel <- ans[[1]][[1]]
for (i in seq(along=(ans[[2]]))) ## ans[[2]] is a list of lists of
## pointers to effects. Each list corresponds to one
@@ -455,7 +433,13 @@
z$maxlikeTargets <- rowSums(ans)
z$maxlikeTargets2 <- ans
z$mult <- x$mult
- z$nrunMH <- z$mult * sum(z$maxlikeTargets[z$effects$basicRate])
+ z$nrunMH <-
+ z$mult * colSums(z$maxlikeTargets2[z$effects$basicRate, ,
+ drop=FALSE ])
+ z$nrunMH < pmax(z$nrunMH, 2)
+ ## make the number pretty
+ z$nrunMH <- ifelse (z$nrunMH > 100,
+ 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
@@ -549,7 +533,7 @@
}
else ## set up the initial chains in the sub processes
{
- ans <- .Call("mlInitializeSubProcesses",
+ ans <- .Call("mlInitializeSubProcesses",
PACKAGE=pkgname, pData, pModel,
simpleRates, z$probs, z$prmin, z$prmib,
x$minimumPermutationLength,
@@ -1851,4 +1835,39 @@
effects
}
-
+##@updateTheta siena07 Copy theta values from previous fit
+updateTheta <- function(effects, prevAns)
+{
+ if (!inherits(effects, "data.frame"))
+ {
+ stop("effects is not a data.frame")
+ }
+ if (!inherits(prevAns, "sienaFit"))
+ {
+ stop("prevAns is not an RSiena fit object")
+ }
+ prevEffects <- prevAns$requestedEffects
+ prevEffects$initialValue <- prevAns$theta
+ if (prevAns$cconditional)
+ {
+ condEffects <- attr(prevAns$f, "condEffects")
+ condEffects$initialValue <- prevAns$rate
+ prevEffects <- rbind(prevEffects, condEffects)
+ }
+ oldlist <- apply(prevEffects, 1, function(x)
+ paste(x[c("name", "shortName",
+ "type", "groupName",
+ "interaction1", "interaction2",
+ "period")],
+ collapse="|"))
+ efflist <- apply(effects, 1, function(x)
+ paste(x[c("name", "shortName",
+ "type", "groupName",
+ "interaction1", "interaction2",
+ "period")],
+ collapse="|"))
+ use <- efflist %in% oldlist
+ effects$initialValue[use] <-
+ prevEffects$initialValue[match(efflist, oldlist)][use]
+ effects
+}
Modified: pkg/RSiena/R/maxlikec.r
===================================================================
--- pkg/RSiena/R/maxlikec.r 2011-10-27 13:28:48 UTC (rev 177)
+++ pkg/RSiena/R/maxlikec.r 2011-10-27 17:30:41 UTC (rev 178)
@@ -12,15 +12,16 @@
##@maxlikec siena07 ML Simulation Module
maxlikec <- function(z, x, INIT=FALSE, TERM=FALSE, initC=FALSE, data=NULL,
effects=NULL, profileData=FALSE, prevAns=NULL,
- returnDeps=FALSE, returnChains=FALSE, byGroup=FALSE,
- returnDataFrame=FALSE)
+ returnChains=FALSE, byGroup=FALSE,
+ returnDataFrame=FALSE, byWave=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)
+ profileData=profileData, returnDeps=FALSE)
z$returnDataFrame <- returnDataFrame
z$returnChains <- returnChains
+ z$byWave <- byWave
if (initC)
{
return(NULL)
@@ -40,14 +41,14 @@
######################################################################
## retrieve stored information
f <- FRANstore()
- if (z$Phase == 2)
- {
- returnDeps <- FALSE
- }
- else
- {
- returnDeps <- z$returnDeps
- }
+ ##if (z$Phase == 2)
+ ##{
+ ## returnDeps <- FALSE
+ ##}
+ ##else
+ ##{
+ ## returnDeps <- z$returnDeps
+ ##}
callGrid <- z$callGrid
## z$int2 is the number of processors if iterating by period, so 1 means
## we are not. Can only parallelize by period at the moment.
@@ -61,9 +62,9 @@
{
theta <- z$theta
}
- ans <- .Call('mlPeriod', PACKAGE=pkgname, z$Deriv, f$pData,
+ ans <- .Call("mlPeriod", PACKAGE=pkgname, z$Deriv, f$pData,
f$pModel, f$myeffects, theta,
- returnDeps, 1, 1, z$nrunMH, z$addChainToStore,
+ 1, 1, z$nrunMH, z$addChainToStore,
z$needChangeContributions, z$returnDataFrame,
z$returnChains)
ans[[6]] <- list(ans[[6]])
@@ -74,24 +75,26 @@
ans[[9]] <- list(ans[[9]])
ans[[10]] <- list(ans[[10]])
}
- }
+ }
else
{
if (z$int2 == 1)
{
- anss <- apply(callGrid, 1, doMLModel, z$Deriv, z$thetaMat,
- returnDeps, z$nrunMH, z$addChainToStore,
+ 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)
}
else
{
use <- 1:(min(nrow(callGrid), z$int2))
- anss <- parRapply(z$cl[use], callGrid, doMLModel, z$Deriv,
- z$thetaMat,
- returnDeps, z$nrunMH, z$addChainToStore,
+ anss <- parRapply(z$cl[use], cbind(callGrid, 1:nrow(callGrid)),
+ doMLModel, z$Deriv, z$thetaMat,
+ z$nrunMH, z$addChainToStore,
z$needChangeContributions,
- z$returnDataFrame, z$returnChains, byGroup, z$theta)
+ z$returnDataFrame, z$returnChains, byGroup,
+ z$theta)
}
## reorganize the anss so it looks like the normal one
ans <- NULL
@@ -117,8 +120,7 @@
## {
## ans[[6]] <- NULL
## }
- ## browser()
- if (returnChains)
+ if (z$returnChains)
{
##TODO put names on these?
fff <- lapply(anss, function(x) x[[6]][[1]])
@@ -129,9 +131,9 @@
ans[[8]] <- lapply(anss, "[[", 8)
ans[[9]] <- lapply(anss, "[[", 9)
ans[[10]] <- lapply(anss, "[[", 10)
- if (!byGroup)
- {
- ans[[8]] <- Reduce("+", ans[[8]]) ## accepts
+ if (!byGroup)
+ {
+ ans[[8]] <- Reduce("+", ans[[8]]) ## accepts
ans[[9]] <- Reduce("+", ans[[9]]) ## rejects
ans[[10]] <- Reduce("+", ans[[10]]) ## aborts
}
@@ -140,12 +142,14 @@
fra <- -t(ans[[1]]) ##note sign change
FRANstore(f)
-
- ## if (returnDeps)
- ## sims <- ans[[6]]
- ## else
- sims <- NULL
-
+ #if (returnDeps || z$returnChains)
+ #{
+ # sims <- ans[[6]]
+ #}
+ #else
+ #{
+ sims <- NULL
+ #}
##print(length(sims[[1]]))
##if (returnDeps) ## this is not finished yet!
##{
@@ -184,7 +188,7 @@
}
##@doMLModel Maximum likelihood
-doMLModel <- function(x, Deriv, thetaMat, returnDeps, nrunMH, addChainToStore,
+doMLModel <- function(x, Deriv, thetaMat, nrunMH, addChainToStore,
needChangeContributions, returnDataFrame, returnChains,
byGroup, theta)
{
@@ -198,8 +202,8 @@
theta <- theta
}
.Call("mlPeriod", PACKAGE=pkgname, Deriv, f$pData,
- f$pModel, f$myeffects, theta, returnDeps,
- as.integer(x[1]), as.integer(x[2]), nrunMH, addChainToStore,
+ f$pModel, f$myeffects, theta,
+ as.integer(x[1]), as.integer(x[2]), nrunMH[x[3]], addChainToStore,
needChangeContributions, returnDataFrame, returnChains)
}
@@ -209,9 +213,14 @@
## of the matrix. Need to be put into a symmetric matrix.
## Tricky part is getting the rates in the right place!
## Note that we do not yet deal with rate effects other than the basic
- dff <- matrix(0, nrow=z$pp, ncol=z$pp)
+ dff <- as(Matrix(0, z$pp, z$pp), "dsyMatrix")
nPeriods <- length(derivList)
- dff2 <- array(0, dim=c(nPeriods, z$pp, z$pp))
+ dff2 <- vector("list", nPeriods)
+ if (z$byWave)
+ {
+ tmp <- as(Matrix(0, z$pp, z$pp), "dsyMatrix")
+ dff2[] <- tmp
+ }
for (period in 1:nPeriods)
{
dffraw <- derivList[[period]]
@@ -250,7 +259,10 @@
rawsub <- rawsub + nonRates * nonRates
dffPeriod <- dffPeriod + t(dffPeriod)
diag(dffPeriod) <- diag(dffPeriod) / 2
- dff2[period , , ] <- dff2[period, , ] - dffPeriod
+ if (z$byWave)
+ {
+ dff2[[period]] <- dff2[[period]] - dffPeriod
+ }
dff <- dff - dffPeriod
}
}
Modified: pkg/RSiena/R/phase1.r
===================================================================
--- pkg/RSiena/R/phase1.r 2011-10-27 13:28:48 UTC (rev 177)
+++ pkg/RSiena/R/phase1.r 2011-10-27 17:30:41 UTC (rev 178)
@@ -42,20 +42,7 @@
endNit <- z$n1 + int - (z$n1 - firstNit) %% int
}
z$n1 <- endNit
- z$sf <- matrix(0, nrow = z$n1, ncol = z$pp)
- z$sf2 <- array(0, dim=c(z$n1, f$observations - 1, z$pp))
- if (!x$maxlike & !z$FinDiff.method)
- {
- z$ssc <- array(0, dim=c(z$n1, f$observations - 1, z$pp))
- }
- else
- {
- z$sdf <- array(0, dim=c(z$n1, z$pp, z$pp))
- z$sdf2 <- array(0, dim=c(z$n1, f$observations -1, z$pp, z$pp))
- }
- z$accepts <- matrix(0, nrow=z$n1, ncol=7)
- z$rejects <- matrix(0, nrow=z$n1, ncol=7)
- z$npos <- rep(0, z$pp)
+ z <- createSiena07stores(z, z$n1, f)
z$writefreq <- 10
z$DerivativeProblem <- FALSE
z$Deriv <- !z$FinDiff.method ## can do both in phase 3 but not here!
@@ -132,7 +119,8 @@
use <- !z$fixed & npos < 5
z$epsilon[use] <- pmin(100.0 * z$scale[use], z$epsilon[use])
z$epsilon[use] <- pmax(0.1 * z$scale[use], z$epsilon[use])
- Report(c("New epsilon =", paste(" ", z$epsilon[use], collapse=""), ".\n"), sep="", cf, fill=80)
+ Report(c("New epsilon =", paste(" ", z$epsilon[use], collapse=""),
+ ".\n"), sep="", cf, fill=80)
if (z$repeatsForEpsilon <= 4)
{
Report("Change value of epsilon and restart Phase 1.\n", cf)
@@ -168,9 +156,8 @@
##@doPhase1it siena07 does 1 iteration in Phase 1
doPhase1it<- function(z, x, zsmall, xsmall, ...)
{
- int <- z$int
DisplayIteration(z)
- if (int == 1)
+ if (z$int == 1)
{
zz <- x$FRAN(zsmall, xsmall)
if (!zz$OK)
@@ -186,64 +173,10 @@
{
zz <- clusterCall(z$cl, usesim, zsmall, xsmall)
z$n <- z$n + z$int
- z$phase1Its <- z$phase1Its + int
+ z$phase1Its <- z$phase1Its + z$int
# browser()
}
- ## browser()
- if (int == 1)
- {
- fra <- colSums(zz$fra)
- fra <- fra - z$targets
- fra2 <- zz$fra
- z$sf[z$nit, ] <- fra
- z$sf2[z$nit, , ] <- zz$fra
- z$sims[[z$nit]] <- zz$sims
- z$chain[[z$nit]] <- zz$chain
- fra <- fra + z$targets
- }
- else
- {
- for (i in 1:int)
- {
- fra <- colSums(zz[[i]]$fra)
- fra <- fra - z$targets
- z$sf[z$nit + (i - 1), ] <- fra
- z$sf2[z$nit + (i - 1), , ] <- zz[[i]]$fra
- z$sims[[z$nit + (i - 1)]] <- zz[[i]]$sims
- }
- fra2 <- t(sapply(zz, function(x)x$fra))
- dim(fra2) <- c(int, nrow(zz[[1]]$fra), z$pp)
- fra <- t(sapply(zz, function(x) colSums(x$fra)))
- }
- if (x$maxlike)
- {
- z$sdf[z$nit, , ] <- zz$dff
- z$sdf2[z$nit, , , ] <- zz$dff2
- z$accepts[z$nit, ] <- zz$accepts
- z$rejects[z$nit, ] <- zz$rejects
- }
- else if (z$FinDiff.method)
- {
- z <- FiniteDifferences(z, x, fra , fra2, ...)
- z$sdf[z$nit:(z$nit + (z$int - 1)), , ] <- z$sdf0
- z$sdf2[z$nit:(z$nit + (z$int - 1)), , , ] <- z$sdf02
- }
- else
- {
- if (int==1)
- {
- if (!is.null(zz[['sc']]))
- z$ssc[z$nit , ,] <- zz$sc
- }
- else
- {
- for (i in 1:int)
- {
- if (!is.null(zz[[i]][['sc']]))
- z$ssc[z$nit + (i - 1), , ] <- zz[[i]]$sc
- }
- }
- }
+ z <- updateSiena07stores(z, zz, x)
CheckBreaks()
if (UserInterruptFlag() || UserRestartFlag())
{
@@ -263,7 +196,7 @@
progress <- val / z$pb$pbmax * 100
if (z$nit <= 5 || z$nit %% z$writefreq == 0 || z$nit %%5 == 0 ||
x$maxlike || z$FinDiff.method ||
- (int > 1 && z$nit %% z$writefreq < int))
+ (z$int > 1 && z$nit %% z$writefreq < z$int))
{
# Report(c('Phase', z$Phase, 'Iteration ', z$nit, '\n'))
if (is.batch())
@@ -274,7 +207,7 @@
else
{
DisplayTheta(z)
- DisplayDeviations(z, fra)
+ DisplayDeviations(z, z$sf[z$nit, ])
}
}
#browser()
@@ -409,6 +342,9 @@
z$phase1sdf <- z$sdf
z$phase1sdf2 <- z$sdf2
z$phase1scores <- z$ssc
+ z$phase1accepts <- z$accepts
+ z$phase1rejects <- z$rejects
+ z$phase1aborts <- z$aborts
##browser()
z
}
@@ -418,8 +354,8 @@
{
if (z$FinDiff.method || x$maxlike)
{
- dfra <- t(apply(z$sdf, c(2, 3), mean))
- }
+ dfra <- t(as.matrix(Reduce("+", z$sdf) / length(z$sdf)))
+ }
else
{
##note that warnings is never set as this piece of code is not executed
@@ -430,7 +366,7 @@
z$jacobianwarn1 <- rep(FALSE, z$pp)
if (any(diag(dfra) <= 0))
- {
+ {
for (i in 1 : z$pp)
{
if (dfra[i, i] < 0)
@@ -485,23 +421,26 @@
}
##@FiniteDifferences siena07 Does the extra iterations for finite differences
-FiniteDifferences <- function(z, x, fra, fra2, ...)
+FiniteDifferences <- function(z, x, fra, fra2)
{
int <- z$int
fras <- array(0, dim = c(int, z$pp, z$pp))
- fras2 <- array(0, dim = c(int, z$f$observations - 1, z$pp, z$pp))
+ if (z$byWave)
+ {
+ fras2 <- array(0, dim = c(int, z$f$observations - 1, z$pp, z$pp))
+ }
+ else
+ {
+ fras2 <- NULL
+ }
xsmall <- NULL
- ##browser()
for (i in 1 : z$pp)
{
- # zdummy <- z[c('theta', 'Deriv', 'cconditional', 'FinDiff.method',
- # 'int2', 'cl')]
zdummy <- makeZsmall(z)
if (z$Phase == 3 || !z$fixed[i])
{
zdummy$theta[i] <- z$theta[i] + z$epsilon[i]
}
- ## assign('.Random.seed',randomseed,pos=1)
if (int == 1)
{
zz <- x$FRAN(zdummy, xsmall, INIT=FALSE, fromFiniteDiff=TRUE)
@@ -515,23 +454,27 @@
{
zz <- clusterCall(z$cl, usesim, zdummy, xsmall,
INIT=FALSE, fromFiniteDiff=TRUE)
- #browser()
}
if (int == 1)
{
fras[1, i, ] <- colSums(zz$fra) - fra
- fras2[1, , i, ] <- zz$fra - fra2
+ if (z$byWave)
+ {
+ fras2[1, , i, ] <- zz$fra - fra2
+ }
}
else
{
for (j in 1 : int)
{
fras[j, i, ] <- colSums(zz[[j]]$fra) - fra[j, ]
- fras2[j, , i, ] <- zz[[j]]$fra - fra2[j, , ]
+ if (z$byWave)
+ {
+ fras2[j, , i, ] <- zz[[j]]$fra - fra2[j, , ]
+ }
}
}
}
- ##browser()
if (z$Phase == 1 && z$nit <= 10)
{
for (ii in 1: min(10 - z$nit + 1, int))
@@ -541,7 +484,10 @@
}
}
z$sdf0 <- fras / rep(rep(z$epsilon, each=int), z$pp)
- z$sdf02 <- fras2 / rep(rep(z$epsilon, each=int * dim(fras2)[2]), z$pp)
+ if (z$byWave)
+ {
+ z$sdf02 <- fras2 / rep(rep(z$epsilon, each=int * dim(fras2)[2]), z$pp)
+ }
z
}
##@derivativeFromScoresAndDeviations siena07 create dfra from scores and deviations
@@ -591,5 +537,118 @@
zsmall$needChangeContributions <- z$needChangeContributions
zsmall$callGrid <- z$callGrid
zsmall$thetaMat <- z$thetaMat
+ zsmall$byWave <- z$byWave
zsmall
}
+
+##@createSiena07Stores siena07 set up the storage areas used in phase 1 and 3
+createSiena07stores <- function(z, nIterations, f)
+{
+ z$sf <- matrix(0, nrow = nIterations , ncol = z$pp)
+ z$sf2 <- array(0, dim=c(nIterations, f$observations - 1, z$pp))
+ if (!z$maxlike && !z$FinDiff.method)
+ {
+ z$ssc <- array(0, dim=c(nIterations, f$observations - 1, z$pp))
+ }
+ else
+ {
+ z$sdf <- vector("list", nIterations)
+ z$sdf2 <- vector("list", nIterations)
+ }
+ ## misdat steps are separated out giving 9 types
+ z$accepts <- array(0, dim=c(nIterations, z$nDependentVariables, 9))
+ z$rejects <- array(0, dim=c(nIterations, z$nDependentVariables, 9))
+ z$aborts <- array(0, dim=c(nIterations, z$nDependentVariables, 9))
+ z$npos <- rep(0, z$pp)
+ if (!is.null(z$cconditional) && z$cconditional)
+ {
+ z$ntim <- matrix(NA, nrow=nIterations, ncol=f$observations - 1)
+ }
+ z$sims <- vector("list", nIterations)
+ z
+}
+##@updateSiena07Stores siena07 store data in phase 1 and 3
+updateSiena07stores <- function(z, zz, x)
+{
+ int <- z$int
+ if (int == 1)
+ {
+ fra <- colSums(zz$fra)
+ fra <- fra - z$targets
+ fra2 <- zz$fra
+ z$sf[z$nit, ] <- fra
+ z$sf2[z$nit, , ] <- zz$fra
+ z$sims[[z$nit]] <- zz$sims
+ z$chain[[z$nit]] <- zz$chain
+ fra <- fra + z$targets
+ }
+ else
+ {
+ for (i in 1:int)
+ {
+ fra <- colSums(zz[[i]]$fra)
+ fra <- fra - z$targets
+ z$sf[z$nit + (i - 1), ] <- fra
+ z$sf2[z$nit + (i - 1), , ] <- zz[[i]]$fra
+ z$sims[[z$nit + (i - 1)]] <- zz[[i]]$sims
+ }
+ fra2 <- t(sapply(zz, function(x)x$fra))
+ dim(fra2) <- c(int, nrow(zz[[1]]$fra), z$pp)
+ fra <- t(sapply(zz, function(x) colSums(x$fra)))
+ }
+ if (x$maxlike)
+ {
+ z$sdf[[z$nit]] <- zz$dff
+ z$sdf2[[z$nit]] <- zz$dff2
+ z$accepts[z$nit, , ] <- zz$accepts
+ z$rejects[z$nit, , ] <- zz$rejects
+ z$aborts[z$nit, , ] <- zz$aborts
+ }
+ else if (z$FinDiff.method)
+ {
+ z <- FiniteDifferences(z, x, fra, fra2)
+ for (i in 0:(z$int - 1))
+ {
+ z$sdf[[z$nit + i]] <- z$sdf0[i + 1, , ]
+ if (z$byWave)
+ {
+ z$sdf2[[z$nit + i]] <- z$sdf02[i + 1, , , ]
+ }
+ }
+ }
+ else
+ {
+ if (int==1)
+ {
+ if (!is.null(zz[['sc']]))
+ {
+ z$ssc[z$nit , ,] <- zz$sc
+ }
+ }
+ else
+ {
+ for (i in 1:int)
+ {
+ if (!is.null(zz[[i]][['sc']]))
+ {
+ z$ssc[z$nit + (i - 1), , ] <- zz[[i]]$sc
+ }
+ }
+ }
+ }
+ if ((!x$maxlike) && z$cconditional && z$Phase == 3)
+ {
+ if (int == 1)
+ {
+ z$ntim[z$nit, ] <- zz$ntim0
+ }
+ else
+ {
+ for (i in 1:int)
+ {
+ z$ntim[z$nit + (i-1), ] <- zz[[i]]$ntim0
+ }
+ }
+ }
+ z
+}
Modified: pkg/RSiena/R/phase3.r
===================================================================
--- pkg/RSiena/R/phase3.r 2011-10-27 13:28:48 UTC (rev 177)
+++ pkg/RSiena/R/phase3.r 2011-10-27 17:30:41 UTC (rev 178)
@@ -22,24 +22,8 @@
z$returnDeps <- z$returnDepsStored
if (x$checktime) z$ctime <- proc.time()[3]
+
## fix up iteration numbers if using multiple processors
- # if (10 %% int == 0)
- # {
- # firstNit <- 10
- # }
- # else
- # {
- # firstNit <- 10 + int - 10 %% int
- # }
- # if ((x$n3 - firstNit) %% int == 0)
- # {
- # endNit <- x$n3
- # }
- # else
- # {
- # endNit <- x$n3 + int - (x$n3 - firstNit) %% int
- # }
- # cat("endNit", endNit, "n3 ", x$n3, "\n")
if (x$n3 %% int > 0)
{
endNit <- x$n3 + int - x$n3 %%int
@@ -48,17 +32,8 @@
{
endNit <- x$n3
}
- z$n3 <- endNit
- z$sf <- matrix(0, nrow = z$n3, ncol = z$pp)
- z$sf2 <- array(0, dim = c(z$n3, f$observations - 1, z$pp))
- z$accepts3 <- matrix(0, nrow=z$n3, ncol=7)
- z$rejects3 <- matrix(0, nrow=z$n3, ncol=7)
- z$aborted3 <- matrix(0, nrow=z$n3, ncol=7)
- if (!is.null(z$cconditional) && z$cconditional)
- {
- z$ntim <- matrix(NA, nrow=z$n3, ncol=f$observations - 1)
- }
- z$sims <- vector("list", z$n3)
+ z$n3 <- endNit
+
## revert to original requested method for phase 3 unless symmetric
if (z$FinDiff.method && !x$FinDiff.method &&
(!is.null(z$FinDiffBecauseSymmetric)) && z$FinDiffBecauseSymmetric)
@@ -77,15 +52,7 @@
else
Report('Estimation of derivatives by the LR method (type 1).\n\n', outf)
- if (!x$maxlike & !z$FinDiff.method)
- {
- z$ssc <- array(0, dim = c(z$n3, f$observations - 1, z$pp))
- }
- else
- {
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/rsiena -r 178
More information about the Rsiena-commits
mailing list