[Rsiena-commits] r108 - in pkg: RSiena RSiena/R RSiena/data RSiena/inst/doc RSiena/man RSiena/src RSiena/src/data RSiena/src/model RSiena/src/model/effects RSiena/src/model/ml RSiena/src/model/variables RSiena/tests RSienaTest RSienaTest/R RSienaTest/doc RSienaTest/inst/doc RSienaTest/man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Jun 22 00:12:15 CEST 2010
Author: ripleyrm
Date: 2010-06-22 00:12:12 +0200 (Tue, 22 Jun 2010)
New Revision: 108
Added:
pkg/RSiena/src/model/effects/AltersCovariateAverageEffect.cpp
pkg/RSiena/src/model/effects/AltersCovariateAverageEffect.h
pkg/RSiena/src/model/effects/CovariateAndNetworkBehaviorEffect.cpp
pkg/RSiena/src/model/effects/CovariateAndNetworkBehaviorEffect.h
pkg/RSiena/src/model/effects/CovariateDistance2AlterEffect.cpp
pkg/RSiena/src/model/effects/CovariateDistance2AlterEffect.h
pkg/RSiena/src/model/effects/CovariateDistance2NetworkEffect.cpp
pkg/RSiena/src/model/effects/CovariateDistance2NetworkEffect.h
pkg/RSiena/src/model/effects/CovariateDistance2SimilarityEffect.cpp
pkg/RSiena/src/model/effects/CovariateDistance2SimilarityEffect.h
Modified:
pkg/RSiena/DESCRIPTION
pkg/RSiena/R/bayes.r
pkg/RSiena/R/effects.r
pkg/RSiena/R/effectsDocumentation.r
pkg/RSiena/R/effectsMethods.r
pkg/RSiena/R/globals.r
pkg/RSiena/R/maxlikec.r
pkg/RSiena/R/phase1.r
pkg/RSiena/R/phase2.r
pkg/RSiena/R/phase3.r
pkg/RSiena/R/print01Report.r
pkg/RSiena/R/printInitialDescription.r
pkg/RSiena/R/sienaDataCreate.r
pkg/RSiena/R/sienaTimeTest.r
pkg/RSiena/R/simstatsc.r
pkg/RSiena/changeLog
pkg/RSiena/data/allEffects.csv
pkg/RSiena/inst/doc/effects.pdf
pkg/RSiena/inst/doc/s_man400.pdf
pkg/RSiena/man/RSiena-package.Rd
pkg/RSiena/man/coDyadCovar.Rd
pkg/RSiena/man/print.sienaEffects.Rd
pkg/RSiena/man/siena01Gui.Rd
pkg/RSiena/man/siena07.Rd
pkg/RSiena/man/sienaDataCreate.Rd
pkg/RSiena/man/sienaDataCreateFromSession.Rd
pkg/RSiena/man/sienaModelCreate.Rd
pkg/RSiena/man/sienaNet.Rd
pkg/RSiena/man/sienaNodeSet.Rd
pkg/RSiena/man/sienaTimeTest.Rd
pkg/RSiena/man/simstats0c.Rd
pkg/RSiena/man/varDyadCovar.Rd
pkg/RSiena/src/data/BehaviorLongitudinalData.cpp
pkg/RSiena/src/data/BehaviorLongitudinalData.h
pkg/RSiena/src/data/Covariate.cpp
pkg/RSiena/src/data/Covariate.h
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/effects/AllEffects.h
pkg/RSiena/src/model/effects/BehaviorEffect.cpp
pkg/RSiena/src/model/effects/BehaviorEffect.h
pkg/RSiena/src/model/effects/CovariateDependentNetworkEffect.cpp
pkg/RSiena/src/model/effects/CovariateDependentNetworkEffect.h
pkg/RSiena/src/model/effects/EffectFactory.cpp
pkg/RSiena/src/model/ml/BehaviorChange.cpp
pkg/RSiena/src/model/ml/BehaviorChange.h
pkg/RSiena/src/model/ml/Chain.cpp
pkg/RSiena/src/model/ml/Chain.h
pkg/RSiena/src/model/ml/MLSimulation.cpp
pkg/RSiena/src/model/ml/MiniStep.cpp
pkg/RSiena/src/model/ml/MiniStep.h
pkg/RSiena/src/model/ml/NetworkChange.cpp
pkg/RSiena/src/model/ml/NetworkChange.h
pkg/RSiena/src/model/variables/BehaviorVariable.cpp
pkg/RSiena/src/model/variables/BehaviorVariable.h
pkg/RSiena/src/model/variables/DependentVariable.h
pkg/RSiena/src/model/variables/NetworkVariable.cpp
pkg/RSiena/src/model/variables/NetworkVariable.h
pkg/RSiena/src/siena07.cpp
pkg/RSiena/tests/parallel.R
pkg/RSiena/tests/parallel.Rout.save
pkg/RSienaTest/DESCRIPTION
pkg/RSienaTest/R/bayes.r
pkg/RSienaTest/R/effectsDocumentation.r
pkg/RSienaTest/R/effectsMethods.r
pkg/RSienaTest/R/print01Report.r
pkg/RSienaTest/R/printInitialDescription.r
pkg/RSienaTest/R/sienaDataCreate.r
pkg/RSienaTest/R/simstatsc.r
pkg/RSienaTest/changeLog
pkg/RSienaTest/doc/RSienaDeveloper.tex
pkg/RSienaTest/doc/s_man400.tex
pkg/RSienaTest/inst/doc/effects.pdf
pkg/RSienaTest/inst/doc/s_man400.pdf
pkg/RSienaTest/man/RSiena-package.Rd
pkg/RSienaTest/man/coDyadCovar.Rd
pkg/RSienaTest/man/print.sienaEffects.Rd
pkg/RSienaTest/man/siena01Gui.Rd
pkg/RSienaTest/man/siena07.Rd
pkg/RSienaTest/man/sienaDataCreate.Rd
pkg/RSienaTest/man/sienaDataCreateFromSession.Rd
pkg/RSienaTest/man/sienaModelCreate.Rd
pkg/RSienaTest/man/sienaNet.Rd
pkg/RSienaTest/man/sienaNodeSet.Rd
pkg/RSienaTest/man/simstats0c.Rd
pkg/RSienaTest/man/varDyadCovar.Rd
Log:
New covariate distance 2 effects, few bug fixes, amendments to help files.
Modified: pkg/RSiena/DESCRIPTION
===================================================================
--- pkg/RSiena/DESCRIPTION 2010-06-21 17:00:38 UTC (rev 107)
+++ pkg/RSiena/DESCRIPTION 2010-06-21 22:12:12 UTC (rev 108)
@@ -1,8 +1,8 @@
Package: RSiena
Type: Package
Title: Siena - Simulation Investigation for Empirical Network Analysis
-Version: 1.0.11.105
-Date: 2010-06-18
+Version: 1.0.11.108
+Date: 2010-06-21
Author: Various
Depends: R (>= 2.9.0), xtable
Imports: Matrix
Modified: pkg/RSiena/R/bayes.r
===================================================================
--- pkg/RSiena/R/bayes.r 2010-06-21 17:00:38 UTC (rev 107)
+++ pkg/RSiena/R/bayes.r 2010-06-21 22:12:12 UTC (rev 108)
@@ -10,7 +10,7 @@
## ****************************************************************************/
##@bayes algorithm fit a Bayesian model
bayes <- function(data, effects, model, nwarm=100, nmain=100, nrunMHBatches=20,
- nrunMH=100)
+ nrunMH=100, save=TRUE)
{
createStores <- function()
{
@@ -138,14 +138,14 @@
## initialise
Report(openfiles=TRUE, type="n") #initialise with no file
- require(lattice)
- ## require(MASS)
- ## dev.new()
- ## histPlot = dev.cur()
- dev.new()
- thetaplot = dev.cur()
- dev.new()
- acceptsplot = dev.cur()
+ if (!save)
+ {
+ require(lattice)
+ dev.new()
+ thetaplot = dev.cur()
+ dev.new()
+ acceptsplot = dev.cur()
+ }
z <- NULL
z$FinDiff.method <- FALSE
z$maxlike <- TRUE
@@ -190,16 +190,9 @@
z <- storeData()
numm <- z$numm
- if (ii %% 10 == 0) ## do some plots
+ if (ii %% 10 == 0 && !save) ## do some plots
{
cat('main after ii',ii,numm, '\n')
- ## dev.set(histPlot)
- ## par(mfrow=c(2,3))
- ## truehist(z$lambdas)
- ## truehist(z$betas[, 1])
- ## truehist(z$candidates[, 1])
- ## truehist(z$betas[, 2])
- ## truehist(z$candidates[, 2])
dev.set(thetaplot)
thetadf <- data.frame(z$lambdas, z$betas)
acceptsdf <- data.frame(z$MHproportions, z$acceptances)
Modified: pkg/RSiena/R/effects.r
===================================================================
--- pkg/RSiena/R/effects.r 2010-06-21 17:00:38 UTC (rev 107)
+++ pkg/RSiena/R/effects.r 2010-06-21 22:12:12 UTC (rev 108)
@@ -684,6 +684,16 @@
substituteNames(covarBehObjEffects[2, ],
zName=names(xx$depvars)[j]))
}
+ if ((types[j] =="oneMode" &&
+ attr(xx$depvars[[j]], 'nodeSet') == nodeSet)
+ || (types[j] == "bipartite" &&
+ attr(xx$depvars[[j]], 'nodeSet')[2] == nodeSet))
+ {
+ covObjEffects <-
+ rbind(covObjEffects,
+ substituteNames(covarBehObjEffects[3, ],
+ zName=names(xx$depvars)[j]))
+ }
}
}
# if (!is.null(covObjEffects))
Modified: pkg/RSiena/R/effectsDocumentation.r
===================================================================
--- pkg/RSiena/R/effectsDocumentation.r 2010-06-21 17:00:38 UTC (rev 107)
+++ pkg/RSiena/R/effectsDocumentation.r 2010-06-21 22:12:12 UTC (rev 108)
@@ -13,7 +13,9 @@
effectsDocumentation <- function(type="html", display=type=="html",
filename="effects")
{
- x <- allEffects[, c(1:2,4:7, 15, 23)]
+ x <- allEffects[, c("effectGroup", "effectName", "shortName",
+ "endowment", "interaction1", "interaction2",
+ "parm", "interactionType")]
storage.mode(x$parm) <- "integer"
names(x)[4] <- "endow?"
names(x)[5] <- "inter1"
Modified: pkg/RSiena/R/effectsMethods.r
===================================================================
--- pkg/RSiena/R/effectsMethods.r 2010-06-21 17:00:38 UTC (rev 107)
+++ pkg/RSiena/R/effectsMethods.r 2010-06-21 22:12:12 UTC (rev 108)
@@ -10,7 +10,7 @@
## *
## ****************************************************************************/
##@print.sienaEffects Methods
-print.sienaEffects <- function(x, fileName=NULL, ...)
+print.sienaEffects <- function(x, fileName=NULL, includeOnly=TRUE,...)
{
if (!inherits(x, "sienaEffects"))
stop("not a legitimate Siena effects object")
@@ -25,9 +25,12 @@
nDependents <- length(unique(x$name))
userSpecifieds <- x$shortName[x$include] %in% c("unspInt", "behUnspInt")
endowments <- !x$type[x$include] %in% c("rate", "eval")
-
- specs <- x[x$include, c("name", "effectName", "include", "fix", "test",
+ specs <- x[, c("name", "effectName", "include", "fix", "test",
"initialValue", "parm")]
+ if (includeOnly)
+ {
+ specs <- specs[x$include, ]
+ }
if (nDependents == 1)
{
specs <- specs[, -1]
@@ -46,10 +49,15 @@
}
specs[, "initialValue"] <- format(round(specs$initialValue,digits=5),
width=10)
- row.names(specs) <- 1:nrow(specs)
-
- print(as.matrix(specs), quote=FALSE)
-
+ if (nrow(specs) > 0)
+ {
+ row.names(specs) <- 1:nrow(specs)
+ print(as.matrix(specs), quote=FALSE)
+ }
+ else
+ {
+ print.data.frame(specs)
+ }
}
else
{
Modified: pkg/RSiena/R/globals.r
===================================================================
--- pkg/RSiena/R/globals.r 2010-06-21 17:00:38 UTC (rev 107)
+++ pkg/RSiena/R/globals.r 2010-06-21 22:12:12 UTC (rev 108)
@@ -23,6 +23,7 @@
x <- x
beverbose <- verbose
besilent <- silent
+ noReportFile <- FALSE
function(txt, dest, fill=FALSE, sep=" ", hdest, openfiles=FALSE,
closefiles=FALSE, type=c("a", "w", "n"), projname="Siena" ,
verbose=FALSE, silent=FALSE)
@@ -40,6 +41,10 @@
{
x$outf <<- file(paste(projname, ".out", sep=""), open="a")
}
+ else if (type == "n")
+ {
+ noReportFile <<- TRUE
+ }
}
else if (closefiles)
@@ -69,7 +74,10 @@
}
else
{
- cat(txt, file = x[[hdest]], fill = fill, sep = sep)
+ if (!noReportFile)
+ {
+ cat(txt, file = x[[hdest]], fill = fill, sep = sep)
+ }
}
}
else
@@ -85,12 +93,18 @@
{
if (is.null(x[[deparse(substitute(dest))]]))
{
- cat(txt, fill=fill, sep=sep)
+ if (!besilent)
+ {
+ cat(txt, fill=fill, sep=sep)
+ }
}
else
{
- cat(txt, file=x[[deparse(substitute(dest))]],
- fill=fill, sep=sep)
+ if (!noReportFile)
+ {
+ cat(txt, file=x[[deparse(substitute(dest))]],
+ fill=fill, sep=sep)
+ }
}
}
}
Modified: pkg/RSiena/R/maxlikec.r
===================================================================
--- pkg/RSiena/R/maxlikec.r 2010-06-21 17:00:38 UTC (rev 107)
+++ pkg/RSiena/R/maxlikec.r 2010-06-21 22:12:12 UTC (rev 108)
@@ -33,7 +33,10 @@
f$pModel <- NULL
f$pData <- NULL
FRANstore(NULL) ## clear the stored object
- PrintReport(z, x)
+ if (is.null(z$print) || z$print)
+ {
+ PrintReport(z, x)
+ }
if (sum(z$test))
{
z$fra <- colMeans(z$sf, na.rm=TRUE)
@@ -41,7 +44,10 @@
z <- c(z, ans)
TestOutput(z, x)
}
- dimnames(z$dfra)[[1]] <- as.list(z$requestedEffects$shortName)
+ if (!is.null(z$dfra))
+ {
+ dimnames(z$dfra)[[1]] <- as.list(z$requestedEffects$shortName)
+ }
return(z)
}
######################################################################
@@ -50,14 +56,14 @@
## retrieve stored information
f <- FRANstore()
## browser()
- if (z$Phase == 2)
- {
- returnDeps <- FALSE
- }
- else
- {
- returnDeps <- f$returnDeps
- }
+ ##if (z$Phase == 2)
+ ##{
+ ## returnDeps <- FALSE
+ ##}
+ ##else
+ ##{
+ returnDeps <- z$returnDeps
+ ##}
## create a grid of periods with group names in case want to parallelize
## using this
groupPeriods <- attr(f, "groupPeriods")
@@ -68,10 +74,10 @@
## we are not
if (z$int2==1 || nrow(callGrid) == 1)
{
- ## for now!
- ans <- .Call('mlPeriod', PACKAGE=pkgname, z$Deriv, f$pData,
- f$pModel, f$myeffects, f$pMLSimulation, z$theta,
- returnDeps, 1, 1)
+ ## for now!
+ ans <- .Call('mlPeriod', PACKAGE=pkgname, z$Deriv, f$pData,
+ f$pModel, f$myeffects, f$pMLSimulation, z$theta,
+ returnDeps, 1, 1, z$nrunMH)
}
else
{
@@ -131,9 +137,9 @@
dffraw <- ans[[7]]
dff <- matrix(0, nrow=z$pp, ncol=z$pp)
start <- 1
+ rawsub <- 1
for (i in 1:length(f$myeffects))
{
- rawsub <- 1
rows <- nrow(f$myeffects[[i]])
nRates <- sum(f$myeffects[[i]]$type == 'rate')
nonRates <- rows - nRates
@@ -145,8 +151,10 @@
start : (start + nonRates - 1)] <-
dffraw[rawsub:(rawsub + nonRates * nonRates - 1)]
start <- start + nonRates
+ rawsub <- rawsub + nonRates * nonRates
}
- dff[upper.tri(dff)] <- dff[lower.tri(dff)]
+ dff <- dff + t(dff)
+ diag(dff) <- diag(dff) / 2
dff <- -dff
}
else
@@ -156,7 +164,18 @@
# browser()
list(fra = fra, ntim0 = NULL, feasible = TRUE, OK = TRUE,
- sims=sims, dff = dff, chain = ans[[6]], accepts=ans[[8]],
+ sims=sims, dff = dff, chain = list(ans[[6]]), accepts=ans[[8]],
rejects= ans[[9]])
}
+dist2full<-function(dis) {
+
+ n<-attr(dis,"Size")
+
+ full<-matrix(0,n,n)
+
+ full[lower.tri(full)]<-dis
+
+ full+t(full)
+
+}
Modified: pkg/RSiena/R/phase1.r
===================================================================
--- pkg/RSiena/R/phase1.r 2010-06-21 17:00:38 UTC (rev 107)
+++ pkg/RSiena/R/phase1.r 2010-06-21 22:12:12 UTC (rev 108)
@@ -568,5 +568,7 @@
zsmall$cconditional <- z$cconditional
zsmall$condvar <- z$condvar
zsmall$pp <- z$pp
+ zsmall$nrunMH <- z$nrunMH
+ zsmall$returnDeps <- z$returnDeps
zsmall
}
Modified: pkg/RSiena/R/phase2.r
===================================================================
--- pkg/RSiena/R/phase2.r 2010-06-21 17:00:38 UTC (rev 107)
+++ pkg/RSiena/R/phase2.r 2010-06-21 22:12:12 UTC (rev 108)
@@ -63,7 +63,7 @@
z
}
##@proc2subphase siena07 Do one subphase of phase 2
-proc2subphase<- function(z, x, subphase, ...)
+proc2subphase <- function(z, x, subphase, useAverage=TRUE, ...)
{
## init subphase of phase 2
z <- AnnouncePhase(z, x, subphase)
@@ -87,14 +87,15 @@
z$time1 <- proc.time()[3]
z$thav <- z$theta
z$thavn <- 1
- ## cat(z$thav, z$theta, '\n')
+ ## cat(z$thav, z$theta, '\n')
z$prod0 <- rep(0, z$pp)
z$prod1 <- rep(0, z$pp)
## ###############################################
## do the iterations for this repeat of this subphase
## ##############################################
+ ##z <- doIterationsCopy(z, x, subphase, ...)
z <- doIterations(z, x, subphase, ...)
- ## if (z$nit == 50) browser()
+ ## if (z$nit == 50) browser()
if (!z$OK || UserInterruptFlag() || UserRestartFlag() ||
EarlyEndPhase2Flag())
{
@@ -155,7 +156,20 @@
' = ', format(subphaseTime, nsmall=4, digits=4),
'\n', sep=''), lf)
}
- z$theta <- z$thav / z$thavn #(z$nit + 1)
+ if (useAverage)
+ {
+ z$theta <- z$thav / z$thavn #(z$nit + 1)
+ }
+ else
+ {
+ ## use regession
+ cat(z$thav / z$thavn, '\n') #(z$nit + 1)
+ mylm <- lm(z$sf[1:z$nit, ] ~ z$thetaStore[1:z$nit, ])
+ coefs <- coef(mylm)[-1, ]
+ newvals <- solve(t(coefs), - mylm$coef[1, ])
+ z$theta <- newvals
+ cat(z$theta, '\n')
+ }
DisplayThetaAutocor(z)
## cat('it',z$nit,'\n')
##recalculate autocor using -1 instead of -2 as error
@@ -192,6 +206,7 @@
ac <- 0
xsmall <- NULL
zsmall <- makeZsmall(z)
+ z$returnDeps <- FALSE
repeat
{
z$n <- z$n+1
@@ -340,3 +355,446 @@
}
z
}
+##@doIterationsCopy siena07 Do all iterations for 1 repeat of 1 subphase of phase 2
+doIterationsCopy <- function(z, x, subphase, numberIterations=10,
+ ...)
+{
+ ## designed to try things out!
+ if (z$cconditional)
+ {
+ stop("cannot use this routine with conditional estimation")
+ }
+ z$nit <- 0
+ ac <- 0
+ xsmall <- NULL
+ zsmall <- makeZsmall(z)
+ z <- setUpPhase2Storage(z, numberIterations)
+ repeat
+ {
+ z$n <- z$n+1
+ z$nit <- z$nit + 1
+ if (subphase == 1 && z$nit > 1)
+ z$time1 <- proc.time()[[3]]
+ if (subphase == 1 && z$nit > 10)
+ {
+ time1 <- proc.time()[[3]] - z$time1
+ if (time1 > 1e-5)
+ {
+ z$writefreq <- max(1, round(20.0 / time1))
+ }
+ else
+ {
+ z$writefreq <- 20
+ }
+ ## if (is.batch())
+ ## {
+ ## z$writefreq <- z$writefreq * 2 ##compensation for it
+ ## ## running faster with no tcl/tk
+ ## }
+ z$writefreq <- roundfreq(z$writefreq)
+ z$writefreq <- 1
+ }
+ if ((z$nit <= 10) || (z$nit %% z$writefreq ==0))
+ {
+ DisplayIteration(z)
+ if (is.batch())
+ {
+ val <- getProgressBar(z$pb)
+ increment <- ifelse(z$nit <= 10, 1, z$writefreq)
+ Report(paste('Phase ', z$Phase, ' Subphase ', subphase,
+ ' Iteration ', z$nit,' Progress: ',
+ round((increment + val) /
+ z$pb$pbmax * 100),
+ '%\n', sep = ''))
+ z$pb <- setProgressBar(z$pb, val + increment)
+ }
+ else
+ {
+ if (z$nit > 1)
+ {
+ DisplayDeviations(z, fra)
+ }
+ if (z$nit %% z$writefreq == 0)
+ {
+ val <- getProgressBar(z$pb)
+ z$pb <-setProgressBar(z$pb, val + z$writefreq)
+ }
+ }
+ }
+ zsmall$nit <- z$nit
+
+ zsmall$addChainToStore <- FALSE
+ zsmall$needChangeContributions <- FALSE
+
+ if (z$int == 1) ## not using a cluster
+ {
+ fra <- z$targets - z$targets
+ for (i in 1:numberIterations)
+ {
+ zz <- x$FRAN(zsmall, xsmall, returnDeps=TRUE)
+ ## browser()
+ fra0 <- colSums(zz$fra) - z$targets
+ fra <- fra + fra0
+ if (!zz$OK)
+ {
+ z$OK <- zz$OK
+ break
+ }
+ z <- storeChainsAndFra(z, zz, i, fra0)
+ }
+ if (!z$OK)
+ {
+ break
+ }
+ fra <- fra / numberIterations
+ z <- calculateLikelihoods(z)
+ }
+ else
+ {
+ zz <- clusterCall(z$cl, usesim, zsmall, xsmall)
+ fra <- sapply(zz, function(x) colSums(x$fra) - z$targets)
+ dim(fra) <- c(z$pp, z$int)
+ fra <- rowMeans(fra)
+ zz$OK <- sapply(zz, function(x) x$OK)
+ if (!all(zz$OK))
+ {
+ z$OK <- FALSE
+ break
+ }
+ }
+ if (x$maxlike)
+ {
+ z$phase2fras[subphase, ,z$nit] <- fra
+ ## z$rejectprops[subphase, , z$nit] <- zz$rejectprop
+ }
+ if (z$nit %% 2 == 1)
+ {
+ prev.fra <- fra
+ }
+ else
+ {
+ z$prod0 <- z$prod0 + fra * fra
+ z$prod1 <- z$prod1 + fra * prev.fra
+ ac <- ifelse (z$prod0 > 1e-12, z$prod1 / z$prod0, -2)
+ z$maxacor <- max(-99, ac[!z$fixed]) ##note -2 > -99
+ z$minacor <- min(1, ac[(!z$fixed) & ac > -1.0])
+ z$ac <- ac
+ if (z$nit %% z$writefreq == 0)
+ {
+ DisplayThetaAutocor(z)
+ }
+ }
+ z <- doChangeStep(z, x, fra)
+ zsmall$theta <- z$theta
+
+ if (x$maxlike && !is.null(x$moreUpdates) && x$moreUpdates > 0)
+ {
+ z <- doMoreUpdates(z, x, x$moreUpdates * subphase)
+ zsmall$theta <- z$theta
+ }
+ ## importance sampling steps
+ importSub <- 1
+ repeat
+ {
+ ## browser()
+ z <- predictOutcomes(z)
+ #cat(z$varWeights,'\n')
+ if (is.na(z$varWeights) ||
+ z$varWeights > var(c(rep(0, numberIterations-1), 1))/10)
+ {
+ break
+ }
+ z <- doChangeStep(z, x, z$predictedStats)
+ zsmall$theta <- z$theta
+ importSub <- importSub + 1
+ if (importSub > 25)
+ {
+ break
+ }
+ }
+ if (zsmall$addChainToStore)
+ {
+ clearStoredChains()
+ }
+ ##check for user interrupt
+ ## browser()
+ CheckBreaks()
+ if (UserInterruptFlag() || UserRestartFlag() || EarlyEndPhase2Flag())
+ {
+ break
+ }
+ ## do we stop?
+ if ( (z$nit >= z$n2min && z$maxacor < 1e-10)
+ || (z$nit >= z$n2max) || (z$nit >= 50 && z$minacor < -0.8 &&
+ z$repeatsubphase < z$maxrepeatsubphase))
+ {
+ break
+ }
+ }
+ z
+}
+
+##@setUpPhase2Storage siena07 create stores and values in nonstandard Phase2
+setUpPhase2Storage <- function(z, niter)
+{
+ nGroup <- z$f$nGroup
+ z$chains <- lapply(1:nGroup, function(x)
+ lapply(1:z$groupPeriods[x], function(y)
+ vector("list", niter)))
+ z$lik0 <- rep(0, niter * sum(z$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)
+ z
+}
+
+##@storeChainaAndFra siena07 update step for storage in non-standard phase 2
+storeChainsAndFra <- function(z, zz, i, fra)
+{
+ for (ii in 1:z$nGroup)
+ {
+ for (jj in 1:z$groupPeriods[ii])
+ {
+ z$chains[[ii]][[jj]][[i]] <- zz$chain[[ii]][[jj]]
+ }
+ }
+ z$iterFra[i, ] <- fra
+ z$thetaStore[z$nit, ] <- z$theta
+ z$sf[z$nit, ] <- fra
+ z
+}
+
+##@calculateLikelihoods siena07 for use in non-standard phase 2
+calculateLikelihoods <- function(z)
+{
+ storeSub <- 1
+ for (ii in 1:z$nGroup)
+ {
+ for (jj in 1:z$groupPeriods[ii])
+ {
+ chains <- z$chains[[ii]][[jj]]
+ for (chain in chains)
+ {
+ if (length(chain) == 0)
+ {
+ stop("no events with this theta")
+ }
+ z$lik0[storeSub] <-
+ getLikelihoodPhase2(chain, z$nactors[[ii]],
+ z$theta[z$rateParameterPosition[[ii]][[jj]]])
+ storeSub <- storeSub + 1
+ }
+ }
+ }
+ z
+}
+
+##@getLikelihoodPhase2 siena07 for use in non-standard phase 2
+getLikelihoodPhase2 <- function(chain, nactors, lambda)
+{
+ 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))
+ loglik <- loglik - sum(nactors * lambda) + sum(nc * log(lambda))
+ loglik
+}
+
+##@predictOutcomes siena07 for use in non-standard phase 2
+predictOutcomes <- function(z)
+{
+ ## now the likelihood for the new theta
+ ps <- getProbabilitiesPhase2(z$chain, z$theta, z$nactors,
+ z$rateParameterPosition, z$cl)$lik
+ ps <- exp(unlist(ps) - z$lik0)
+ ps <- ps / sum(ps)
+ z$predictedStats <- apply(z$iterFra, 2, function(x) sum(x * ps))
+ z$varWeights <- var(ps)
+ z
+}
+##@doChangeStep siena07 for use in non-standard phase 2, or outside siena07
+doChangeStep <- function(z, x, fra)
+{
+ ## limit change. Reporting is delayed to
+ ## end of phase.
+ ## browser()
+ if (x$diag)## !maxlike at present
+ {
+ maxrat<- max(ifelse(z$sd, abs(fra)/ z$sd, 1.0))#### check this
+ if (maxrat > x$maxmaxrat)
+ {
+ maxrat <- x$maxmaxrat / maxrat
+ z$truncated[z$nit] <- TRUE
+ }
+ else
+ maxrat <- 1.0
+ fchange<- z$gain * fra * maxrat / diag(z$dfra)
+ }
+ else
+ {
+ fchange <- as.vector(z$gain * fra %*% z$dinv)
+ }
+ ## browser()
+ fchange[z$fixed] <- 0.0
+ ## check positivity restriction
+ if (!is.null(z$nit))
+ {
+ z$positivized[z$nit, ] <- z$posj & (fchange > z$theta)
+ }
+ fchange <- ifelse(z$posj & (fchange > z$theta), z$theta * 0.5, fchange)
+
+ z$theta <- z$theta - fchange
+ if (!is.null(z$thav))
+ {
+ z$thav <- z$thav + z$theta
+ z$thavn <- z$thavn + 1
+ }
+ z
+}
+
+##@clearStoredChains algorithms Clears storage used for chains in C.
+clearStoredChains <- function()
+{
+ f <- FRANstore()
+ .Call("clearStoredChains", PACKAGE=pkgname, f$pModel)
+}
+
+##@getProbabilitiesPhase2 algorithms Recalculates change contributions
+getProbabilitiesPhase2 <- function(chain, theta, nactors, rateParameterPosition,
+ cl=NULL)
+{
+ f <- FRANstore()#
+ getScores <- FALSE
+ ##print(theta)
+ for (i in 1:length(chain)) # group
+ {
+ for (j in 1:length(chain[[i]])) # period
+ {
+ if (!is.null(cl) )
+ {
+ tmp <- parLapply(cl, chain[[i]][[j]],
+ function(x, i, j, theta, getScores, k, n)
+ {
+ f <- FRANstore()
+ resp <- .Call("getChainProbabilitiesList",
+ PACKAGE = pkgname, x,
+ f$pData, f$pModel, as.integer(i),
+ as.integer(j), f$myeffects,
+ theta, getScores)
+ lik <-
+ getLikelihoodPhase2(resp[[1]], nactors[[i]],
+ theta[k])
+ if (getScores)
+ {
+ sc <- resp[[2]]
+ }
+ else
+ {
+ sc <- NULL
+ }
+ list(lik=lik, sc=sc)
+ }, i=i, j=j, theta=theta, getScores=getScores,
+ k=rateParameterPosition[[i]][[j]],
+ n=nactors[[i]]
+ )
+ lik <- sapply(tmp, function(x)x[[1]])
+ if (getScores)
+ {
+ sc <- t(sapply(tmp, function(x)x[[2]]))
+ }
+ else
+ {
+ sc <- NULL
+ }
+ }
+ else
+ {
+ lik <- rep(0, length(chain[[i]][[j]]))
+ sc <- matrix(0, nrow=length(chain[[i]][[j]]),
+ ncol=length(theta))
+ for (k in 1:length(chain[[i]][[j]])) # one chain
+ {
+ resp <- .Call("getChainProbabilitiesList",
+ PACKAGE = pkgname,
+ chain[[i]][[j]][[k]],
+ f$pData, f$pModel, as.integer(i),
+ as.integer(j), f$myeffects, theta,
+ getScores)
+
+ lik[k] <-
+ getLikelihoodPhase2(resp[[1]], nactors[[i]],
+ theta[rateParameterPosition[[i]][[j]]])
+ if (getScores)
+ {
+ sc[k,] <- resp[[2]]
+ }
+ }
+ }
+ }
+ }
+ list(lik=lik, sc=sc)
+}
+##@initForAlgorithms siena07 stores values for use in nonstandard Phase2
+initForAlgorithms <- function(z)
+{
+ if (z$cconditional)
+ {
+ return(z)
+ }
+ nGroup <- z$f$nGroup
+ z$nDependentVariables <- length(z$f$depNames)
+ atts <- attributes(z$f)
+ z$groupPeriods <- atts$groupPeriods - 1
+ netnames <- names(z$f[[1]]$depvars)
+ z$nactors <- lapply(1:nGroup, function(i, periods, data)
+ {
+ tmp <- sapply(data[[i]]$depvars, function(x)
+ dim(x)[1])
+ tmp <- tmp[match(netnames, names(data[[i]]$depvars))]
+ tmp
+ }, periods=z$groupPeriods, data=z$f[1:nGroup]
+ )
+ z$rateParameterPosition <-
+ lapply(1:nGroup, function(i, periods, data)
+ {
+ lapply(1:periods[i], function(j)
+ {
+ rateEffects <-
+ z$effects[z$effects$basicRate &
+ z$effects$period == j &
+ z$effects$group == i,]
+ rateEffects <-
+ rateEffects[match(netnames,
+ rateEffects$name), ]
+ tmp <- as.numeric(row.names(rateEffects))
+ names(tmp) <- netnames
+ tmp
+ }
+ )
+ }, periods=z$groupPeriods, data=z$f[1:nGroup]
+ )
+ z$evalParameterPosition <-
+ lapply(netnames, function(x)
+ {
+ as.numeric(row.names(z$effects[z$effects$name == x &
+ z$effect$type =="eval", ]))
+ }
+ )
+ names(z$evalParameterPosition) <- netnames
+ z$endowParameterPosition <-
+ lapply(netnames, function(x)
+ {
+ as.numeric(row.names(z$effects[z$effects$name == x &
+ z$effect$type =="endow", ]))
+ }
+ )
+ z$basicRate <- z$effects$basicRate
+ z$nGroup <- nGroup
+ z
+}
Modified: pkg/RSiena/R/phase3.r
===================================================================
--- pkg/RSiena/R/phase3.r 2010-06-21 17:00:38 UTC (rev 107)
+++ pkg/RSiena/R/phase3.r 2010-06-21 22:12:12 UTC (rev 108)
@@ -19,6 +19,7 @@
DisplayTheta(z)
z$Phase <- 3
int <- z$int
+ z$returnDeps <- z$returnDepsStored
if (x$checktime) z$ctime <- proc.time()[3]
## fix up iteration numbers if using multiple processors
Modified: pkg/RSiena/R/print01Report.r
===================================================================
--- pkg/RSiena/R/print01Report.r 2010-06-21 17:00:38 UTC (rev 107)
+++ pkg/RSiena/R/print01Report.r 2010-06-21 22:12:12 UTC (rev 108)
@@ -33,7 +33,8 @@
}
}
##@reportDataObject internal print01Report
- reportDataObject <- function(x, periodFromStart=0, multi=FALSE)
+ reportDataObject <- function(x, periodFromStart=0, multi=FALSE,
+ session=session)
{
##@reportStart internal print01Report
reportStart <- function()
@@ -357,7 +358,6 @@
{
filename <-
session$Filename[session$Name == netname]
-
Report(c(" was read from file ", filename, ".\n"),
sep="", outf)
}
@@ -823,14 +823,16 @@
paste("Subproject ", i, ": <", names(data)[i], ">",
sep="", collapse="")
)
- reportDataObject(data[[i]], periodFromStart, multi=TRUE)
+ thisSession <- session[session$Group == names(data)[i],]
+ reportDataObject(data[[i]], periodFromStart, multi=TRUE,
+ session=thisSession)
periodFromStart <- periodFromStart + data[[i]]$observations
}
}
else
{
Heading(1, outf, "Data input.")
- reportDataObject(data[[1]], 0, multi=FALSE)
+ reportDataObject(data[[1]], 0, multi=FALSE, session=session)
}
atts <- attributes(data)
nets <- atts$types != "behavior"
@@ -1092,6 +1094,15 @@
length(atts$vCovars) > 0)
{
netnames <- atts$netnames
+ vCovarSim2 <-
+ lapply(data, function(x)
+ lapply(x$vCovars, function(y) attr(y, "simMeans")))
+ behSim2 <-
+ lapply(data, function(x)
+ lapply(x$depvars, function(y) attr(y, "simMeans")))
+ cCovarSim2 <-
+ lapply(data, function(x)
+ lapply(x$cCovars, function(y) attr(y, "simMeans")))
if (nData > 1)
{
vCovarSim <-
@@ -1176,9 +1187,104 @@
}
}
}
+ ## report on alter similarity means
+ Report(c("\nFor the distance-2 similarity variable calculated from",
+ "each actor \ncovariate, the mean is also subtracted.\n",
+ "These means are:\n"), outf)
+ if (nData == 1) ## ie we may have constant covariates
+ {
+ sims <- cCovarSim2[[1]]
+ for (i in seq(along=sims))
+ {
+ if (atts$cCovarPoszvar[i])
+ {
+ Report(c("Distance-2 Similarity ",
+ format(atts$cCovars[i], width=12),
+ ':', format(names(sims[[i]]), width=12,
+ justify="right"), '\n'), sep="", outf)
+ Report(c(rep(" ", 35),
+ format(round(sims[[i]], 4), width=12, nsmall=4),
+ '\n'), sep="", outf)
+ }
+ }
+ }
+ for (i in seq(along=atts$netnames))
+ {
+ if (atts$types[i] == "behavior" && atts$bPoszvar[i])
+ {
+ if (nData > 1)
+ {
+ thisSim <- lapply(behSim2, function(x)x[[netnames[i]]])
+ Report(c("Distance-2 Similarity ",
+ format(atts$netnames[i], width=12), ":\n"),
+ sep="", outf)
+ Report(c(rep(" ", 24),
+ format(names(thisSim[[1]]),
+ width=12, justify="right"), "\n"), sep="",
+ outf)
+ mystr <- format(paste(" Subproject ", 1:nData, " <",
+ atts$names, "> ", sep=""))
+ for (j in seq(along=thisSim))
+ {
+ Report(c(mystr[j], format(round(thisSim[[j]], 4),
+ nsmall=4, width=12), "\n"),
+ sep="", outf)
+ }
+ Report("\n", outf)
+ }
+ else
+ {
+ sims <- behSim2[[1]]
+ Report(c("Distance-2 Similarity ", format(atts$netnames[i],
+ width=12),
+ ':', format(names(sims[[i]]), width=12), '\n'),
+ sep="", outf)
+ Report(c(rep(" ", 35), format(round(sims[[i]], 4), width=12,
+ nsmall=4), '\n'), sep="", outf)
+ }
+ }
+ }
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/rsiena -r 108
More information about the Rsiena-commits
mailing list