[Rsiena-commits] r337 - in pkg: RSiena RSiena/R RSiena/man RSienaTest RSienaTest/R RSienaTest/doc RSienaTest/man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Dec 5 22:44:36 CET 2018
Author: tomsnijders
Date: 2018-12-05 22:44:35 +0100 (Wed, 05 Dec 2018)
New Revision: 337
Added:
pkg/RSienaTest/doc/RateEffects.txt
pkg/RSienaTest/doc/notes_sienaBayes.txt
Modified:
pkg/RSiena/ChangeLog
pkg/RSiena/DESCRIPTION
pkg/RSiena/R/effects.r
pkg/RSiena/R/effectsMethods.r
pkg/RSiena/R/initializeFRAN.r
pkg/RSiena/R/iwlsm.R
pkg/RSiena/R/sienaprint.r
pkg/RSiena/man/RSiena-package.Rd
pkg/RSiena/man/effectsDocumentation.Rd
pkg/RSiena/man/includeEffects.Rd
pkg/RSiena/man/setEffect.Rd
pkg/RSiena/man/sienaAlgorithmCreate.Rd
pkg/RSienaTest/ChangeLog
pkg/RSienaTest/DESCRIPTION
pkg/RSienaTest/R/bayesTest.r
pkg/RSienaTest/R/effectsMethods.r
pkg/RSienaTest/R/initializeFRAN.r
pkg/RSienaTest/R/sienaBayes.r
pkg/RSienaTest/R/sienaprint.r
pkg/RSienaTest/doc/RSiena.bib
pkg/RSienaTest/doc/RSiena_Manual.pdf
pkg/RSienaTest/doc/RSiena_Manual.tex
pkg/RSienaTest/man/RSiena-package.Rd
pkg/RSienaTest/man/extract.sienaBayes.Rd
pkg/RSienaTest/man/getEffects.Rd
pkg/RSienaTest/man/includeEffects.Rd
pkg/RSienaTest/man/setEffect.Rd
pkg/RSienaTest/man/siena07.Rd
pkg/RSienaTest/man/sienaAlgorithmCreate.Rd
pkg/RSienaTest/man/sienaBayes.Rd
Log:
Version 1.2-14: new sienaBayes version, some extra docs
Modified: pkg/RSiena/ChangeLog
===================================================================
--- pkg/RSiena/ChangeLog 2018-10-30 11:22:36 UTC (rev 336)
+++ pkg/RSiena/ChangeLog 2018-12-05 21:44:35 UTC (rev 337)
@@ -1,3 +1,19 @@
+2018-12-05 R-Forge Revision 337, packages version 1.2-14.
+Changes in RSiena and RSienaTest:
+ * Check for length of parameter mult for ML estimation (initializeFRAN),
+ extend explanation of mult in help page for sienaAlgorithmCreate.
+Changes in RSienaTest:
+ * extract.sienaBayes corrected.
+ * sienaBayes: subphases for initial global parameter estimation decreased
+ to 2;
+ prior variances for fixed parameters introduced;
+ after initialization by MoM estimation of multi-group data set,
+ a precision weighted mean (called "Kelley") of this estimate and the
+ prior mean is used for initial parameter values;
+ prewarming phase introduced before improveMH;
+ in case prevBayes is used, also parameter nImproveMH in the function
+ call of sienaBayes supersedes this parameter in the prevBayes object.
+
2018-10-29 R-Forge Revision 336, packages version 1.2-13.
Changes in RSiena and RSienaTest:
* Correct error in print01Report that occurred for changing dyadic covariates
@@ -2,3 +18,3 @@
given as lists of sparse matrices.
- * Correct error in sienaGroupCreate for centered actor covariates.
+ * Correct error in sienaGroupCreate for non-centered actor covariates.
* Drop dependence on utils (the packages did not depend on it).
@@ -34,6 +50,12 @@
reported timing changed to elapsed system time;
some reordering of parameters.
This also required some corresponding changes in sienaprint.r.
+ * DESCRIPTION: Added 'Encoding: UTF-8' because sienaGOF.r complained about
+ non-ASCII.
+ * Data.cpp: Add a new 'SimNetworkData' type: This is for variables that
+ behave like dependent variable for effects but are not changed by the
+ simulation. Currently used for the primary setting network.
+ * PrimaryLayer: A network layer holding the full primary setting as network.
2018-05-13 CRAN for RSiena, version 1.2-12.
Modified: pkg/RSiena/DESCRIPTION
===================================================================
--- pkg/RSiena/DESCRIPTION 2018-10-30 11:22:36 UTC (rev 336)
+++ pkg/RSiena/DESCRIPTION 2018-12-05 21:44:35 UTC (rev 337)
@@ -2,8 +2,8 @@
Package: RSiena
Type: Package
Title: Siena - Simulation Investigation for Empirical Network Analysis
-Version: 1.2-13
-Date: 2018-10-29
+Version: 1.2-14
+Date: 2018-12-05
Author: Ruth Ripley, Krists Boitmanis, Tom A.B. Snijders, Felix Schoenenberger
Depends: R (>= 2.15.0)
Imports: Matrix, tcltk, lattice, parallel, MASS, methods
Modified: pkg/RSiena/R/effects.r
===================================================================
--- pkg/RSiena/R/effects.r 2018-10-30 11:22:36 UTC (rev 336)
+++ pkg/RSiena/R/effects.r 2018-12-05 21:44:35 UTC (rev 337)
@@ -373,8 +373,8 @@
objEffects$type == "eval",
c('include', "initialValue", "untrimmedValue")] <-
list(TRUE, starts$degree, starts$untrimmed)
- objEffects[objEffects$shortName=='transTriads' &
- objEffects$type=='eval','include'] <- TRUE
+# objEffects[objEffects$shortName=='transTriads' &
+# objEffects$type=='eval','include'] <- TRUE
}
else
{
Modified: pkg/RSiena/R/effectsMethods.r
===================================================================
--- pkg/RSiena/R/effectsMethods.r 2018-10-30 11:22:36 UTC (rev 336)
+++ pkg/RSiena/R/effectsMethods.r 2018-12-05 21:44:35 UTC (rev 337)
@@ -129,10 +129,13 @@
}
if (includeRandoms)
{
+ nfeff <- sum((!x$randomEffects) & x$include & (!x$fix) & (!x$basicRate))
nreff <- sum(x$randomEffects & x$include & (!x$fix))
nrate <- sum(x$basicRate & x$include & (x$group==2) & (!x$fix))
if (sum(x$group != 1) > 0) # else there is only one group, and counting should be different.
{
+ cat('Length of priorSigEta for sienaBayes, if used, should be',
+ nfeff, '.\n')
cat('Dimensions of priorMu and priorSigma for sienaBayes should be',
nreff, '+', nrate, '=', nreff+nrate,'.\n')
}
Modified: pkg/RSiena/R/initializeFRAN.r
===================================================================
--- pkg/RSiena/R/initializeFRAN.r 2018-10-30 11:22:36 UTC (rev 336)
+++ pkg/RSiena/R/initializeFRAN.r 2018-12-05 21:44:35 UTC (rev 337)
@@ -703,6 +703,13 @@
z$maxlikeTargets <- rowSums(ans)
z$maxlikeTargets2 <- ans
z$mult <- x$mult
+ length.nrunMH <- length(colSums(z$maxlikeTargets2[z$requestedEffects$basicRate,, drop=FALSE ]))
+ if ((length(z$mult) >= 2) &&
+ (length(z$mult) != length.nrunMH))
+ {
+ stop(paste('Length of parameter mult in the algorithm object is incorrect',
+ '(should be 1 or', length.nrunMH, ').'))
+ }
z$nrunMH <- floor(z$mult *
colSums(z$maxlikeTargets2[z$requestedEffects$basicRate,
, drop=FALSE ]))
Modified: pkg/RSiena/R/iwlsm.R
===================================================================
--- pkg/RSiena/R/iwlsm.R 2018-10-30 11:22:36 UTC (rev 336)
+++ pkg/RSiena/R/iwlsm.R 2018-12-05 21:44:35 UTC (rev 337)
@@ -17,7 +17,7 @@
#/******************************************************************************
# * SIENA: Simulation Investigation for Empirical Network Analysis
# *
-# * Web: http://www.stats.ox.ac.uk/~snidjers/siena
+# * Web: http://www.stats.ox.ac.uk/~snijders/siena
# *
# * File: iwlsm.r
# *
Modified: pkg/RSiena/R/sienaprint.r
===================================================================
--- pkg/RSiena/R/sienaprint.r 2018-10-30 11:22:36 UTC (rev 336)
+++ pkg/RSiena/R/sienaprint.r 2018-12-05 21:44:35 UTC (rev 337)
@@ -1427,12 +1427,14 @@
}
else
{
+ ncr <- max(sapply(x$effectName, nchar)) + 1
+ codestring <- paste("%-", ncr, "s", sep="") # the minus signs leads to right padding
cat("Bayesian estimation.\n")
cat("Prior distribution:\n")
cat("\nMu ")
for (i in seq(along=x$priorMu))
{
- cat(x$effectName[x$varyingParametersInGroup][i],
+ cat(sprintf(codestring, x$effectName[x$varyingParametersInGroup][i]),
sprintf("%8.4f", x$priorMu[i,1]),"\n ")
}
cat("\nSigma ")
@@ -1445,6 +1447,32 @@
{
cat("\nKappa ",sprintf("%8.4f", x$priorKappa),"\n")
}
+ if (!is.null(x$anyPriorEta))
+ {
+ cat("\nEta ")
+ if (x$anyPriorEta)
+ {
+ cat("\nFor the fixed parameters, prior variance:\n")
+ var.eta <- rep(NA, length(x$set2prior))
+ var.eta[x$set2prior] <- 1/(2*x$priorPrecEta)
+ for (i in seq(along=x$set2prior))
+ {
+ cat(sprintf(codestring, x$effectName[!x$varyingParametersInGroup][i]))
+ if (x$set2prior[i])
+ {
+ cat(sprintf("%8.4f", var.eta[i]),"\n ")
+ }
+ else
+ {
+ cat(" (constant prior) \n")
+ }
+ }
+ }
+ else
+ {
+ cat("\nFor the fixed parameters, constant prior.\n")
+ }
+ }
# cat("\nFor the basic rate parameters, ")
# cat("the prior is on the square root scale.\n\n")
}
@@ -1452,14 +1480,19 @@
{
cat("\nBasic rates parameters are treated as incidental parameters.\n\n")
}
- cat("\nAlgorithm specifications were nwarm =",x$nwarm,", nmain =", x$nmain,
+ cat("\nAlgorithm specifications were ")
+ if (!is.null(x$nprewarm))
+ {
+ cat(" nprewarm =",x$nprewarm,",")
+ }
+ cat(" nwarm =",x$nwarm,", nmain =", x$nmain,
", nrunMHBatches =", x$nrunMHBatches,
", nImproveMH =", x$nImproveMH,
",\n nSampVarying =", x$nSampVarying, ", nSampConst =", x$nSampConst,
", mult =", x$mult, ".\n")
if (!is.null(nfirst))
{
- cat("For the results, nwarm is superseded by nfirst = ", nfirst, ".")
+ cat("For these results, nwarm is superseded by nfirst = ", nfirst, ".")
}
if (ntot < x$nwarm + x$nmain)
{
@@ -1486,7 +1519,7 @@
cat(sprintf("%4.2f",
colMeans(x$ThinBayesAcceptances[first:ntot,])/x$nrunMHBatches),
fill=TRUE,"\n")
- cat("This should ideally be close to 0.25.\n")
+ cat("This should ideally be between 0.15 and 0.50.\n")
}
print.sienaBayesFit(x, nfirst)
if (ntot > first+2)
Modified: pkg/RSiena/man/RSiena-package.Rd
===================================================================
--- pkg/RSiena/man/RSiena-package.Rd 2018-10-30 11:22:36 UTC (rev 336)
+++ pkg/RSiena/man/RSiena-package.Rd 2018-12-05 21:44:35 UTC (rev 337)
@@ -44,8 +44,8 @@
\tabular{ll}{
Package: \tab RSiena\cr
Type: \tab Package\cr
- Version: \tab 1.2-13\cr
- Date: \tab 2018-10-29\cr
+ Version: \tab 1.2-14\cr
+ Date: \tab 2018-12-05\cr
Depends: \tab R (>= 3.0.0)\cr
Imports: \tab Matrix\cr
Suggests: \tab tcltk, network, codetools, lattice, MASS, parallel,
Modified: pkg/RSiena/man/effectsDocumentation.Rd
===================================================================
--- pkg/RSiena/man/effectsDocumentation.Rd 2018-10-30 11:22:36 UTC (rev 336)
+++ pkg/RSiena/man/effectsDocumentation.Rd 2018-12-05 21:44:35 UTC (rev 337)
@@ -32,7 +32,7 @@
a table, either latex or html. This table presents all the available
effects present in this version of RSiena, not delimited
by a particular data set. The default file name is "effects.tex" or
- "effects.pdf", respectively. The latter file is also shown in the browser
+ "effects.html", respectively. The latter file is also shown in the browser
when requesting\cr
\code{RShowDoc("effects", package="RSiena")}
Modified: pkg/RSiena/man/includeEffects.Rd
===================================================================
--- pkg/RSiena/man/includeEffects.Rd 2018-10-30 11:22:36 UTC (rev 336)
+++ pkg/RSiena/man/includeEffects.Rd 2018-12-05 21:44:35 UTC (rev 337)
@@ -15,8 +15,7 @@
\item{myeff}{a Siena effects object as created by \code{\link{getEffects}}
}
\item{\dots}{
-short names to identify the effects which should be included or excluded.
-}
+short names to identify the effects which should be included or excluded.}
\item{include}{Boolean. default TRUE, but can be switched to FALSE to
turn off an effect.}
\item{name}{Name of network for which effects are being
@@ -44,19 +43,22 @@
to interactions between effects, but to dependence of effects on
other variables in the data set.
The arguments should identify the effects completely.
+
+ The function \code{includeEffects} operates as an interface
+ setting the "include" column on selected rows of the effects object,
+ to the value requested (TRUE or FALSE).
The value of \code{myeff$initialValue} can be set by function
\code{\link{setEffect}}.
+ The function \code{\link{setEffect}} can operate on the effects object
+ in a more detailed way, e.g., setting the value of \code{myeff$initialValue},
+ but applies to one effect at the time.
A list of all effects available in a given effects object (e.g.,
\code{myeff}), including their names of dependent variables,
effect names, short names,
and values of interaction1 and interaction2 (if any),
is obtained by executing \code{\link{effectsDocumentation}(myeff)}.
-
- The function \code{}includeEffects operates as an interface
- setting the "include" column on selected rows of the effects object,
- to the value requested (TRUE or FALSE).
}
\value{
An updated version of the input effects object, with the include, test, and fix
Modified: pkg/RSiena/man/setEffect.Rd
===================================================================
--- pkg/RSiena/man/setEffect.Rd 2018-10-30 11:22:36 UTC (rev 336)
+++ pkg/RSiena/man/setEffect.Rd 2018-12-05 21:44:35 UTC (rev 337)
@@ -74,17 +74,18 @@
}
\item{group}{
Number of group if basic rate.
+ Only relevant for \code{\link{sienaGroup}} data sets.
}
-\item{character}{Boolean: is the short name a character string or not.}
+\item{character}{Boolean: whether the short name is a character string.}
}
\details{
The arguments \code{shortName}, \code{name}, \code{type},
\code{interaction1}, \code{interaction2}, \code{effect1}, \code{effect2},
-and \code{effect3} should identify the effects completely.
-The column elements of the output effects object for \code{parm},
-\code{fix}, \code{test}, \code{randomEffects},
-\code{initialValue}, and \code{timeDummy}
-will be set to the values requested.}
+\code{effect3}, and \code{group} should identify the effects completely.
+
+The column elements of the output effects object for \code{parameter},
+\code{fix}, \code{test}, \code{randomEffects}, \code{initialValue},
+and \code{timeDummy} will be set to the values requested.}
\value{
An object of class \code{\link{sienaEffects}} or
\code{\link{sienaGroupEffects}}. This will be an updated version of the
Modified: pkg/RSiena/man/sienaAlgorithmCreate.Rd
===================================================================
--- pkg/RSiena/man/sienaAlgorithmCreate.Rd 2018-10-30 11:22:36 UTC (rev 336)
+++ pkg/RSiena/man/sienaAlgorithmCreate.Rd 2018-12-05 21:44:35 UTC (rev 337)
@@ -144,7 +144,11 @@
\item{mult}{Multiplication factor for maximum likelihood and Bayes. Number of
steps per iteration is set to this multiple of the total distance
between the observations at start and finish of the wave.
- Decreasing \code{mult} below a certain value has no further effect.}
+ Decreasing \code{mult} below a certain value has no further effect.\cr
+ This can be either a number (which needs to be positive) or a vector
+ of numbers, of length equal to the number of basic rate parameters in the
+ model, i.e., the number of periods times the number of dependent
+ variables.}
\item{simOnly}{Logical: If TRUE, then the calculation of the covariance
matrix and standard errors of the estimates at the end of
Phase 3 of the estimation algorithm in function siena07 is skipped.
Modified: pkg/RSienaTest/ChangeLog
===================================================================
--- pkg/RSienaTest/ChangeLog 2018-10-30 11:22:36 UTC (rev 336)
+++ pkg/RSienaTest/ChangeLog 2018-12-05 21:44:35 UTC (rev 337)
@@ -1,3 +1,19 @@
+2018-12-05 R-Forge Revision 337, packages version 1.2-14.
+Changes in RSiena and RSienaTest:
+ * Check for length of parameter mult for ML estimation (initializeFRAN),
+ extend explanation of mult in help page for sienaAlgorithmCreate.
+Changes in RSienaTest:
+ * extract.sienaBayes corrected.
+ * sienaBayes: subphases for initial global parameter estimation decreased
+ to 2;
+ prior variances for fixed parameters introduced;
+ after initialization by MoM estimation of multi-group data set,
+ a precision weighted mean (called "Kelley") of this estimate and the
+ prior mean is used for initial parameter values;
+ prewarming phase introduced before improveMH;
+ in case prevBayes is used, also parameter nImproveMH in the function
+ call of sienaBayes supersedes this parameter in the prevBayes object.
+
2018-10-29 R-Forge Revision 336, packages version 1.2-13.
Changes in RSiena and RSienaTest:
* Correct error in print01Report that occurred for changing dyadic covariates
@@ -2,3 +18,3 @@
given as lists of sparse matrices.
- * Correct error in sienaGroupCreate for centered actor covariates.
+ * Correct error in sienaGroupCreate for non-centered actor covariates.
* Drop dependence on utils (the packages did not depend on it).
Modified: pkg/RSienaTest/DESCRIPTION
===================================================================
--- pkg/RSienaTest/DESCRIPTION 2018-10-30 11:22:36 UTC (rev 336)
+++ pkg/RSienaTest/DESCRIPTION 2018-12-05 21:44:35 UTC (rev 337)
@@ -2,8 +2,8 @@
Package: RSienaTest
Type: Package
Title: Siena - Simulation Investigation for Empirical Network Analysis
-Version: 1.2-13
-Date: 2018-10-29
+Version: 1.2-14
+Date: 2018-12-05
Author: Ruth Ripley, Krists Boitmanis, Tom A.B. Snijders, Felix Schoenenberger
Depends: R (>= 2.15.0)
Imports: Matrix, tcltk, lattice, parallel, MASS, RUnit, methods
Modified: pkg/RSienaTest/R/bayesTest.r
===================================================================
--- pkg/RSienaTest/R/bayesTest.r 2018-10-30 11:22:36 UTC (rev 336)
+++ pkg/RSienaTest/R/bayesTest.r 2018-12-05 21:44:35 UTC (rev 337)
@@ -283,8 +283,14 @@
stop('all elements of zlist must be sienaBayesFit objects')
}
- #browser()
niter <- sapply(zlist,function(z){dim(z$ThinPosteriorMu)[1]})
+ if (length(niter) > 1)
+ {
+ if (var(niter) > 1e-6)
+ {
+ stop("Length of these sienaBayes objects is not constant")
+ }
+ }
nit <- niter[1] - nfirst + 1
if (extracted %in% c("all", "rates")){
@@ -307,7 +313,8 @@
if (extracted %in% c("all", "objective", "varying")){
# mu
- nind <- sum(zlist[[1]]$varyingParametersInGroup)
+ varyings <- zlist[[1]]$varyingInEstimated
+ nind <- sum(varyings)
res2 <- array(dim = c(nit, length(zlist), (2*nind)))
if (nind <= 0)
{
@@ -316,10 +323,10 @@
}
else
{
- EffName <- getNames(zlist[[1]])[zlist[[1]]$varyingParametersInGroup]
+ EffName <- getNames(zlist[[1]])[varyings]
for (h in 1:length(zlist)){
res2[,h,1:nind] <-
- zlist[[h]]$ThinPosteriorMu[(niter[h]-nit+1):niter[h], ,drop=FALSE]
+ zlist[[h]]$ThinPosteriorMu[(niter[h]-nit+1):niter[h], , drop=FALSE]
postsig <- zlist[[h]]$ThinPosteriorSigma
postsigg <- matrix(NA,dim(postsig)[1],nind)
if (sdLog){
Modified: pkg/RSienaTest/R/effectsMethods.r
===================================================================
--- pkg/RSienaTest/R/effectsMethods.r 2018-10-30 11:22:36 UTC (rev 336)
+++ pkg/RSienaTest/R/effectsMethods.r 2018-12-05 21:44:35 UTC (rev 337)
@@ -129,10 +129,13 @@
}
if (includeRandoms)
{
+ nfeff <- sum((!x$randomEffects) & x$include & (!x$fix) & (!x$basicRate))
nreff <- sum(x$randomEffects & x$include & (!x$fix))
nrate <- sum(x$basicRate & x$include & (x$group==2) & (!x$fix))
if (sum(x$group != 1) > 0) # else there is only one group, and counting should be different.
{
+ cat('Length of priorSigEta for sienaBayes, if used, should be',
+ nfeff, '.\n')
cat('Dimensions of priorMu and priorSigma for sienaBayes should be',
nreff, '+', nrate, '=', nreff+nrate,'.\n')
}
Modified: pkg/RSienaTest/R/initializeFRAN.r
===================================================================
--- pkg/RSienaTest/R/initializeFRAN.r 2018-10-30 11:22:36 UTC (rev 336)
+++ pkg/RSienaTest/R/initializeFRAN.r 2018-12-05 21:44:35 UTC (rev 337)
@@ -703,6 +703,13 @@
z$maxlikeTargets <- rowSums(ans)
z$maxlikeTargets2 <- ans
z$mult <- x$mult
+ length.nrunMH <- length(colSums(z$maxlikeTargets2[z$requestedEffects$basicRate,, drop=FALSE ]))
+ if ((length(z$mult) >= 2) &&
+ (length(z$mult) != length.nrunMH))
+ {
+ stop(paste('Length of parameter mult in the algorithm object is incorrect',
+ '(should be 1 or', length.nrunMH, ').'))
+ }
z$nrunMH <- floor(z$mult *
colSums(z$maxlikeTargets2[z$requestedEffects$basicRate,
, drop=FALSE ]))
Modified: pkg/RSienaTest/R/sienaBayes.r
===================================================================
--- pkg/RSienaTest/R/sienaBayes.r 2018-10-30 11:22:36 UTC (rev 336)
+++ pkg/RSienaTest/R/sienaBayes.r 2018-12-05 21:44:35 UTC (rev 337)
@@ -18,12 +18,13 @@
sienaBayes <- function(data, effects, algo, saveFreq=100,
initgainGlobal=0.1, initgainGroupwise = 0.02, initfgain=0.2,
gamma=0.05, initML=FALSE,
+ priorSigEta=NULL,
priorMu=NULL, priorSigma=NULL, priorDf=NULL, priorKappa=NULL,
priorRatesFromData=2,
frequentist=FALSE, incidentalBasicRates=FALSE,
reductionFactor=0.5,
delta=1e-4,
- nwarm=50, nmain=250, nrunMHBatches=20,
+ nprewarm=50, nwarm=50, nmain=250, nrunMHBatches=20,
nSampVarying=1, nSampConst=1, nSampRates=0,
nImproveMH=100, targetMHProb=0.25,
lengthPhase1=round(nmain/5), lengthPhase3=round(nmain/5),
@@ -33,13 +34,13 @@
nbrNodes=1, clusterType=c("PSOCK", "FORK"),
getDocumentation=FALSE)
{
- ##@createStores internal sienaBayes Bayesian set up stores
- createStores <- function()
- {
- # npar <- length(z$theta)
+ ##@createStores internal sienaBayes Bayesian set up stores
+ createStores <- function()
+ {
+ # npar <- length(z$theta)
numberRows <- nmain * nrunMHBatches
- # z$candidates <<- array(NA, dim=c(numberRows, z$nGroup, npar))
- # z$acceptances <<- matrix(NA, nrow=numberRows, ncol=z$nGroup)
+ # z$candidates <<- array(NA, dim=c(numberRows, z$nGroup, npar))
+ # z$acceptances <<- matrix(NA, nrow=numberRows, ncol=z$nGroup)
if (storeAll)
{
z$acceptances <<- matrix(NA, nrow=numberRows, ncol=z$nGroup)
@@ -77,12 +78,12 @@
# z$ThinP3EtaScores <<- matrix(NA, z$lengthPhase3, z$p2)
} # end createStores
- ##@storeData internal sienaBayes; put data in stores
- storeData <- function()
- {
- start <- z$sub + 1
- nrun <- ncol(zm$accepts)
- end <- start + nrun - 1
+ ##@storeData internal sienaBayes; put data in stores
+ storeData <- function()
+ {
+ start <- z$sub + 1
+ nrun <- ncol(zm$accepts)
+ end <- start + nrun - 1
if (storeAll)
{
z$acceptances[start:end, ] <<- zm$accepts
@@ -92,36 +93,36 @@
}
# Also store a thinned version of the process,
# containing only the final values in the main steps:
- z$ThinPosteriorMu[z$iimain, ] <<- zm$posteriorMu[nrun,]
- z$ThinPosteriorEta[z$iimain, ] <<- zm$posteriorEta[nrun,]
- z$ThinPosteriorSigma[z$iimain, , ] <<- zm$PosteriorSigma[nrun,,]
+ z$ThinPosteriorMu[z$iimain, ] <<- zm$posteriorMu[nrun,]
+ z$ThinPosteriorEta[z$iimain, ] <<- zm$posteriorEta[nrun,]
+ z$ThinPosteriorSigma[z$iimain, , ] <<- zm$PosteriorSigma[nrun,,]
z$ThinBayesAcceptances[z$iimain, ] <<- zm$BayesAcceptances
if ((z$incidentalBasicRates) & (!is.null(zm$thescores)))
{
z$ThinScores[z$iimain, ] <<- zm$thescores
}
- for (group in 1:z$nGroup)
- {
+ for (group in 1:z$nGroup)
+ {
# Drop those elements of the parameters array
# that are NA because for a given group they refer
# to rate parameters of other groups.
z$ThinParameters[z$iimain, group, ] <<-
z$thetaMat[group, !is.na(z$thetaMat[group, ])]
- }
+ }
if (storeAll)
{
z$MHacceptances[start:end, , , ] <<- zm$MHaccepts
z$MHrejections[start:end, , , ] <<- zm$MHrejects
z$MHproportions[start:end, , , ] <<- zm$MHaccepts /
(zm$MHaccepts + zm$MHrejects)
- }
+ }
# Dimension of the following wrong for more than 2 waves:
# z$ThinEtaScores[z$iimain, ,] <<-
# t(zm$ans.last$fra)[z$set2,,drop=FALSE]
# z$ThinEtaHessian[z$iimain, ,] <<-
# as.double(zm$ans.last$dff[z$set2,z$set2,drop=FALSE])
- z$sub <<- z$sub + nrun
- } # end storeData
+ z$sub <<- z$sub + nrun
+ } # end storeData
##@byM internal sienaBayes; for output of vector to screen in nice lines
byM <- function(x,M=10)
@@ -136,9 +137,9 @@
}
##@improveMH internal sienaBayes; find scale factors
- improveMH <- function(tiny=1.0e-15, totruns=100, target=targetMHProb,
+ improveMH <- function(tiny=1.0e-15, totruns=100, target=targetMHProb,
maxiter=20, tolerance=totruns/20, getDocumentation=FALSE)
- {
+ {
# A rough stochastic approximation algorithm.
# Steps when at least two iterations of "totruns" runs resulted,
# for each coordinate of actual, in "desired" acceptances
@@ -149,19 +150,19 @@
return(tt)
}
desired <- trunc(target*totruns)
- iter <- 0
+ iter <- 0
nearGoal <- rep(FALSE, z$nGroup+2)
farFromGoal <- rep(TRUE, z$nGroup+2)
pastSuccesses <- rep(0, z$nGroup+2)
groups <- c(rep(TRUE, z$nGroup), FALSE, FALSE)
cat('improveMH\n')
cat('Desired acceptances', desired, '.\n')
- repeat
- {
- iter <- iter + 1
- MCMCcycle(nrunMH=totruns, nSampVar=1, nSampCons=1,
+ repeat
+ {
+ iter <- iter + 1
+ MCMCcycle(nrunMH=totruns, nSampVar=1, nSampCons=1,
nSampRate=0, change=FALSE, bgain=-1)
- actual <- zm$BayesAcceptances
+ actual <- zm$BayesAcceptances
## actual is a vector of the expected number of acceptances by group
## and for the non-varying parameters, in the MH change of
## parameters out of nrunMHBatches trials
@@ -191,14 +192,14 @@
farMult <- ifelse(actual > 5*desired, 5, farMult)
farMult <- ifelse(5*actual < desired, 0.2, farMult)
actual.previous <- actual
- pastSuccesses <- ifelse(abs(actual - desired) <= tolerance,
+ pastSuccesses <- ifelse(abs(actual - desired) <= tolerance,
pastSuccesses+1, pastSuccesses)
# The calculation of the update:
multiplier <- ifelse(farFromGoal, farMult,
1 + (actual - desired)/
(sqrt(sqrt(iter)*desired*(totruns - desired))))
multiplier <- pmin(pmax(multiplier, 0.1), 10)
- z$scaleFactors <<-
+ z$scaleFactors <<-
z$scaleFactors * multiplier[groups]
z$scaleFactor0 <<-
z$scaleFactor0 * multiplier[z$nGroup+1]
@@ -209,30 +210,30 @@
cat('\nany(is.na(z$scaleFactors)) in improveMH\n')
browser()
}
- cat('\n',iter, '. ', sprintf("%5.1f", actual),
- '\n multipliers ', sprintf("%5.3f", multiplier),
- '\n scaleFactors ', sprintf("%5.3f", z$scaleFactors),
- sprintf("%5.3f", z$scaleFactor0),
- sprintf("%5.3f", z$scaleFactorEta),'\n')
+ cat('\n',iter, '. ', sprintf("%6.1f", actual),
+ '\n multipliers ', sprintf("%6.4f", multiplier),
+ '\n scaleFactors ', sprintf("%6.4f", z$scaleFactors),
+ sprintf("%6.4f", z$scaleFactor0),
+ sprintf("%6.4f", z$scaleFactorEta),'\n')
flush.console()
- if ((min(pastSuccesses) >= 2) || (iter == maxiter))
- {
- break
- }
- if (any(z$scaleFactors < tiny))
- {
+ if ((min(pastSuccesses) >= 2) || (iter == maxiter))
+ {
+ break
+ }
+ if (any(z$scaleFactors < tiny))
+ {
cat('tiny: ', which(z$scaleFactors < tiny),'\n')
cat('This is a problem: scaleFactors too small.\n')
cat('Try a model with fewer randomly varying effects.\n')
stop('Tiny scaleFactors.')
- }
- }
- cat('fine tuning took ', iter, ' iterations.\n')
+ }
+ }
+ cat('fine tuning took ', iter, ' iterations.\n')
#drop cat('scaleFactors:', z$scaleFactors, z$scaleFactor0,
# z$scaleFactorEta, '\n')
flush.console()
- } # end improveMH
+ } # end improveMH
##@MCMCcycle internal sienaBayes; some loops of
@@ -256,7 +257,7 @@
zm$posteriorEta <<- matrix(0, nrow=nrunMH, ncol=z$p2)
zm$PosteriorSigma <<- array( NA, dim=c(nrunMH, z$p1, z$p1))
zm$BayesAcceptances <<- rep(NA, z$nGroup+2)
- zsmall <<- getFromNamespace("makeZsmall", pkgname)(z)
+ zsmall <<- getFromNamespace("makeZsmall", pkgname)(z)
for (i in 1:nrunMH)
{
@@ -426,7 +427,7 @@
restrictProposal <- rep(TRUE, z$p1) # no restriction
}
}
- # do the fandango
+ # do the fandango
## get a multivariate normal with covariance matrix proposalCov
## multiplied by a scale factor which varies between groups
# propose the changes in the varying parameters
@@ -465,11 +466,11 @@
thetaNew[, z$basicRate] <- ifelse(thetaNew[, z$basicRate] < 0.01,
thetaOld[, z$basicRate], thetaNew[, z$basicRate] )
}
- priorOld <- sapply(1:z$nGroup, function(i)
+ priorOld <- sapply(1:z$nGroup, function(i)
{
tmp <- thetaOld[i, z$set1]
use <- !is.na(tmp)
- dmvnorm2(tmp[use], mean=z$muTemp, evs=z$SigmaTemp.evs,
+ dmvnorm2(tmp[use], mean=z$muTemp, evs=z$SigmaTemp.evs,
sigma.inv=z$SigmaTemp.inv) #density
}
)
@@ -477,7 +478,7 @@
{
tmp <- thetaNew[i, z$set1]
use <- !is.na(tmp)
- dmvnorm2(tmp[use], mean=z$muTemp, evs=z$SigmaTemp.evs,
+ dmvnorm2(tmp[use], mean=z$muTemp, evs=z$SigmaTemp.evs,
sigma.inv=z$SigmaTemp.inv) #density
}
)
@@ -514,7 +515,7 @@
}
else
{# used for improveMH, so that it is better to use probabilities
- # for proposals.accept than actual acceptances
+ # for proposals.accept than actual acceptances
z$thetaMat <<- thetaOld
proposals.accept <- exp(pmax(pmin(proposalProbability,0), -100))
# truncation at -100 to avoid rare cases of production of NAs
@@ -554,8 +555,18 @@
logpOld <- getProbabilitiesFromC(z)[[1]]
z$thetaMat <<- thetaNew
logpNew <- getProbabilitiesFromC(z)[[1]]
- proposalProbability <- sum(logpNew - logpOld)
-# cat("eta proposal ", proposalProbability, ", ", logpNew, ", ", logpOld, "\n")
+ if (z$anyPriorEta)
+ {
+ partEtaOld <- thetaOld[1, z$set2][z$set2prior]
+ partEtaNew <- thetaNew[1, z$set2][z$set2prior]
+ logPriorDif <- sum(z$priorPrecEta * (partEtaOld^2 - partEtaNew^2))
+ proposalProbability <- sum(logpNew - logpOld) + logPriorDif
+ }
+ else
+ {
+ proposalProbability <- sum(logpNew - logpOld)
+ }
+# cat("eta proposal ", proposalProbability, ", ", logpNew, ", ", logpOld, "\n")
constantPar.accept <- (log(runif(1)) < proposalProbability)
if (is.na(constantPar.accept)){constantPar.accept <- FALSE}
if (change)
@@ -568,7 +579,7 @@
}
else
{# used for improveMH, so that it is better to use probabilities
- # for z$accept than actual acceptances
+ # for z$accept than actual acceptances
z$acceptEta <<- exp(min(proposalProbability, 0))
z$thetaMat <<- thetaOld
}
@@ -651,7 +662,7 @@
# The update needs to be applied only to the unmasked part,
# because the masked part plays no role.
z$muTemp[updated] <<- z$muTemp[updated] -
- bgain * correct * (z$muTemp[updated] - muhat)
+ bgain * correct * (z$muTemp[updated] - muhat)
# or with a matrix:
# z$muTemp[updated] - bgain *
# as.vector(z$dmuinv %*% (z$muTemp[updated] - muhat))
@@ -671,7 +682,7 @@
etaScores <- rowSums(scores[z$set2,,drop=FALSE])
#browser()
#cat("difference = ", max(abs(scores.ans.last + (scores))))
-#cat(" getProbabilitiesFromC\n")
+#cat(" getProbabilitiesFromC\n")
#browser()
#prc <- getProbabilitiesFromC(z, getScores=TRUE)
#str(prc,1)
@@ -736,7 +747,7 @@
{
cat("correctMatrix: Eigenvalues Sigma = ", sort(esv), "\n")
# tol <- dd*max(abs(esv))*.Machine$double.eps
-# delt <- 2*tol # factor two is just to make sure the resulting
+# delt <- 2*tol # factor two is just to make sure the resulting
# matrix passes all numerical tests of positive definiteness
# TS: I replace this delta by a larger number
tau <- pmax(0, delt - esv)
@@ -757,7 +768,7 @@
thetaMean <- rep(NA, z$pp)
for (group in 1:z$nGroup)
{
- thetaMean[z$ratePositions[[group]]] <- colMeans(
+ thetaMean[z$ratePositions[[group]]] <- colMeans(
z$ThinParameters[, group, !z$generalParametersInGroup,
drop=FALSE], na.rm=TRUE)
}
@@ -810,9 +821,9 @@
}
cr + 0.05*diag(diag(cr), nrow = dim(x)[2])
}
- ## ###################################### ##
- ## start of function sienaBayes() proper ##
- ## ###################################### ##
+ ## ###################################### ##
+ ## start of function sienaBayes() proper ##
+ ## ###################################### ##
if (getDocumentation != FALSE)
{
if (getDocumentation == TRUE)
@@ -846,14 +857,14 @@
z <- initializeBayes(data, effects, algo, nbrNodes,
initgainGlobal=initgainGlobal,
initgainGroupwise=initgainGroupwise, initML=initML,
+ priorSigEta=priorSigEta,
priorMu=priorMu, priorSigma=priorSigma,
priorDf=priorDf, priorKappa=priorKappa,
priorRatesFromData=priorRatesFromData,
frequentist=frequentist,
incidentalBasicRates=incidentalBasicRates,
reductionFactor=reductionFactor, gamma=gamma,
- delta=delta, nmain=nmain, nwarm=nwarm,
- nImproveMH=nImproveMH,
+ delta=delta, nmain=nmain, nprewarm=nprewarm, nwarm=nwarm,
lengthPhase1=lengthPhase1, lengthPhase3=lengthPhase3,
prevAns=prevAns, usePrevOnly=usePrevOnly,
silentstart=silentstart, clusterType=clusterType)
@@ -874,16 +885,18 @@
z$sub <- 0
zm <- list()
zsmall <- list()
+
+ z$nImproveMH <- nImproveMH
if (nrunMHBatches >= 2)
- {
- z$nrunMHBatches <- nrunMHBatches
- }
+ {
+ z$nrunMHBatches <- nrunMHBatches
+ }
else
- {
- cat("nrunMHBatches increased to 2.\n")
- z$nrunMHBatches <- 2
- nrunMHBatches <- 2
- }
+ {
+ cat("nrunMHBatches increased to 2.\n")
+ z$nrunMHBatches <- 2
+ nrunMHBatches <- 2
+ }
z$nSampVarying <- min(nSampVarying, 1)
if (nSampVarying < 1)
{
@@ -903,15 +916,33 @@
class(z) <- "sienaBayesFit"
z$correctSigma <- 0
- # Determine multiplicative constants for proposal distribution
ctime1 <- proc.time()
cat(round((ctime1-ctime)[3]),' seconds elapsed\n')
+
+ # Pre-warming phase
+ bgain <- initfgain
+ for (ii in 1:nprewarm)
+ {
+ MCMCcycle(nrunMH=nrunMHBatches, nSampVar=nSampVarying,
+ nSampCons=nSampConst, nSampRate=nSampRates,
+ bgain=bgain)
+ cat('Pre-warming step',ii,'(',nprewarm,')\n')
+ cat("Accepts ",sum(zm$BayesAcceptances),"/",
+ z$nGroup*nrunMHBatches,"\n")
+ flush.console()
+ }
+ print('end of pre-warming')
+ ctime3p<- proc.time()
+ cat('pre-warming took', round((ctime3p-ctime1)[3]),'seconds elapsed.\n')
+ flush.console()
+
+ # Determine multiplicative constants for proposal distribution
improveMH(totruns=nImproveMH, target=targetMHProb)
ctime2<- proc.time()
cat('improveMH', round((ctime2-ctime1)[3]),' seconds elapsed for improveMH.\n')
flush.console()
}
- else # (!is.null(prevBayes))
+ else # (!is.null(prevBayes))
{
if (inherits(prevBayes, "sienaBayesFit"))
{
@@ -946,6 +977,11 @@
z$nSampConst <- nSampConst
cat("nSampConst changed to ", nSampConst, ".\n")
}
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/rsiena -r 337
More information about the Rsiena-commits
mailing list