[Rsiena-commits] r189 - in pkg: RSiena RSiena/R RSiena/data RSiena/inst/doc RSiena/man RSiena/po RSiena/src RSiena/src/model RSiena/src/model/ml RSienaTest RSienaTest/R RSienaTest/data RSienaTest/doc RSienaTest/inst/doc RSienaTest/inst/examples RSienaTest/man RSienaTest/po RSienaTest/src RSienaTest/src/model/ml
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Dec 14 14:27:35 CET 2011
Author: ripleyrm
Date: 2011-12-14 14:27:34 +0100 (Wed, 14 Dec 2011)
New Revision: 189
Added:
pkg/RSiena/po/R-RSiena.pot
pkg/RSienaTest/po/R-RSienaTest.pot
Removed:
pkg/RSiena/man/maxlikec.Rd
pkg/RSienaTest/man/maxlikec.Rd
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/print01Report.r
pkg/RSiena/R/printDataReport.r
pkg/RSiena/R/robmon.r
pkg/RSiena/R/sienaeffects.r
pkg/RSiena/R/sienaprint.r
pkg/RSiena/R/simstatsc.r
pkg/RSiena/changeLog
pkg/RSiena/data/allEffects.csv
pkg/RSiena/inst/doc/RSiena_Manual.pdf
pkg/RSiena/man/sienaModelCreate.Rd
pkg/RSiena/man/simstats0c.Rd
pkg/RSiena/src/model/Model.cpp
pkg/RSiena/src/model/ml/MLSimulation.cpp
pkg/RSiena/src/siena07models.cpp
pkg/RSienaTest/DESCRIPTION
pkg/RSienaTest/R/algorithms.r
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/print01Report.r
pkg/RSienaTest/R/printDataReport.r
pkg/RSienaTest/R/robmon.r
pkg/RSienaTest/R/sienaeffects.r
pkg/RSienaTest/R/sienaprint.r
pkg/RSienaTest/R/simstatsc.r
pkg/RSienaTest/changeLog
pkg/RSienaTest/data/allEffects.csv
pkg/RSienaTest/doc/RSienaDeveloper.tex
pkg/RSienaTest/doc/RSiena_Manual.tex
pkg/RSienaTest/inst/doc/RSiena_Manual.pdf
pkg/RSienaTest/inst/examples/runalg.r
pkg/RSienaTest/man/sienaModelCreate.Rd
pkg/RSienaTest/man/simstats0c.Rd
pkg/RSienaTest/src/model/ml/MLSimulation.cpp
pkg/RSienaTest/src/siena07models.cpp
Log:
Minor bug fixes and tidying code, added .pot files
Modified: pkg/RSiena/DESCRIPTION
===================================================================
--- pkg/RSiena/DESCRIPTION 2011-12-07 20:32:36 UTC (rev 188)
+++ pkg/RSiena/DESCRIPTION 2011-12-14 13:27:34 UTC (rev 189)
@@ -1,8 +1,8 @@
Package: RSiena
Type: Package
Title: Siena - Simulation Investigation for Empirical Network Analysis
-Version: 1.0.12.187
-Date: 2011-12-06
+Version: 1.0.12.189
+Date: 2011-12-14
Author: Various
Depends: R (>= 2.10.0)
Imports: Matrix
Modified: pkg/RSiena/R/bayes.r
===================================================================
--- pkg/RSiena/R/bayes.r 2011-12-07 20:32:36 UTC (rev 188)
+++ pkg/RSiena/R/bayes.r 2011-12-14 13:27:34 UTC (rev 189)
@@ -201,7 +201,7 @@
z$thetaMat[!z$accept, ] <<- thetaOld[!z$accept, ]
}
## print(thetaNew)
- }
+ }
## ################################
## start of function proper
## ################################
@@ -226,7 +226,7 @@
}
ctime <- proc.time()[1]
- z <- initializeBayes(data, effects, model, nbrNodes, priorSigma,
+ z <- initializeBayes(data, effects, model, nbrNodes, priorSigma,
prevAns=prevAns, clusterType=clusterType)
createStores()
@@ -270,7 +270,7 @@
tseriesplot = dev.cur()
dev.new()
tseriesratesplot = dev.cur()
- }
+ }
for (ii in 1:nwarm)
{
@@ -311,9 +311,9 @@
thetaNames<- paste(z$effects$name[!z$basicRate],
z$effects$shortName[!z$basicRate], sep=".")
rateNames <- paste(z$effects$name[basicRate],
- z$effects$shortName[basicRate],
- z$effects$period[basicRate],
- z$effects$group[basicRate], sep=".")
+ z$effects$shortName[basicRate],
+ z$effects$period[basicRate],
+ z$effects$group[basicRate], sep=".")
names(ratesdf) <- rateNames
ratesdf <- cbind(Group=thetadf[, 1, drop=FALSE], ratesdf)
names(thetadf)[-1] <- make.names(thetaNames, unique=TRUE)
@@ -424,7 +424,9 @@
z$FRAN <- getFromNamespace(model$FRANname, pkgname)
z <- z$FRAN(z, model, INIT=TRUE, data=data, effects=effects,
prevAns=prevAns)
- is.batch(TRUE)
+ z$basicRate <- z$effects$basicRate
+ z$nGroup <- z$f$nGroup
+ is.batch(TRUE)
WriteOutTheta(z)
Modified: pkg/RSiena/R/initializeFRAN.r
===================================================================
--- pkg/RSiena/R/initializeFRAN.r 2011-12-07 20:32:36 UTC (rev 188)
+++ pkg/RSiena/R/initializeFRAN.r 2011-12-14 13:27:34 UTC (rev 189)
@@ -10,10 +10,13 @@
# *****************************************************************************/
##@initializeFRAN siena07 reformat data and send to C. get targets.
initializeFRAN <- function(z, x, data, effects, prevAns=NULL, initC,
- profileData=FALSE,
- returnDeps)
+ profileData=FALSE, returnDeps=FALSE,
+ returnChains=FALSE, byGroup=FALSE,
+ returnDataFrame=FALSE, byWave=FALSE,
+ returnLoglik=FALSE, onlyLoglik=FALSE)
{
- ## fix up the interface so can call from outside robmon framework
+ z$effectsName <- deparse(substitute(effects))
+ ## fix up the interface so can call from outside robmon framework
if (is.null(z$FinDiff.method))
{
z$FinDiff.method <- FALSE
@@ -306,7 +309,7 @@
pData, lapply(f, function(x)x$dyvCovars))
ans <-.Call('ExogEvent', PACKAGE=pkgname,
pData, lapply(f, function(x)x$exog))
- ## split the names of the constraints
+ ## split the names of the constraints
higher <- attr(f, "allHigher")
disjoint <- attr(f, "allDisjoint")
atLeastOne <- attr(f, "allAtLeastOne")
@@ -479,13 +482,13 @@
CONDTARGET <- NULL
}
- simpleRates <- TRUE
- if (any(!z$effects$basicRate & z$effects$type =="rate"))
- {
+ simpleRates <- TRUE
+ if (any(!z$effects$basicRate & z$effects$type =="rate"))
+ {
## browser()
- simpleRates <- FALSE
- }
- z$simpleRates <- simpleRates
+ simpleRates <- FALSE
+ }
+ z$simpleRates <- simpleRates
ans <- .Call("setupModelOptions", PACKAGE=pkgname,
pData, pModel, MAXDEGREE, CONDVAR, CONDTARGET,
@@ -568,7 +571,7 @@
if (!initC)
{
z$f <- f
- z <- initForAlgorithms(z)
+ ##z <- initForAlgorithms(z)
z$periodNos <- attr(data, "periodNos")
z$f$myeffects <- NULL
z$f$myCompleteEffects <- NULL
@@ -586,7 +589,19 @@
z$returnDeps <- returnDeps
z$returnDepsStored <- returnDeps
z$observations <- f$observations
- z
+ z$returnChains <- returnChains
+ z$byGroup <- byGroup
+ z$byWave <- byWave
+ z$returnDataFrame <- returnDataFrame
+ z$nDependentVariables <- length(z$f$depNames)
+ if (initC)
+ {
+ NULL
+ }
+ else
+ {
+ z
+ }
}
##@createEdgeLists siena07 Reformat data for C++
createEdgeLists<- function(mat, matorig, bipartite)
Modified: pkg/RSiena/R/maxlikec.r
===================================================================
--- pkg/RSiena/R/maxlikec.r 2011-12-07 20:32:36 UTC (rev 188)
+++ pkg/RSiena/R/maxlikec.r 2011-12-14 13:27:34 UTC (rev 189)
@@ -7,51 +7,17 @@
# *
# * Description: This module contains the code for simulating the process,
# * communicating with C++. For use with maximum likelihood method, so
-# * never conditional or from finite differences, or parallel testing!
+# * never conditional or from finite differences, or parallel testing.
# *****************************************************************************/
##@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, returnLoglik=FALSE, onlyLoglik=FALSE)
+maxlikec <- function(z, x, data=NULL, effects=NULL,
+ returnChains=FALSE, byGroup=FALSE, byWave=FALSE,
+ returnDataFrame=FALSE,
+ returnLoglik=FALSE, onlyLoglik=FALSE)
{
- if (INIT || initC) ## initC is to initialise multiple C processes in
- {
- z <- initializeFRAN(z, x, data, effects, prevAns, initC,
- profileData=profileData, returnDeps=FALSE)
- z$returnDataFrame <- returnDataFrame
- z$returnChains <- returnChains
- z$byWave <- byWave
- if (initC)
- {
- return(NULL)
- }
- else
- {
- return(z)
- }
- }
- if (TERM)
- {
- z <- terminateFRAN(z, x)
- return(z)
- }
- ######################################################################
- ## iteration entry point
- ######################################################################
## retrieve stored information
f <- FRANstore()
- ##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.
if (nrow(callGrid) == 1)
{
if (byGroup)
@@ -69,17 +35,17 @@
z$returnChains, returnLoglik, onlyLoglik)
if (!onlyLoglik)
{
- ans[[6]] <- list(ans[[6]])
- ans[[7]] <- list(ans[[7]])
- if (byGroup)
- {
- ans[[8]] <- list(ans[[8]])
- ans[[9]] <- list(ans[[9]])
- ans[[10]] <- list(ans[[10]])
- }
- }
- else
- {
+ ans[[6]] <- list(ans[[6]])
+ ans[[7]] <- list(ans[[7]])
+ if (byGroup)
+ {
+ ans[[8]] <- list(ans[[8]])
+ ans[[9]] <- list(ans[[9]])
+ ans[[10]] <- list(ans[[10]])
+ }
+ }
+ else
+ {
ans[[2]] <- list(ans[[2]])
ans[[3]] <- list(ans[[3]])
ans[[4]] <- list(ans[[4]])
@@ -88,6 +54,8 @@
}
else
{
+ ## z$int2 is the number of processors if iterating by period, so 1 means
+ ## we are not. Can only parallelize by period withmaxlike.
if (z$int2 == 1)
{
anss <- apply(cbind(callGrid, 1:nrow(callGrid)),
@@ -111,27 +79,27 @@
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 (z$returnChains)
- {
- fff <- lapply(anss, function(x) x[[6]][[1]])
- fff <- split(fff, callGrid[, 1 ]) ## split by group
- ans[[6]] <- fff
- }
- ans[[7]] <- lapply(anss, "[[", 7) ## derivative
- ans[[8]] <- lapply(anss, "[[", 8)
- ans[[9]] <- lapply(anss, "[[", 9)
- ans[[10]] <- lapply(anss, "[[", 10)
- if (!byGroup)
- {
- ans[[8]] <- Reduce("+", ans[[8]]) ## accepts
- ans[[9]] <- Reduce("+", ans[[9]]) ## rejects
- ans[[10]] <- Reduce("+", ans[[10]]) ## aborts
- }
+ ans[[1]] <- sapply(anss, "[[", 1) ## statistics
+ ans[[2]] <- NULL ## scores
+ ans[[3]] <- NULL ## seeds
+ ans[[4]] <- NULL ## ntim
+ ans[[5]] <- NULL # randomseed
+ if (z$returnChains)
+ {
+ fff <- lapply(anss, function(x) x[[6]][[1]])
+ fff <- split(fff, callGrid[, 1 ]) ## split by group
+ ans[[6]] <- fff
+ }
+ ans[[7]] <- lapply(anss, "[[", 7) ## derivative
+ ans[[8]] <- lapply(anss, "[[", 8)
+ ans[[9]] <- lapply(anss, "[[", 9)
+ ans[[10]] <- lapply(anss, "[[", 10)
+ if (!byGroup)
+ {
+ ans[[8]] <- Reduce("+", ans[[8]]) ## accepts
+ ans[[9]] <- Reduce("+", ans[[9]]) ## rejects
+ ans[[10]] <- Reduce("+", ans[[10]]) ## aborts
+ }
ans[[11]] <- sapply(anss, "[[", 11)
}
else ##onlyLoglik is always byGroup (bayes)
@@ -160,17 +128,17 @@
{
fra <- -t(ans[[1]]) ##note sign change
- list(fra = fra, ntim0 = NULL, feasible = TRUE, OK = TRUE,
+ list(fra = fra, ntim0 = NULL, feasible = TRUE, OK = TRUE,
sims=NULL, dff = dff, dff2=dff2,
- chain = ans[[6]], accepts=ans[[8]],
+ chain = ans[[6]], accepts=ans[[8]],
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,
Modified: pkg/RSiena/R/phase1.r
===================================================================
--- pkg/RSiena/R/phase1.r 2011-12-07 20:32:36 UTC (rev 188)
+++ pkg/RSiena/R/phase1.r 2011-12-14 13:27:34 UTC (rev 189)
@@ -342,7 +342,7 @@
}
if (int == 1)
{
- zz <- x$FRAN(zdummy, xsmall, INIT=FALSE, fromFiniteDiff=TRUE)
+ zz <- x$FRAN(zdummy, xsmall, fromFiniteDiff=TRUE)
if (!zz$OK)
{
z$OK <- zz$OK
@@ -352,7 +352,7 @@
else
{
zz <- clusterCall(z$cl, simstats0c, zdummy, xsmall,
- INIT=FALSE, fromFiniteDiff=TRUE)
+ fromFiniteDiff=TRUE)
}
if (int == 1)
{
Modified: pkg/RSiena/R/phase2.r
===================================================================
--- pkg/RSiena/R/phase2.r 2011-12-07 20:32:36 UTC (rev 188)
+++ pkg/RSiena/R/phase2.r 2011-12-14 13:27:34 UTC (rev 189)
@@ -26,8 +26,6 @@
z$phase2fras <- array(0, dim=c(4, z$pp, 1000))
z$rejectprops <- array(0, dim=c(4, 4, 1000))
}
- ## int <- 1
- ## f <- FRANstore()
z$Phase <- 2
z$writefreq <- 1
if (!is.batch())
@@ -39,7 +37,6 @@
z$Deriv <- FALSE
z$sd <- sqrt(apply(z$sf, 2, function(x) sum(x^2) / nrow(z$sf) - mean(x)^2))
z$sd[z$fixed] <- 0
- ## browser()
Report(paste('\nPhase 2 has', x$nsub, 'subphases.\n'), cf)
z$gain <- 2 * x$firstg
if (x$nsub <= 0)
@@ -88,7 +85,7 @@
## ###############################################
## do the iterations for this repeat of this subphase
## ##############################################
- ##z <- doIterationsCopy(z, x, subphase, ...)
+ ##z <- doIterationsCopy(z, x, subphase, ...) removed as out of sync
z <- doIterations(z, x, subphase, ...)
## if (z$nit == 50) browser()
if (!z$OK || UserInterruptFlag() || UserRestartFlag() ||
@@ -218,11 +215,6 @@
{
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)
}
if ((z$nit <= 10) || (z$nit %% z$writefreq ==0))
@@ -234,8 +226,7 @@
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),
+ round((increment + val) / z$pb$pbmax * 100),
'%\n', sep = ''))
z$pb <- setProgressBar(z$pb, val + increment)
}
@@ -256,7 +247,6 @@
if (z$int == 1) ## not parallel runs at this level
{
zz <- x$FRAN(zsmall, xsmall)
- ## browser()
fra <- colSums(zz$fra) - z$targets
if (!zz$OK)
{
@@ -277,11 +267,7 @@
break
}
}
- if (x$maxlike)
- {
- # z$phase2fras[subphase, ,z$nit] <- fra
- # z$rejectprops[subphase, , z$nit] <- zz$rejectprop
- }
+ ## setup check for end of phase. either store or calculate
if (z$nit %% 2 == 1)
{
prev.fra <- fra
@@ -301,10 +287,9 @@
}
## limit change. Reporting is delayed to
## end of phase.
- ## browser()
- if (x$diag)## !maxlike at present
+ if (x$diag) ## !maxlike at present
{
- maxrat<- max(ifelse(z$sd, abs(fra)/ z$sd, 1.0))#### check this
+ maxrat<- max(ifelse(z$sd, abs(fra)/ z$sd, 1.0))
if (maxrat > x$maxmaxrat)
{
maxrat <- x$maxmaxrat / maxrat
@@ -320,7 +305,6 @@
{
fchange <- as.vector(z$gain * fra %*% z$dinv)
}
- ## browser()
fchange[z$fixed] <- 0.0
## check positivity restriction
z$positivized[z$nit, ] <- z$posj & (fchange > z$theta)
@@ -335,7 +319,6 @@
zsmall$theta <- z$theta
}
##check for user interrupt
- ## browser()
CheckBreaks()
if (UserInterruptFlag() || UserRestartFlag() || EarlyEndPhase2Flag())
{
@@ -351,446 +334,4 @@
}
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, simstats0c, 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$f$groupPeriods[x] - 1), function(y)
- vector("list", niter)))
- 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)
- 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$f$groupPeriods[ii] - 1))
- {
- 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$f$groupPeriods[ii] - 1))
- {
- 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
- groupPeriods <- attr(z$f, "groupPeriods")
- z$nDependentVariables <- length(z$f$depNames)
- ## atts <- attributes(z$f)
- 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=groupPeriods - 1, 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=groupPeriods - 1, 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/print01Report.r
===================================================================
--- pkg/RSiena/R/print01Report.r 2011-12-07 20:32:36 UTC (rev 188)
+++ pkg/RSiena/R/print01Report.r 2011-12-14 13:27:34 UTC (rev 189)
@@ -287,25 +287,49 @@
Report(c("\nFor observation moment ",
k + periodFromStart,
", number of missing values ",
- "are:\nNodes\n"),
+ "are:\n"),
sep="", outf)
- tmp <- format(cbind(1:atts$netdims[1],
- missrow, misscol))
- Report(tmp[, 1], fill=60, outf)
- Report("missing in rows\n", outf)
- Report(tmp[, 2], fill=60, outf)
- Report("missing in columns\n", outf)
- Report(tmp[, 3], fill=60, outf)
- Report(c("Total number of missing data: ",
- sum(missrow),
- ", corresponding to a fraction of ",
- format(round(sum(missrow)/atts$netdims[1] /
- (atts$netdims[1] - 1), 3), nsmall=3),
- ".\n"), sep="", outf)
- if (k > 1)
- Report(c("In reported in- and outdegrees,",
- "missings are not counted.\n"), outf)
- Report("\n", outf)
+ if (attr(depvar, "type") != "bipartite")
+ {
+ Report("Nodes\n", outf)
+ tmp <- format(cbind(1:atts$netdims[1],
+ missrow, misscol))
+ Report(tmp[, 1], fill=60, outf)
+ Report("missing in rows\n", outf)
+ Report(tmp[, 2], fill=60, outf)
+ Report("missing in columns\n", outf)
+ Report(tmp[, 3], fill=60, outf)
+ mult <- atts$netdims[1] - 1
+ }
+ else
+ {
+ Report("Senders\n", outf)
+ tmp <- format(cbind(1:atts$netdims[1],
+ missrow))
+ Report(tmp[, 1], fill=60, outf)
+ Report("missing in rows\n", outf)
+ Report(tmp[, 2], fill=60, outf)
+ tmp <- format(cbind(1:atts$netdims[2],
+ misscol))
+ Report("Receivers\n", outf)
+ Report(tmp[, 1], fill=60, outf)
+ Report("missing in columns\n", outf)
+ Report(tmp[, 2], fill=60, outf)
+ mult <- atts$netdims[2]
+ }
+ Report(c("Total number of missing data: ",
+ sum(missrow),
+ ", corresponding to a fraction of ",
+ format(round(sum(missrow)/
+ atts$netdims[1] /
+ mult, 3),
+ nsmall=3),
+ ".\n"), sep="", outf)
+
+ if (k > 1)
+ Report(c("In reported in- and outdegrees,",
+ "missings are not counted.\n"), outf)
+ Report("\n", outf)
}
else
{
Modified: pkg/RSiena/R/printDataReport.r
===================================================================
--- pkg/RSiena/R/printDataReport.r 2011-12-07 20:32:36 UTC (rev 188)
+++ pkg/RSiena/R/printDataReport.r 2011-12-14 13:27:34 UTC (rev 189)
@@ -11,6 +11,11 @@
##@DataReport siena07 Print report
DataReport <- function(z, x, f)
{
+ if (!is.null(z$effectsName))
+ {
+ Report(c("Effects object used:", z$effectsName, "\n"), outf)
+ }
+
if (z$maxlike)
{
Report(c(z$nrunMH,
@@ -32,7 +37,7 @@
observations <- attr(f, "observations") ##note this is total number of
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/rsiena -r 189
More information about the Rsiena-commits
mailing list