[Rsiena-commits] r156 - in pkg: RSiena RSiena/R RSiena/inst/doc RSiena/man RSiena/src RSiena/src/model RSiena/src/model/ml RSiena/src/model/variables RSienaTest RSienaTest/R RSienaTest/doc RSienaTest/inst/doc RSienaTest/man RSienaTest/src RSienaTest/src/model RSienaTest/src/model/ml RSienaTest/src/model/variables
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Jun 13 16:32:33 CEST 2011
Author: ripleyrm
Date: 2011-06-13 16:32:32 +0200 (Mon, 13 Jun 2011)
New Revision: 156
Modified:
pkg/RSiena/DESCRIPTION
pkg/RSiena/NAMESPACE
pkg/RSiena/R/bayes.r
pkg/RSiena/R/effects.r
pkg/RSiena/R/initializeFRAN.r
pkg/RSiena/R/maxlikec.r
pkg/RSiena/R/phase1.r
pkg/RSiena/R/sienaprint.r
pkg/RSiena/changeLog
pkg/RSiena/inst/doc/RSiena_Manual.pdf
pkg/RSiena/man/RSiena-package.Rd
pkg/RSiena/src/model/Model.cpp
pkg/RSiena/src/model/Model.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/MLSimulation.h
pkg/RSiena/src/model/variables/DependentVariable.cpp
pkg/RSiena/src/model/variables/DependentVariable.h
pkg/RSiena/src/siena07internals.cpp
pkg/RSiena/src/siena07internals.h
pkg/RSiena/src/siena07models.cpp
pkg/RSiena/src/siena07models.h
pkg/RSiena/src/siena07setup.cpp
pkg/RSiena/src/siena07utilities.cpp
pkg/RSiena/src/siena07utilities.h
pkg/RSienaTest/DESCRIPTION
pkg/RSienaTest/NAMESPACE
pkg/RSienaTest/R/bayes.r
pkg/RSienaTest/R/effects.r
pkg/RSienaTest/R/initializeFRAN.r
pkg/RSienaTest/R/maxlikec.r
pkg/RSienaTest/R/phase1.r
pkg/RSienaTest/R/sienaprint.r
pkg/RSienaTest/changeLog
pkg/RSienaTest/doc/RSiena_Manual.tex
pkg/RSienaTest/doc/bayes.tex
pkg/RSienaTest/inst/doc/RSiena_Manual.pdf
pkg/RSienaTest/man/RSiena-package.Rd
pkg/RSienaTest/src/model/Model.cpp
pkg/RSienaTest/src/model/Model.h
pkg/RSienaTest/src/model/ml/Chain.cpp
pkg/RSienaTest/src/model/ml/Chain.h
pkg/RSienaTest/src/model/ml/MLSimulation.cpp
pkg/RSienaTest/src/model/ml/MLSimulation.h
pkg/RSienaTest/src/model/variables/DependentVariable.cpp
pkg/RSienaTest/src/model/variables/DependentVariable.h
pkg/RSienaTest/src/siena07internals.cpp
pkg/RSienaTest/src/siena07internals.h
pkg/RSienaTest/src/siena07models.cpp
pkg/RSienaTest/src/siena07models.h
pkg/RSienaTest/src/siena07setup.cpp
pkg/RSienaTest/src/siena07utilities.cpp
pkg/RSienaTest/src/siena07utilities.h
Log:
Bayesian module
Modified: pkg/RSiena/DESCRIPTION
===================================================================
--- pkg/RSiena/DESCRIPTION 2011-06-12 16:27:55 UTC (rev 155)
+++ pkg/RSiena/DESCRIPTION 2011-06-13 14:32:32 UTC (rev 156)
@@ -1,12 +1,12 @@
Package: RSiena
Type: Package
Title: Siena - Simulation Investigation for Empirical Network Analysis
-Version: 1.0.12.155
-Date: 2011-06-11
+Version: 1.0.12.156
+Date: 2011-06-13
Author: Various
Depends: R (>= 2.10.0), xtable
Imports: Matrix
-Suggests: tcltk, snow, rlecuyer, network, codetools, lattice
+Suggests: tcltk, snow, rlecuyer, network, codetools, lattice, MASS
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-06-12 16:27:55 UTC (rev 155)
+++ pkg/RSiena/NAMESPACE 2011-06-13 14:32:32 UTC (rev 156)
@@ -42,3 +42,4 @@
S3method(summary, sienaEffects)
S3method(print, summary.sienaEffects)
S3method(edit, sienaEffects)
+S3method(print, chains.data.frame)
Modified: pkg/RSiena/R/bayes.r
===================================================================
--- pkg/RSiena/R/bayes.r 2011-06-12 16:27:55 UTC (rev 155)
+++ pkg/RSiena/R/bayes.r 2011-06-13 14:32:32 UTC (rev 156)
@@ -10,24 +10,21 @@
## ****************************************************************************/
##@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, prevAns=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$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
}
storeData <- function()
@@ -37,115 +34,96 @@
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$acceptances[, start:end] <- z$accepts
+ z$candidates[start:end,, ] <- z$parameters
z$posteriorTot <- z$posteriorTot + colSums(z$parameters)
- for (i in npar)
+ for (group in 1:z$nGroup)
{
- z$posteriorMII <- z$posteriorMII +
- outer(z$parameters[i, ], z$parameters[i, ])
+ for (i in npar)
+ {
+ z$posteriorMII[group, , ] <- z$posteriorMII[group, ,] +
+ outer(z$parameters[i, group, ], z$parameters[i, group,])
+ }
}
- z$MHacceptances[start:end, ] <- z$MHaccepts
- z$MHrejections[start:end, ] <- z$MHrejects
- z$MHproportions[start:end, ] <-z$MHaccepts/ (z$MHaccepts + z$MHrejects)
+ z$MHacceptances[start:end, , ] <- z$MHaccepts
+ z$MHrejections[start:end, , ] <- z$MHrejects
+ z$MHproportions[start:end, , ] <- z$MHaccepts /
+ (z$MHaccepts + z$MHrejects)
z$sub <- z$sub + nrun
z
}
- ## we have removed the rejected values of betas and now want to fill in
- ## the gaps by duplicating the previous value. Need to use last one from
- ## previous or original theta if we have NA in first row.
- carryForward <- function(parameters, betas)
- {
- npar <- nrow(parameters)
- nbeta <- nrow(betas)
- if (npar < nbeta)
- {
- parameters <- rbind(betas[(nbeta - npar), ], parameters)
- }
- else
- {
- parameters <- rbind(z$theta[!z$basicRate], parameters)
- }
- parameters <-
- apply(parameters, 2, function(x)
- {
- x.pos <- which(!is.na(x))
- if (length(x.pos) == 0 || x.pos[1] != 1)
- {
- x.pos <- c(1, x.pos)
- }
- x[rep(x.pos, c(diff(x.pos),
- length(x) - x.pos[length(x.pos)] + 1))]
- }
- )
- parameters[-1, ]
- }
improveMH <- function(z, x, tiny=1.0e-15, desired=40, maxiter=100,
tolerance=15)
{
rescaleCGD <- function(iter)
{
- if (actual > desired)
- {
- u <- 2 - ((iter - actual) / (iter - desired))
- }
- else
- {
- u <- 1 / (2 - (actual / desired))
- }
- if (abs(actual - desired) <= tolerance)
- {
- number <<- number + 1
- if (number == 2) success <<- TRUE
- }
- else
- {
- number <<- 0
- }
+ u <- ifelse (actual > desired,
+ 2 - ((iter - actual) / (iter - desired)),
+ 1 / (2 - (actual / desired)))
+ number <<- ifelse(abs(actual - desired) <= tolerance,
+ number + 1, 0 )
+ success <<- number >= 2
u
}
iter <- 0
- number <- 0
- success <- FALSE
+ number <- rep(0, z$nGroup)
+ success <- rep(FALSE, z$nGroup)
repeat
{
iter <- iter + 1
- z <- MCMCcycle(z, nrunMH=1, nrunMHBatches=100)
+ z <- MCMCcycle(z, nrunMH=1, nrunMHBatches=100, change=FALSE)
actual <- z$BayesAcceptances ## acceptances
- ans <- rescaleCGD(100 * (z$observations - 1))
- z$scaleFactor <- z$scaleFactor * ans
- if (success | iter == maxiter)
+ ans <- rescaleCGD(100)
+ update <- number < 3
+ z$scaleFactors[update] <- z$scaleFactors[update] * ans[update]
+ cat(actual, ans, z$scaleFactors, '\n')
+ if (all(success) | iter == maxiter)
{
break
}
- if (z$scaleFactor < tiny)
+ if (any(z$scaleFactors < tiny))
{
cat('scalefactor < tiny\n')
browser()
}
}
cat('fine tuning took ', iter, ' iterations. Scalefactor:',
- z$scaleFactor, '\n')
- z
+ z$scaleFactors, '\n')
+ z
}
## ################################
## start of function proper
## ################################
- ## 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,
+ prevAns=prevAns)
z <- createStores(z)
z$sub <- 0
+ if (is.null(z$dfra) && is.null(dfra))
+ {
+ z <- getDFRA(z, 10)
+ }
+ else
+ {
+ if (!is.null(dfra))
+ {
+ z$dfra <- dfra
+ }
+ lambda <- z$theta[z$basicRate]
+ dfra <- z$dfra
+ z$dfra[z$basicRate, ] <- z$dfra[z$basicRate,] * lambda
+ z$dfra[, z$basicRate] <- z$dfra[, z$basicRate] * lambda
+ diag(z$dfra)[z$basicRate] <- lambda *
+ diag(dfra)[z$basicRate] + lambda * lambda *
+ colMeans(z$sf)[z$basicRate]
+ print(z$dfra)
+ z$dfra <- solve(z$dfra)
+ print(z$dfra)
+
+
+ }
z <- improveMH(z)
if (!save)
@@ -154,13 +132,16 @@
dev.new()
thetaplot = dev.cur()
dev.new()
- acceptsplot = dev.cur()
- }
+ ratesplot = dev.cur()
+ dev.new()
+ tseriesplot = dev.cur()
+ dev.new()
+ tseriesratesplot = dev.cur()
+ }
for (ii in 1:nwarm)
{
z <- MCMCcycle(z, nrunMH=4, nrunMHBatches=20)
- numm <- z$numm
}
for (ii in 1:nmain)
@@ -168,31 +149,38 @@
z <- MCMCcycle(z, nrunMH=nrunMH, nrunMHBatches=nrunMHBatches)
z <- storeData()
- numm <- z$numm
-
if (ii %% 10 == 0 && !save) ## do some plots
{
- cat('main after ii',ii,numm, '\n')
+ cat('main after ii', ii, '\n')
dev.set(thetaplot)
+ thetadf <-
+ lapply(1:z$nGroup, function(i)
+ {
+ data.frame(Group=rep(i, ii * nrunMHBatches),
+ z$candidates[1:(ii * nrunMHBatches), i, ])
+ }
+ )
+ thetadf <- do.call(rbind, thetadf)
basicRate <- z$basicRate
- lambdas <- z$lambdas
- lambdas[lambdas == 0] <- NA
- thetadf <- data.frame(lambdas, z$betas)
+ ##thetadf <- data.frame(z$candidates)
acceptsdf <- data.frame(z$MHproportions,
z$acceptances)
- lambdaNames <- paste(z$effects$name[basicRate],
- z$effects$shortName[basicRate],
- z$effects$period[basicRate],
- z$effects$group[basicRate], sep=".")
- betaNames <- paste(z$effects$name[!basicRate],
- z$effects$shortName[!basicRate], sep=".")
- names(thetadf) <- make.names(c(lambdaNames, betaNames),
- unique=TRUE)
+ ratesdf <- thetadf[, -1][, z$basicRate]
+ thetadf <- cbind(Group=thetadf[, 1], thetadf[, -1][, !z$basicRate])
+ 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=".")
+ names(ratesdf) <- rateNames
+ ratesdf <- cbind(Group=thetadf[, 1], ratesdf)
+ names(thetadf)[-1] <- make.names(thetaNames, unique=TRUE)
names(acceptsdf) <- c("InsDiag", "CancDiag", "Permute", "InsPerm",
"DelPerm", "InsMissing", "DelMissing",
"BayesAccepts")
- varnames <- paste(names(thetadf), sep="", collapse= " + ")
- varcall <- paste("~ ", varnames, sep="", collapse="")
+ varnames <- paste(names(thetadf)[-1], sep="", collapse= " + ")
+ varcall <- paste("~ ", varnames, " | Group", sep="", collapse="")
print(histogram(as.formula(varcall), data=thetadf, scales="free",
outer=TRUE, breaks=NULL, type="density",
panel=function(x, ...)
@@ -201,91 +189,154 @@
panel.densityplot(x, darg=list(na.rm=TRUE), ...)
}
))
- dev.set(acceptsplot)
- varnames <- paste(names(acceptsdf), sep="", collapse= " + ")
- varcall <- paste("~ ", varnames, sep="", collapse="")
- print(histogram(as.formula(varcall), data=acceptsdf,
- scales=list(x="same", y="free"),
+ dev.set(ratesplot)
+ varnames <- paste(names(ratesdf)[-1], sep="", collapse= " + ")
+ varcall <- paste("~ ", varnames, sep="", collapse="")
+ print(histogram(as.formula(varcall), data=ratesdf, scales="free",
outer=TRUE, breaks=NULL, type="density",
panel=function(x, ...)
{
panel.histogram(x, ...)
panel.densityplot(x, darg=list(na.rm=TRUE), ...)
- }))
+ }
+ ))
+ varnames <- paste(names(thetadf)[-1], sep="", collapse= " + ")
+ varcall <- paste(varnames, "~ 1:", ii * nrunMHBatches * z$nGroup,
+ " | Group", sep="", collapse="")
+ dev.set(tseriesplot)
+ print(xyplot(as.formula(varcall), data=thetadf, scales="free",
+ outer=TRUE))
+ varnames <- paste(names(ratesdf)[-1], sep="", collapse= " + ")
+ varcall <- paste(varnames, "~ 1:", ii * nrunMHBatches * z$nGroup,
+ sep="", collapse="")
+ dev.set(tseriesratesplot)
+ print(xyplot(as.formula(varcall), data=ratesdf, scales="free",
+ outer=TRUE))
+ ## dev.set(acceptsplot)
+ ## varnames <- paste(names(acceptsdf), sep="", collapse= " + ")
+ ## varcall <- paste("~ ", varnames, sep="", collapse="")
+ ## print(histogram(as.formula(varcall), data=acceptsdf,
+ ## scales=list(x="same", y="free"),
+ ## outer=TRUE, breaks=NULL, type="density",
+ ## panel=function(x, ...)
+ ## {
+ ## panel.histogram(x, ...)
+ ## panel.densityplot(x, darg=list(na.rm=TRUE), ...)
+ ## }))
}
}
+ z$FRAN <- NULL
z
}
-MCMCcycle <- function(z, nrunMH, nrunMHBatches)
+MCMCcycle <- function(z, nrunMH, nrunMHBatches, change=TRUE)
{
- f <- FRANstore() ## retrieve info
+ 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))
+ storeNrunMH <- z$nrunMH
+ z$nrunMH <- nrunMH
+ for (i in 1:nrunMHBatches)
+ {
+ ans <- z$FRAN(z, returnChains=TRUE, byGroup=TRUE)
+ z$chain <- ans$chain
+ z <- sampleParameters(z, change)
+ z$accepts[, i] <- z$accept
+ z$parameters[i, , ] <- z$thetaMat
+ z$MHaccepts[i, ,] <-
+ t(do.call(cbind,
+ tapply(ans$accepts, factor(z$callGrid[, 1]),
+ function(x)Reduce("+", x))))
+ z$MHrejects[i, , ] <-
+ t(do.call(cbind, tapply(ans$rejects, factor(z$callGrid[, 1]),
+ function(x)Reduce("+", x))))
+ z$MHaborts[i, , ] <- t(do.call(cbind,
+ tapply(ans$aborts,
+ factor(z$callGrid[, 1]),
+ function(x)Reduce("+", x))))
+ }
+ z$BayesAcceptances <- rowSums(z$accepts)
+ z$nrunMH <- storeNrunMH
+ z
+}
+sampleParameters <- function(z, change=TRUE)
+{
+ ## get a multivariate normal with covariance matrix dfra multiplied by a
+ ## scale factor which varies between groups
+ require(MASS)
+ thetaChanges <- t(sapply(1:z$nGroup, function(i)
+ {
+ tmp <- z$thetaMat[i, ]
+ use <- !is.na(z$thetaMat[i, ])
+ tmp[use] <-
+ mvrnorm(1, mu=rep(0, sum(use)),
+ Sigma=z$scaleFactors[i] *
+ z$dfra[use, use])
+ tmp
+ }
+ ))
- ## 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))))
+ thetaOld <- z$thetaMat
+ thetaOld[, z$basicRate] <- log(thetaOld[, z$basicRate])
+ thetaNew <- thetaOld + thetaChanges
- ## 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)
+ priorOld <- sapply(1:z$nGroup, function(i)
+ {
+ tmp <- thetaOld[i, ]
+ use <- !is.na(tmp)
+ dmvnorm(tmp[use], mean=rep(0, sum(use)),
+ sigma=z$priorSigma[use, use])
+ }
+ )
+ priorNew <- sapply(1:z$nGroup, function(i)
+ {
+ tmp <- thetaNew[i, ]
+ use <- !is.na(tmp)
+ dmvnorm(tmp[use], mean=rep(0, sum(use)),
+ sigma=z$priorSigma[use, use])
+ }
+ )
+ logpOld <- getProbs(z, z$chain)
+
+ storeNrunMH <- z$nrunMH
+ z$nrunMH <- 0
+ thetaNew[, z$basicRate] <- exp(thetaNew[, z$basicRate])
+ z$thetaMat <- thetaNew
+ resp <- z$FRAN(z, returnChains=TRUE, byGroup=TRUE)
+ logpNew <- getProbs(z, resp$chain)
+
+ proposalProbability <- priorNew - priorOld + logpNew - logpOld
+
+ z$accept <- log(runif(length(proposalProbability))) < proposalProbability
+ thetaOld[, z$basicRate] <- exp(thetaOld[, z$basicRate])
+ if (!change)
{
- ans <- .Call("MCMCcycle", PACKAGE=pkgname, f$pData, f$pModel,
- f$myeffects, as.integer(1), as.integer(1),
- z$scaleFactor, nrunMH, nrunMHBatches)
+ z$thetaMat <- thetaOld
}
else
{
- if (z$int2 == 1)
- {
- anss <- apply(callGrid, 1, doMCMCcycle, z$scaleFactor,
- nrunMH, nrunMHBatches)
- }
- 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]]
+ ##print(z$thetaMat)
+ z$thetaMat[!z$accept, ] <- thetaOld[!z$accept, ]
+ }
+ z$nrunMH <- storeNrunMH
+ ## cat(thetaNew, priorNew, priorOld, logpNew, logpOld,
+ ## exp(proposalProbability),
+ ## z$accept,'\n')
z
}
-doMCMCcycle <- function(x, scaleFactor, nrunMH, nrunMHBatches)
-{
- group <- x[1]
- period <- x[2]
- f <- FRANstore()
- .Call("MCMCcycle", PACKAGE=pkgname, f$pData, f$pModel,
- f$myeffects, as.integer(period),
- as.integer(group),
- scaleFactor, nrunMH, nrunMHBatches)
-}
-initializeBayes <- function(data, effects, model, nbrNodes, seed)
+initializeBayes <- function(data, effects, model, nbrNodes, seed, priorSigma,
+ prevAns)
{
## initialise
set.seed(seed)
Report(openfiles=TRUE, type="n") #initialise with no file
z <- NULL
+ z$Phase <- 1
+ z$Deriv <- FALSE
z$FinDiff.method <- FALSE
z$maxlike <- TRUE
model$maxlike <- TRUE
@@ -298,10 +349,10 @@
{
set.seed(model$randomSeed)
}
- model$FRAN <- getFromNamespace(model$FRANname, pos=grep("RSiena",
- search())[1])
- z <- model$FRAN(z, model, INIT=TRUE, data=data, effects=effects)
-
+ z$FRAN <- getFromNamespace(model$FRANname, pos=grep("RSiena",
+ search())[1])
+ z <- z$FRAN(z, model, INIT=TRUE, data=data, effects=effects,
+ prevAns=prevAns)
is.batch(TRUE)
WriteOutTheta(z)
@@ -317,14 +368,133 @@
clusterCall(z$cl, storeinFRANstore, FRANstore())
clusterCall(z$cl, FRANstore)
clusterCall(z$cl, initializeFRAN, z, model,
- initC = TRUE, profileData=FALSE, returnDeps=FALSE)
+ initC = TRUE, profileData=FALSE, returnDeps=FALSE)
clusterSetupRNG(z$cl,
seed = as.integer(runif(6, max=.Machine$integer.max)))
}
- z$numm <- 20
- z$scaleFactor <- 1
+ z$scaleFactors <- rep(1, z$nGroup)
+ ## z$returnDataFrame <- TRUE # chains come back as data frames not lists
+ z$returnChains <- TRUE
+ if (is.null(priorSigma))
+ {
+ z$priorSigma <- diag(z$pp) * 10000
+ }
+ else
+ {
+ z$priorSigma <- priorSigma
+ }
+ z$ratePositions <- lapply(z$rateParameterPosition, unlist)
+ for (i in 1:z$nGroup)
+ {
+ use <- rep(FALSE, z$pp)
+ use[z$ratePositions[[i]]] <- TRUE
+ use[!z$basicRate] <- TRUE
+ z$thetaMat[i, !use] <- NA
+ }
+ z$callGrid <- cbind(rep(1:z$nGroup, z$groupPeriods),
+ as.vector(unlist(sapply(z$groupPeriods,
+ function(x) 1:x))))
+ z
+}
-
+getDFRA <- 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$ssc <- matrix(0, nrow=n, ncol=z$pp)
+ z$Deriv <- TRUE
+ for (i in 1:n)
+ {
+ ans <- z$FRAN(z)
+ z$sdf[i, , ] <- ans$dff
+ z$ssc[i, ] <- colSums(ans$fra)
+ }
+ dfra <- t(apply(z$sdf, c(2, 3), mean))
+ ssc <- colMeans(z$ssc)
+ z$dfra <- dfra
+ lambda <- z$theta[z$basicRate]
+ z$dfra[z$basicRate, ] <- z$dfra[z$basicRate,] * lambda
+ z$dfra[, z$basicRate] <- z$dfra[, z$basicRate] * lambda
+ diag(z$dfra)[z$basicRate] <- lambda *
+ diag(dfra)[z$basicRate] + lambda * lambda * ssc[z$basicRate]
+ ##print(z$dfra)
+ z$dfra <- solve(z$dfra)
+ ##print(z$dfra)
+ ##$eS <- eigen(z$dfra, symmetric=TRUE, EISPACK=TRUE)
+ ##z$ev <- z$eS$values
+ ##z$covFactor <- z$eS$vectors %*% diag(sqrt(pmax(z$ev, 0)), z$pp)
z
}
+
+getLikelihood <- function(chain, nactors, lambda, simpleRates)
+{
+ loglik <- 0
+ ncvals <- sapply(chain, function(x)x[[3]])
+ nc <- nactors
+ nc[] <- 0
+ ncvals <- table(ncvals)
+ nc[names(ncvals)] <- ncvals
+ logChoiceProb <- sapply(chain, function(x)x[[9]])
+ logOptionSetProb <- sapply(chain, function(x)x[[8]])
+ loglik <- sum(logChoiceProb) + sum(logOptionSetProb)
+ #print(sum(logOptionSetProb))
+ if (simpleRates)
+ {
+ loglik <- loglik - sum(nactors * lambda) + sum(nc * log(lambda))# -
+ # sum(lfactorial(nc))
+ }
+ else
+ {
+ mu <- attr(chain, "mu")
+ sigma <- sqrt(attr(chain, "sigma2"))
+ finalReciprocalRate <- attr(chain, "finalReciprocalRate")
+ loglik <- loglik + dnorm(1, mu, sigma, log=TRUE) +
+ log(finalReciprocalRate)
+ }
+ loglik
+}
+
+getProbs <- function(z, chain)
+{
+ sapply(1: length(chain), function(i)
+ {
+ groupChain <- chain[[i]]
+ sum(sapply(1:length(groupChain), function(j)
+ {
+ periodChain <- groupChain[[j]]
+ theta1 <- z$thetaMat[i, ]
+ k <- z$rateParameterPosition[[i]][[j]]
+ getLikelihood(periodChain, z$nactors[[i]], theta1[k],
+ z$simpleRates)
+ }
+ ))
+ }
+ )
+}
+
+flattenChains <- function(zz)
+{
+ for (i in 1:length(zz)) ##group
+ {
+ for (j in 1:length(zz[[i]])) ## period
+ {
+ attr(zz[[i]][[j]], "group") <- i
+ attr(zz[[i]][[j]], "period") <- j
+ }
+ }
+ zz <- do.call(c, zz)
+ zz
+}
+
+dmvnorm <- function(x, mean , sigma)
+{
+ if (is.vector(x))
+ {
+ x <- matrix(x, ncol=length(x))
+ }
+ distval <- mahalanobis(x, center=mean, cov=sigma)
+ logdet <- sum(log(eigen(sigma, symmetric=TRUE, only.values=TRUE)$values))
+ -(ncol(x) * log(2 * pi) + logdet + distval) / 2
+}
Modified: pkg/RSiena/R/effects.r
===================================================================
--- pkg/RSiena/R/effects.r 2011-06-12 16:27:55 UTC (rev 155)
+++ pkg/RSiena/R/effects.r 2011-06-13 14:32:32 UTC (rev 156)
@@ -696,12 +696,12 @@
objEffects <-
objEffects[!objEffects$shortName == "density", ]
}
- if (attr(xx$depvars[[i]],'uponly') || attr(xx$depvars[[i]],
- 'downonly'))
- {
- 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
Modified: pkg/RSiena/R/initializeFRAN.r
===================================================================
--- pkg/RSiena/R/initializeFRAN.r 2011-06-12 16:27:55 UTC (rev 155)
+++ pkg/RSiena/R/initializeFRAN.r 2011-06-13 14:32:32 UTC (rev 156)
@@ -455,6 +455,14 @@
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)
+ ## create a grid of periods with group names in case want to
+ ## parallelize using this
+ groupPeriods <- attr(f, "groupPeriods")
+ z$callGrid <- cbind(rep(1:nGroup, groupPeriods - 1),
+ as.vector(unlist(sapply(groupPeriods - 1,
+ function(x) 1:x))))
}
}
@@ -488,7 +496,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 +582,16 @@
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/RSiena/R/maxlikec.r
===================================================================
--- pkg/RSiena/R/maxlikec.r 2011-06-12 16:27:55 UTC (rev 155)
+++ pkg/RSiena/R/maxlikec.r 2011-06-13 14:32:32 UTC (rev 156)
@@ -12,12 +12,15 @@
##@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,
+ returnDataFrame=FALSE)
{
if (INIT || initC) ## initC is to initialise multiple C processes in phase3
{
z <- initializeFRAN(z, x, data, effects, prevAns, initC,
profileData=profileData, returnDeps=returnDeps)
+ z$returnDataFrame <- returnDataFrame
+ z$returnChains <- returnChains
if (initC)
{
return(NULL)
@@ -45,22 +48,24 @@
{
returnDeps <- z$returnDeps
}
- ## create a grid of periods with group names in case want to parallelize
- ## using this
- groupPeriods <- attr(f, "groupPeriods")
- callGrid <- cbind(rep(1:f$nGroup, groupPeriods - 1),
- as.vector(unlist(sapply(groupPeriods - 1,
- function(x) 1:x))))
+ callGrid <- z$callGrid
## 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)
{
+ if (byGroup)
+ {
+ theta <- z$thetaMat[1,]
+ }
+ else
+ {
+ theta <- z$theta
+ }
ans <- .Call('mlPeriod', PACKAGE=pkgname, z$Deriv, f$pData,
- f$pModel, f$myeffects, z$theta,
+ f$pModel, f$myeffects, theta,
returnDeps, 1, 1, z$nrunMH, z$addChainToStore,
z$needChangeContributions, z$returnDataFrame,
- returnChains)
+ z$returnChains)
ans[[6]] <- list(ans[[6]])
ans[[7]] <- list(ans[[7]])
}
@@ -68,18 +73,19 @@
{
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)
+ z$returnChains, byGroup, z$theta)
}
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)
+ z$returnDataFrame, z$returnChains, byGroup, z$theta)
}
## reorganize the anss so it looks like the normal one
ans <- NULL
@@ -108,30 +114,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
+ 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 +153,101 @@
## {
## periodNos <- periodNo:(periodNo + length(sims[[i]][[j]]) - 1)
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/rsiena -r 156
More information about the Rsiena-commits
mailing list