[Rsiena-commits] r154 - in pkg: RSiena RSiena/man RSienaTest RSienaTest/R RSienaTest/doc RSienaTest/man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Jun 5 20:17:56 CEST 2011
Author: ripleyrm
Date: 2011-06-05 20:17:56 +0200 (Sun, 05 Jun 2011)
New Revision: 154
Added:
pkg/RSienaTest/doc/bayes.tex
pkg/RSienaTest/doc/maxlikec.tex
pkg/RSienaTest/doc/missingsEtc.tex
Modified:
pkg/RSiena/DESCRIPTION
pkg/RSiena/changeLog
pkg/RSiena/man/RSiena-package.Rd
pkg/RSienaTest/DESCRIPTION
pkg/RSienaTest/R/bayes.r
pkg/RSienaTest/R/initializeFRAN.r
pkg/RSienaTest/R/maxlikec.r
pkg/RSienaTest/changeLog
pkg/RSienaTest/doc/RSIENAspec.tex
pkg/RSienaTest/doc/simstats0c.tex
pkg/RSienaTest/man/RSiena-package.Rd
Log:
New Bayes routine
Modified: pkg/RSiena/DESCRIPTION
===================================================================
--- pkg/RSiena/DESCRIPTION 2011-06-04 13:03:10 UTC (rev 153)
+++ pkg/RSiena/DESCRIPTION 2011-06-05 18:17:56 UTC (rev 154)
@@ -1,8 +1,8 @@
Package: RSiena
Type: Package
Title: Siena - Simulation Investigation for Empirical Network Analysis
-Version: 1.0.12.153
-Date: 2011-06-04
+Version: 1.0.12.154
+Date: 2011-06-05
Author: Various
Depends: R (>= 2.10.0), xtable
Imports: Matrix
Modified: pkg/RSiena/changeLog
===================================================================
--- pkg/RSiena/changeLog 2011-06-04 13:03:10 UTC (rev 153)
+++ pkg/RSiena/changeLog 2011-06-05 18:17:56 UTC (rev 154)
@@ -1,3 +1,9 @@
+2011-06-05 R-forge revision 154 RSienaTest only
+
+ * doc/missingsEtc.tex, doc/bayes.tex, doc.maxlikec.tex: added
+ * doc/RSienaSpec.tex, doc/simstats0c.tex: updated
+ * R/bayes.r, R/maxlikec.r, R/initializeFRAN.r: new Bayes routine
+
2011-06-04 R-forge revision 153
* R/effects.r, R/initializeFRAN.r, R/print07Report.r, R/sienaprint.r,
Modified: pkg/RSiena/man/RSiena-package.Rd
===================================================================
--- pkg/RSiena/man/RSiena-package.Rd 2011-06-04 13:03:10 UTC (rev 153)
+++ pkg/RSiena/man/RSiena-package.Rd 2011-06-05 18:17:56 UTC (rev 154)
@@ -30,8 +30,8 @@
\tabular{ll}{
Package: \tab RSiena\cr
Type: \tab Package\cr
-Version: \tab 1.0.12.153\cr
-Date: \tab 2011-06-04\cr
+Version: \tab 1.0.12.154\cr
+Date: \tab 2011-06-05\cr
License: \tab GPL-2 \cr
LazyLoad: \tab yes\cr
}
Modified: pkg/RSienaTest/DESCRIPTION
===================================================================
--- pkg/RSienaTest/DESCRIPTION 2011-06-04 13:03:10 UTC (rev 153)
+++ pkg/RSienaTest/DESCRIPTION 2011-06-05 18:17:56 UTC (rev 154)
@@ -1,8 +1,8 @@
Package: RSienaTest
Type: Package
Title: Siena - Simulation Investigation for Empirical Network Analysis
-Version: 1.0.12.153
-Date: 2011-06-04
+Version: 1.0.12.154
+Date: 2011-06-05
Author: Various
Depends: R (>= 2.10.0), xtable
Imports: Matrix
Modified: pkg/RSienaTest/R/bayes.r
===================================================================
--- pkg/RSienaTest/R/bayes.r 2011-06-04 13:03:10 UTC (rev 153)
+++ pkg/RSienaTest/R/bayes.r 2011-06-05 18:17:56 UTC (rev 154)
@@ -10,24 +10,24 @@
## ****************************************************************************/
##@bayes algorithm fit a Bayesian model
bayes <- function(data, effects, model, nwarm=100, nmain=100, nrunMHBatches=20,
- nrunMH=100, save=TRUE, nbrNodes=1, seed=1)
+ nrunMH=100, save=TRUE, nbrNodes=1, seed=1, dfra=NULL,
+ priorSigma=NULL)
{
createStores <- function(z)
{
npar <- length(z$theta)
basicRate <- z$basicRate
- numberRows <- nmain * nrunMHBatches *
- (z$observations - 1)
- z$posteriorTot <- rep(0, npar)
- z$posteriorMII <- matrix(0, nrow=npar, ncol=npar)
- z$lambdadist <- matrix(NA, nrow=numberRows, ncol=sum(basicRate))
- z$lambdas <- matrix(NA, nrow=numberRows, ncol=sum(basicRate))
- z$betas <- matrix(NA, nrow=numberRows, ncol=sum(!basicRate))
- z$candidates <- matrix(NA, nrow=numberRows, ncol=sum(!basicRate))
- z$acceptances <- rep(NA, numberRows)
- z$MHacceptances <- matrix(NA, nrow=numberRows, ncol=7)
- z$MHrejections <- matrix(NA, nrow=numberRows , ncol=7)
- z$MHproportions <- matrix(NA, nrow=numberRows, ncol=7)
+ numberRows <- nmain * nrunMHBatches
+ z$posteriorTot <- matrix(0, nrow=z$nGroup, ncol=npar)
+ z$posteriorMII <- array(0, dim=c(z$nGroup, npar, npar))
+ #z$lambdadist <- matrix(NA, nrow=numberRows, ncol=sum(basicRate))
+ #z$lambdas <- matrix(NA, nrow=numberRows, ncol=sum(basicRate))
+ #z$betas <- matrix(NA, nrow=numberRows, ncol=sum(!basicRate))
+ 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(z$nGroup, numberRows, 7))
+ z$MHrejections <- array(NA, dim=c(z$nGroup, numberRows, 7))
+ z$MHproportions <- array(NA, dim=c(z$nGroup, numberRows, 7))
z
}
storeData <- function()
@@ -37,15 +37,15 @@
basicRate <- z$basicRate
npar <- length(z$theta)
end <- start + nrun - 1
- z$accepts <- as.logical(z$accepts)
z$acceptances[start:end] <- z$accepts
- z$lambdas[start:end, ] <- z$parameters[, basicRate]
- z$lambdadist[start:end, ] <- z$shapes[, basicRate]
- z$candidates[start:end, ] <- z$parameters[, !basicRate]
- z$parameters[!z$accepts, !basicRate] <- NA
- z$parameters[, !basicRate] <-
- carryForward(z$parameters[, !basicRate], z$betas[1:end,])
- z$betas[start:end, ] <- z$parameters[, !basicRate]
+ # z$lambdas[start:end, ] <- z$parameters[, basicRate]
+ # z$lambdadist[start:end, ] <- z$shapes[, basicRate]
+ z$candidates[start:end, ] <- z$parameters
+ # z$parameters[start:end, ] <- NA
+ # z$parameters[, !basicRate] <-
+ # carryForward(z$parameters[, !basicRate], z$betas[1:end,])
+ # z$betas[start:end, ] <- z$parameters[, !basicRate]
+ browser()
z$posteriorTot <- z$posteriorTot + colSums(z$parameters)
for (i in npar)
{
@@ -120,12 +120,12 @@
z <- MCMCcycle(z, nrunMH=1, nrunMHBatches=100)
actual <- z$BayesAcceptances ## acceptances
ans <- rescaleCGD(100 * (z$observations - 1))
- z$scaleFactor <- z$scaleFactor * ans
+ z$scaleFactors <- z$scaleFactors * ans
if (success | iter == maxiter)
{
break
}
- if (z$scaleFactor < tiny)
+ if (z$scaleFactors < tiny)
{
cat('scalefactor < tiny\n')
browser()
@@ -141,11 +141,20 @@
## Should we save and restore theta if we use multiple processors?
## Currently get separate thetas per wave (and then mix them up).
- z <- initializeBayes(data, effects, model, nbrNodes, seed)
+ z <- initializeBayes(data, effects, model, nbrNodes, seed, priorSigma)
z <- createStores(z)
z$sub <- 0
+ if (is.null(dfra))
+ {
+ z <- getDFRA(z, 10)
+ }
+ else
+ {
+ ## maybe some validation here one day
+ z$dfra <- dfra
+ }
z <- improveMH(z)
if (!save)
@@ -214,73 +223,90 @@
}))
}
}
+ z$FRAN <- NULL
z
}
MCMCcycle <- function(z, nrunMH, nrunMHBatches)
{
- f <- FRANstore() ## retrieve info
-
- ## set up a matrix of required group/periods
- groupPeriods <- attr(f, "groupPeriods")
- callGrid <- cbind(rep(1:f$nGroup, groupPeriods - 1),
- as.vector(unlist(sapply(groupPeriods - 1,
- function(x) 1:x))))
-
- ## 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.
- ## browser()
- if (nrow(callGrid) == 1)
+ z$accepts <- matrix(NA, nrow=z$nGroup, nrunMHBatches)
+ z$parameters <- array(NA, dim=c(nrunMHBatches, z$nGroup, z$pp))
+ for (i in 1:nrunMHBatches)
{
- ans <- .Call("MCMCcycle", PACKAGE=pkgname, f$pData, f$pModel,
- f$myeffects, as.integer(1), as.integer(1),
- z$scaleFactor, nrunMH, nrunMHBatches)
- }
- else
- {
- if (z$int2 == 1)
+ for (j in 1:nrunMH)
{
- anss <- apply(callGrid, 1, doMCMCcycle, z$scaleFactor,
- nrunMH, nrunMHBatches)
+ ans <- z$FRAN(z, returnChains=TRUE, byGroup=TRUE)
}
- else
- {
- use <- 1:(min(nrow(callGrid), z$int2))
- anss <- parRapply(z$cl[use], callGrid, doMCMCcycle, z$scaleFactor,
- nrunMH, nrunMHBatches)
- }
- ## reorganize the anss so it looks like the normal one
- ans <- NULL
- ans[[1]] <- c(sapply(anss, "[[", 1))## acceptances
- ans[[2]] <- do.call(rbind, lapply(anss, "[[", 2))
- ans[[3]] <- do.call(rbind, lapply(anss, "[[", 3))
- ans[[4]] <- do.call(rbind, lapply(anss, "[[", 4))
- ans[[5]] <- do.call(rbind, lapply(anss, "[[", 5))
- ans[[6]] <- do.call(c, lapply(anss, "[[", 6))
- }
- ## process the return values
- z$BayesAcceptances <- sum(ans[[1]])
- z$accepts <- ans[[1]]
- z$parameters <- ans[[2]]
- z$shapes <- ans[[3]]
- z$MHaccepts <- ans[[4]]
- z$MHrejects <- ans[[5]]
- z$chains <- ans[[6]]
+ z$chain <- ans$chain
+ z <- sampleParameters(z)
+ z$accepts[, i] <- z$accept
+ z$parameters[i, , ] <- z$thetaMat
+ }
+ z$BayesAcceptances <- sum(z$accepts)
z
}
-
-doMCMCcycle <- function(x, scaleFactor, nrunMH, nrunMHBatches)
+sampleParameters <- function(z)
{
- group <- x[1]
- period <- x[2]
+ #browser()
f <- FRANstore()
- .Call("MCMCcycle", PACKAGE=pkgname, f$pData, f$pModel,
- f$myeffects, as.integer(period),
- as.integer(group),
- scaleFactor, nrunMH, nrunMHBatches)
+ ## get a multivariate normal with covariance matrix dfra multiplied by a
+ ## scale factor which varies between groups
+ require(MASS)
+ require(mvtnorm)
+ thetaChanges <- t(sapply(z$scaleFactors, function(x)
+ mvrnorm(1, mu=rep(0, z$pp),
+ Sigma= x*z$dfra)))
+ priorOld <- apply(z$thetaMat, 1, dmvnorm, mean=rep(0, z$pp),
+ sigma=z$priorSigma,
+ log=TRUE)
+ priorNew<- apply(z$thetaMat + thetaChanges, 1, dmvnorm, mean=rep(0, z$pp),
+ sigma=z$priorSigma, log=TRUE)
+ logpOld <- sapply(z$chain, function(x)
+ {
+ sum(sapply(x, function(x1) sum(x1$LogOptionSetProb +
+ x1$LogChoiceProb)))
+ }
+ )
+ thetas <- z$thetaMat + thetaChanges
+ ans <- lapply(1:length(z$chain), function(i)
+ {
+ chaini <- z$chain[[i]]
+ lapply(1:length(chaini), function(i1)
+ {
+ ch <- chaini[[i1]]
+ .Call("getChainProbabilities", ch,
+ PACKAGE = pkgname,
+ f$pData, f$pModel,
+ as.integer(i),
+ as.integer(i1),
+ f$myeffects, thetas[i, ])
+ }
+ )
+ }
+ )
+ logpNew <- sapply(ans, function(x)
+ {
+ sum(sapply(x, function(x1) sum(x1$LogOptionSetProb +
+ x1$LogChoiceProb)))
+ }
+ )
+ proposalProbability <- priorNew - priorOld + logpNew - logpOld
+ if (log(runif(1)) < proposalProbability)
+ {
+ z$accept <- TRUE
+ z$thetaMat <- thetas
+ }
+ else
+ {
+ z$accept <- FALSE
+ }
+ cat(thetas, priorNew, priorOld, logpNew, logpOld, exp(proposalProbability),
+ z$accept,'\n')
+ z
}
-initializeBayes <- function(data, effects, model, nbrNodes, seed)
+
+initializeBayes <- function(data, effects, model, nbrNodes, seed, priorSigma)
{
## initialise
set.seed(seed)
@@ -298,9 +324,9 @@
{
set.seed(model$randomSeed)
}
- model$FRAN <- getFromNamespace(model$FRANname, pos=grep("RSiena",
+ z$FRAN <- getFromNamespace(model$FRANname, pos=grep("RSiena",
search())[1])
- z <- model$FRAN(z, model, INIT=TRUE, data=data, effects=effects)
+ z <- z$FRAN(z, model, INIT=TRUE, data=data, effects=effects)
is.batch(TRUE)
@@ -324,7 +350,32 @@
z$numm <- 20
z$scaleFactor <- 1
+ z$scaleFactors <- rep(1, z$nGroup)
+ z$returnDataFrame <- TRUE # chains come back as data frames not lists
+ if (is.null(priorSigma))
+ {
+ z$priorSigma <- diag(z$pp) * 10000
+ }
+ else
+ {
+ z$priorSigma <- priorSigma
+ }
+ z
+}
-
+getDFRA <- function(z, n)
+{
+ # do n MLmodelsteps with the initial thetas and get
+ # derivs
+ z$sdf <- array(0, dim=c(n, z$pp, z$pp))
+ z$Phase <- 1
+ z$Deriv <- TRUE
+ for (i in 1:n)
+ {
+ ans <- z$FRAN(z)
+ z$sdf[i, , ] <- ans$dff
+ }
+ z$dfra <- t(apply(z$sdf, c(2, 3), mean))
z
}
+
Modified: pkg/RSienaTest/R/initializeFRAN.r
===================================================================
--- pkg/RSienaTest/R/initializeFRAN.r 2011-06-04 13:03:10 UTC (rev 153)
+++ pkg/RSienaTest/R/initializeFRAN.r 2011-06-05 18:17:56 UTC (rev 154)
@@ -455,6 +455,8 @@
z$maxlikeTargets2 <- ans
z$mult <- x$mult
z$nrunMH <- z$mult * sum(z$maxlikeTargets[z$effects$basicRate])
+ ##thetaMat is to allow different thetas for each group in Bayes
+ z$thetaMat <- matrix(z$theta, nrow=nGroup, ncol=z$pp, byrow=TRUE)
}
}
@@ -488,7 +490,7 @@
if (x$maxlike)
{
- simpleRates <- TRUE ## get a sensible test for this soon!
+ simpleRates <- TRUE
if (any(!z$effects$basicRate & z$effects$type =="rate"))
{
# browser()
@@ -574,14 +576,15 @@
z$f <- f
z <- initForAlgorithms(z)
z$periodNos <- attr(data, "periodNos")
- if (! returnDeps) {
+ z$f$myeffects <- NULL
+ z$f$myCompleteEffects <- NULL
+ if (!returnDeps) {
z$f[1:nGroup] <- NULL
}
}
if (initC || (z$int == 1 && z$int2 == 1 &&
(is.null(z$nbrNodes) || z$nbrNodes == 1)))
{
-
f[1:nGroup] <- NULL
}
FRANstore(f) ## store f in FRANstore
Modified: pkg/RSienaTest/R/maxlikec.r
===================================================================
--- pkg/RSienaTest/R/maxlikec.r 2011-06-04 13:03:10 UTC (rev 153)
+++ pkg/RSienaTest/R/maxlikec.r 2011-06-05 18:17:56 UTC (rev 154)
@@ -12,7 +12,7 @@
##@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)
+ returnDeps=FALSE, returnChains=FALSE, byGroup=FALSE)
{
if (INIT || initC) ## initC is to initialise multiple C processes in phase3
{
@@ -68,7 +68,7 @@
{
if (z$int2 == 1)
{
- anss <- apply(callGrid, 1, doMLModel, z$Deriv, z$theta,
+ anss <- apply(callGrid, 1, doMLModel, z$Deriv, z$thetaMat,
returnDeps, z$nrunMH, z$addChainToStore,
z$needChangeContributions, z$returnDataFrame,
returnChains)
@@ -76,7 +76,8 @@
else
{
use <- 1:(min(nrow(callGrid), z$int2))
- anss <- parRapply(z$cl[use], callGrid, doMLModel, z$Deriv, z$theta,
+ anss <- parRapply(z$cl[use], callGrid, doMLModel, z$Deriv,
+ z$thetaMat,
returnDeps, z$nrunMH, z$addChainToStore,
z$needChangeContributions,
z$returnDataFrame, returnChains)
@@ -108,30 +109,33 @@
## browser()
if (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
}
ans[[7]] <- lapply(anss, "[[", 7) ## derivative
ans[[8]] <- lapply(anss, "[[", 8)
- ans[[8]] <- do.call("+", ans[[8]]) ## accepts
ans[[9]] <- lapply(anss, "[[", 9)
- ans[[9]] <- do.call("+", ans[[9]]) ## rejects
ans[[10]] <- lapply(anss, "[[", 10)
- ans[[10]] <- do.call("+", ans[[10]]) ## aborts
+ if (!byGroup)
+ {
+ ans[[8]] <- Reduce("+", ans[[8]]) ## accepts
+ ans[[9]] <- Reduce("+", ans[[9]]) ## rejects
+ ans[[10]] <- Reduce("+", ans[[10]]) ## aborts
+ }
}
- dff <- ans[[7]]
fra <- -t(ans[[1]]) ##note sign change
FRANstore(f)
- ## if (returnDeps)
- ## sims <- ans[[6]]
- ## else
- sims <- NULL
+ ## if (returnDeps)
+ ## sims <- ans[[6]]
+ ## else
+ sims <- NULL
- #print(length(sims[[1]]))
+ ##print(length(sims[[1]]))
##if (returnDeps) ## this is not finished yet!
##{
## ## attach the names
@@ -144,93 +148,92 @@
## {
## 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
{
- ## current format is a list of vectors of the lower? triangle
- ## 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)
- nPeriods <- length(ans[[7]])
- dff2 <- array(0, dim=c(nPeriods, z$pp, z$pp))
- for (period in 1:nPeriods)
- {
- dffraw <- ans[[7]][[period]]
- ## start indexes row/col for first effect for the variable
- start <- 1
- ## rawsub is subscript in the vector
- rawsub <- 1
- ## f$myeffects is a list of data frames per dependent variable
- ## of selected effects.
- for (i in 1:length(f$myeffects))
- {
- dffPeriod <- matrix(0, nrow=z$pp, ncol=z$pp)
- ## rows is the number of effects for this variable
- rows <- nrow(f$myeffects[[i]])
- ## nRates is the number of rate effects for this variable
- ## at present nRates will be the number of periods.
- nRates <- sum(f$myeffects[[i]]$type == 'rate')
- ## nonRates is the number of rows/cols in the objective function
- ## part of the derivative matrix dff.
- nonRates <- rows - nRates
- ## first put the basic rate for this variable in the right place
- dffPeriod[cbind(start:(start + nRates -1),
- start:(start + nRates - 1))] <-
- dffraw[rawsub:(rawsub + nRates - 1)]
- ##
- ## now the matrix of objective function effects
- ##
- start <- start + nRates
- ##
- rawsub <- rawsub + nRates
- ##
- dffPeriod[start : (start + nonRates - 1 ),
- start : (start + nonRates - 1)] <-
- dffraw[rawsub:(rawsub + nonRates * nonRates - 1)]
- start <- start + nonRates
- rawsub <- rawsub + nonRates * nonRates
- dffPeriod <- dffPeriod + t(dffPeriod)
- diag(dffPeriod) <- diag(dffPeriod) / 2
- dff2[period , , ] <- dff2[period, , ] - dffPeriod
- dff <- dff - dffPeriod
- }
- }
+ resp <- reformatDerivs(z, f, ans[[7]])
+ dff <- resp[[1]]
+ dff2 <- resp[[2]]
}
else
{
dff <- NULL
dff2 <- NULL
}
- # browser()
-#print(length(ans[[6]][[1]][[1]]))
+ ## browser()
+ ##print(length(ans[[6]][[1]][[1]]))
list(fra = fra, ntim0 = NULL, feasible = TRUE, OK = TRUE,
sims=sims, dff = dff, dff2=dff2,
chain = ans[[6]], accepts=ans[[8]],
rejects= ans[[9]], aborts=ans[[10]])
}
-dist2full<-function(dis) {
-
- n<-attr(dis,"Size")
-
- full<-matrix(0,n,n)
-
- full[lower.tri(full)]<-dis
-
- full+t(full)
-
-}
##@doMLModel Maximum likelihood
-doMLModel <- function(x, Deriv, theta, returnDeps, nrunMH, addChainToStore,
+doMLModel <- function(x, Deriv, thetaMat, returnDeps, nrunMH, addChainToStore,
needChangeContributions, returnDataFrame, returnChains)
{
f <- FRANstore()
.Call("mlPeriod", PACKAGE=pkgname, Deriv, f$pData,
- f$pModel, f$myeffects, theta, returnDeps,
+ f$pModel, f$myeffects, thetaMat[x[1], ], returnDeps,
as.integer(x[1]), as.integer(x[2]), nrunMH, addChainToStore,
needChangeContributions, returnDataFrame, returnChains)
}
+
+reformatDerivs <- function(z, f, derivList)
+{
+ ## current format is a list of vectors of the lower? triangle
+ ## 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)
+ nPeriods <- length(derivList)
+ dff2 <- array(0, dim=c(nPeriods, z$pp, z$pp))
+ for (period in 1:nPeriods)
+ {
+ dffraw <- derivList[[period]]
+ ## start indexes row/col for first effect for the variable
+ start <- 1
+ ## rawsub is subscript in the vector
+ rawsub <- 1
+ ## f$myeffects is a list of data frames per dependent variable
+ ## of selected effects.
+ for (i in 1:length(f$myeffects))
+ {
+ dffPeriod <- matrix(0, nrow=z$pp, ncol=z$pp)
+ ## rows is the number of effects for this variable
+ rows <- nrow(f$myeffects[[i]])
+ ## nRates is the number of rate effects for this variable
+ ## at present nRates will be the number of periods.
+ nRates <- sum(f$myeffects[[i]]$type == 'rate')
+ ## nonRates is the number of rows/cols in the objective function
+ ## part of the derivative matrix dff.
+ nonRates <- rows - nRates
+ ## first put the basic rate for this variable in the right place
+ dffPeriod[cbind(start:(start + nRates -1),
+ start:(start + nRates - 1))] <-
+ dffraw[rawsub:(rawsub + nRates - 1)]
+ ##
+ ## now the matrix of objective function effects
+ ##
+ start <- start + nRates
+ ##
+ rawsub <- rawsub + nRates
+ ##
+ dffPeriod[start : (start + nonRates - 1 ),
+ start : (start + nonRates - 1)] <-
+ dffraw[rawsub:(rawsub + nonRates * nonRates - 1)]
+ start <- start + nonRates
+ rawsub <- rawsub + nonRates * nonRates
+ dffPeriod <- dffPeriod + t(dffPeriod)
+ diag(dffPeriod) <- diag(dffPeriod) / 2
+ dff2[period , , ] <- dff2[period, , ] - dffPeriod
+ dff <- dff - dffPeriod
+ }
+ }
+ list(dff, dff2)
+}
+
Modified: pkg/RSienaTest/changeLog
===================================================================
--- pkg/RSienaTest/changeLog 2011-06-04 13:03:10 UTC (rev 153)
+++ pkg/RSienaTest/changeLog 2011-06-05 18:17:56 UTC (rev 154)
@@ -1,3 +1,9 @@
+2011-06-05 R-forge revision 154 RSienaTest only
+
+ * doc/missingsEtc.tex, doc/bayes.tex, doc.maxlikec.tex: added
+ * doc/RSienaSpec.tex, doc/simstats0c.tex: updated
+ * R/bayes.r, R/maxlikec.r, R/initializeFRAN.r: new Bayes routine
+
2011-06-04 R-forge revision 153
* R/effects.r, R/initializeFRAN.r, R/print07Report.r, R/sienaprint.r,
Modified: pkg/RSienaTest/doc/RSIENAspec.tex
===================================================================
--- pkg/RSienaTest/doc/RSIENAspec.tex 2011-06-04 13:03:10 UTC (rev 153)
+++ pkg/RSienaTest/doc/RSIENAspec.tex 2011-06-05 18:17:56 UTC (rev 154)
@@ -310,7 +310,8 @@
\STATE Check that the ends of each interval in each object are not less than 1
or greater than the number of waves and that each line has an even number of
digits.
-\STATE Generate a data frame of events(section \ref{sec:events}), a matrix of active flags(section \ref{sec:active}) and a matrix
+\STATE Generate a data frame of events(section \ref{sec:events}),
+a matrix of activeStart flags(section \ref{sec:active}) and a matrix
of actions(section \ref{sec:actions})\\
\STATE Add these to the object as attributes
\ENDFOR
@@ -325,11 +326,11 @@
\sfn{ }\nnm{actor}\\
\sfn{ }\nnm{time} between 0 and 1
\end{algorithmic}
-\subsubsection{ActiveFlags}
+\subsubsection{ActiveStart Flags}
\label{sec:active}
\begin{algorithmic}
-\STATE Active flags matrix has a row per actor and a column per period\\
-\sfn{ }TRUE if the actor is present at the start of the period
+\STATE activeStart matrix has a row per actor and a column per period\\
+\sfn{ }TRUE if the actor is active at the start of the period
\end{algorithmic}
\subsubsection{Action}
\label{sec:actions}
Added: pkg/RSienaTest/doc/bayes.tex
===================================================================
--- pkg/RSienaTest/doc/bayes.tex (rev 0)
+++ pkg/RSienaTest/doc/bayes.tex 2011-06-05 18:17:56 UTC (rev 154)
@@ -0,0 +1,103 @@
+\documentclass[12pt,a4paper]{article}
+\usepackage[pdftex,dvipsnames]{color}
+\usepackage{graphicx,times}
+\usepackage{amsmath}
+\usepackage{amsfonts}
+\usepackage{amssymb}
+\usepackage{algorithm}
+\usepackage{algorithmic}
+\usepackage[pdfstartview={}]{hyperref}
+\textheight=9.5in
+\textwidth=6.0in
+\topmargin=-0.5in
+\evensidemargin=0in
+\oddsidemargin=0in
+\newcommand{\eps}{\varepsilon}
+\newcommand{\abso}[1]{\;\mid#1\mid\;}
+\renewcommand{\=}{\,=\,}
+\newcommand{\+}{\,+\,}
+% ----------------------------------------------------------------
+\newcommand{\remark}[1]{\par\noindent{\color[named]{ProcessBlue}#1}\par}
+\newcommand{\mcc}[2]{\multicolumn{#1}{c}{#2}}
+\newcommand{\mcp}[2]{\multicolumn{#1}{c|}{#2}}
+\newcommand{\nm}[1]{\textsf{\small #1}}
+\newcommand{\nnm}[1]{\textsf{\small\textit{#1}}}
+\newcommand{\nmm}[1]{\nnm{#1}}
+\newcommand{\R}{{\sf R }}
+\newcommand{\sfn}[1]{\textbf{\texttt{#1}}}
+\newcommand{\Rn}{{\sf R}}
+\newcommand{\rs}{{\sf RSiena}}
+\newcommand{\RS}{{\sf RSiena }}
+\newcommand{\SI}{{\sf Siena3 }}
+\newcommand{\Sn}{{\sf Siena3}}
+% no labels in list of references:
+\makeatletter
+\renewcommand\@biblabel{}
+\makeatother
+
+\hyphenation{Snij-ders Duijn DataSpecification dataspecification dependentvariable ModelSpecification}
+
+% centered section headings with a period after the number;
+% sans serif fonts for section and subsection headings
+\renewcommand{\thesection}{\arabic{section}.}
+\renewcommand{\thesubsection}{\thesection\arabic{subsection}}
+\makeatletter
+ \renewcommand{\section}{\@startsection{section}{1}
+ {0pt}{\baselineskip}{0.5\baselineskip}
+ {\centering\sffamily} }
+ \renewcommand{\subsection}{\@startsection{subsection}{2}
+ {0pt}{0.7\baselineskip}{0.3\baselineskip}
+ {\sffamily} }
+\makeatother
+
+\newcommand{\ts}[1]{\par{\color[named]{Red}TS: #1}\par}
+
+\renewcommand{\baselinestretch}{1.0} %% For line spacing.
+\setlength{\parindent}{0pt}
+\setlength{\parskip}{1ex plus1ex}
+\begin{document}
+
+\title{RSiena: Bayesian models}
+\author{Ruth Ripley}
+\date{}
+\maketitle
+
+\centerline{\emph{\today}}
+\bigskip
+\begin{enumerate}
+\item Input of a covariance matrix for the prior is mandatory: default the
+ identity times 10000.
+\item Input of a covariance matrix for the Bayesian proposal is optional: if
+ null it will be calculated from a few iterations of the MH step.
+\item Default number of iterations for this step if performed is 10.
+\item Could allow clever extraction of dfra from a previous result of siena07 or
+ bayes. Maybe this will work anyway via prevAns.
+\end{enumerate}
+\paragraph{Main algorithm}
+\begin{algorithmic}
+\STATE set up data in C as usual
+\STATE create minimal chain and do burnin
+\STATE Do $n$ groups of \sfn{nrunMH} steps to estimate \nnm{dfra} as in siena07
+\STATE Get scalefactors such that about 40 out of 100 Bayesian proposals after
+single MH steps are accepted. See below for details of
+generation and probabilities for Bayesian proposals.
+
+\end{algorithmic}
+\paragraph{Bayesian proposals}
+\label{sec:prop}
+\begin{algorithmic}
+\FORALL {groups}
+\STATE get a multivariate normal with mean 0 and covariance \nnm{dfra} *
+scale factor for this group
+\STATE Calculate the proposal probability:
+\begin{description}
+\item[prior] Multivariate normal density for the parameters with mean 0
+ and covariance as supplied in the input argument.
+\item[chain] Calculate sum of log probabilities of choice of variable/actor
+plus log choice probabilities.
+\end{description}
+\STATE The log probability of acceptance is then new - old of log prior +
+log chain
+\ENDFOR
+\end{algorithmic}
+\end{document}
Property changes on: pkg/RSienaTest/doc/bayes.tex
___________________________________________________________________
Added: svn:eol-style
+ native
Added: pkg/RSienaTest/doc/maxlikec.tex
===================================================================
--- pkg/RSienaTest/doc/maxlikec.tex (rev 0)
+++ pkg/RSienaTest/doc/maxlikec.tex 2011-06-05 18:17:56 UTC (rev 154)
@@ -0,0 +1,566 @@
+\documentclass[12pt,a4paper]{article}
+\usepackage[pdftex,dvipsnames]{color}
+\usepackage{graphicx,times}
+\usepackage{amsmath}
+\usepackage{amsfonts}
+\usepackage{amssymb}
+\usepackage{algorithm}
+\usepackage{algorithmic}
+\usepackage[pdfstartview={}]{hyperref}
+\textheight=9.5in
+\topmargin=-0.5in
+\newcommand{\eps}{\varepsilon}
+\newcommand{\abso}[1]{\;\mid#1\mid\;}
+\renewcommand{\=}{\,=\,}
+\newcommand{\+}{\,+\,}
+% ----------------------------------------------------------------
+\newcommand{\remark}[1]{\par\noindent{\color[named]{ProcessBlue}#1}\par}
+\newcommand{\mcc}[2]{\multicolumn{#1}{c}{#2}}
+\newcommand{\mcp}[2]{\multicolumn{#1}{c|}{#2}}
+\newcommand{\nm}[1]{\textsf{\small #1}}
+\newcommand{\nnm}[1]{\textsf{\small\textit{#1}}}
+\newcommand{\nmm}[1]{\nnm{#1}}
+\newcommand{\R}{{\sf R }}
+\newcommand{\sfn}[1]{\textbf{\texttt{#1}}}
+\newcommand{\Rn}{{\sf R}}
+\newcommand{\rs}{{\sf RSiena}}
+\newcommand{\RS}{{\sf RSiena }}
+\newcommand{\SI}{{\sf Siena3 }}
+\newcommand{\Sn}{{\sf Siena3}}
+% no labels in list of references:
+\makeatletter
+\renewcommand\@biblabel{}
+\makeatother
+
+\hyphenation{Snij-ders Duijn DataSpecification dataspecification dependentvariable ModelSpecification}
+
+% centered section headings with a period after the number;
+% sans serif fonts for section and subsection headings
+\renewcommand{\thesection}{\arabic{section}.}
+\renewcommand{\thesubsection}{\thesection\arabic{subsection}}
+\makeatletter
+ \renewcommand{\section}{\@startsection{section}{1}
+ {0pt}{\baselineskip}{0.5\baselineskip}
+ {\centering\sffamily} }
+ \renewcommand{\subsection}{\@startsection{subsection}{2}
+ {0pt}{0.7\baselineskip}{0.3\baselineskip}
+ {\sffamily} }
+\makeatother
+
+\newcommand{\ts}[1]{\par{\color[named]{Red}TS: #1}\par}
+
+\renewcommand{\baselinestretch}{1.0} %% For line spacing.
+\setlength{\parindent}{0pt}
+\setlength{\parskip}{1ex plus1ex}
+\raggedright
+\begin{document}
+
+\title{Maximum likelihood/FRAN}
+\author{Ruth Ripley}
+\date{}
+\maketitle
+
+\centerline{\emph{\today}}
+\bigskip
+\section{Introduction}
+
+The \R function \nnm{maxlikec} (so named during development when I had other
+versions which did not call C!) is the function which controls the maximum
+likleihood simulations. It is used as an argument to the function
+\nnm{sienaModelCreate} and its name is stored as the element named \nm{FRAN} of
+the model object. It has three (or 3+) modes: initial, to set up the C++ data
+structure, terminal (to tidy up) and ordinary (to do one complete
+simulation). (There is an extra mode for using multiple processors, which does
+some of the work of an initial call.) It uses C++ functions for most of its
+work.
+
+The interface to C++ and to \nnm{robmon} is more or less the same as for
+\nnm{simstats0c}. The terminal call is the same as for \nnm{simstats0c}.
+
+
+\section{Initial call to maxlikec}
+As for \nnm{simstatsc}, plus an additional call to C++ to initialise the
+chains.
+
+\begin{enumerate}
+\item Copy over the parameters needed for maxmimum liklihood:
+\begin{enumerate}
+\item Values controlling length of permutation
+\begin{description}
+\item[maximumPermutationLength]
+\item[minimumPermutationLength]
+\item[initialPermutationLength]
+\end{description}
+\item Probabilities of the different MH steps
+\begin{description}
+\item[insertDiagonalProbability]
+\item[cancelDiagonalProbability]
+\item[permuteProbability]
+\item[insertPermuteProbability]
+\item[deletePermuteProbability]
+\item[insertRandomMissingProbability]
+\item[deleteRandomMissingProbability]
+\end{description}
+\item Proportion of missing data
+\begin{description}
+\item[prmin]
+\item[prmib]
+\end{description}
+\end{enumerate}
+\item
+\begin{algorithmic}
+\IF{this is the first call, (to main process)}
+\STATE Set up a minimal chain: see section \ref{sec:minch}.
+\STATE Do a pre burnin for each chain: see section \ref{sec:preburnin}
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/rsiena -r 154
More information about the Rsiena-commits
mailing list