[Rsiena-commits] r153 - in pkg: RSiena RSiena/R RSiena/inst/doc RSiena/man RSiena/src RSiena/src/model RSiena/src/model/effects RSiena/src/model/ml RSiena/src/model/variables RSienaTest RSienaTest/R RSienaTest/doc RSienaTest/inst/doc RSienaTest/inst/examples RSienaTest/man RSienaTest/src RSienaTest/src/model RSienaTest/src/model/effects RSienaTest/src/model/effects/generic RSienaTest/src/model/ml RSienaTest/src/model/variables
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat Jun 4 15:03:12 CEST 2011
Author: ripleyrm
Date: 2011-06-04 15:03:10 +0200 (Sat, 04 Jun 2011)
New Revision: 153
Added:
pkg/RSiena/inst/doc/RSiena_Manual.pdf
pkg/RSienaTest/doc/RSiena_Manual.tex
pkg/RSienaTest/inst/doc/RSiena_Manual.pdf
pkg/RSienaTest/src/model/effects/generic/GwespFunction.cpp
pkg/RSienaTest/src/model/effects/generic/GwespFunction.h
Removed:
pkg/RSiena/inst/doc/s_man400.pdf
pkg/RSienaTest/doc/s_man400.tex
pkg/RSienaTest/inst/doc/s_man400.pdf
pkg/RSienaTest/src/model/effects/GwespBBEffect.cpp
pkg/RSienaTest/src/model/effects/GwespBBEffect.h
pkg/RSienaTest/src/model/effects/GwespBFEffect.cpp
pkg/RSienaTest/src/model/effects/GwespBFEffect.h
pkg/RSienaTest/src/model/effects/GwespEffect.cpp
pkg/RSienaTest/src/model/effects/GwespEffect.h
pkg/RSienaTest/src/model/effects/GwespFBEffect.cpp
pkg/RSienaTest/src/model/effects/GwespFBEffect.h
pkg/RSienaTest/src/model/effects/GwespFFEffect.cpp
pkg/RSienaTest/src/model/effects/GwespFFEffect.h
pkg/RSienaTest/src/model/effects/GwespRREffect.cpp
pkg/RSienaTest/src/model/effects/GwespRREffect.h
Modified:
pkg/RSiena/DESCRIPTION
pkg/RSiena/R/Sienatest.r
pkg/RSiena/R/bayes.r
pkg/RSiena/R/effects.r
pkg/RSiena/R/effectsMethods.r
pkg/RSiena/R/initializeFRAN.r
pkg/RSiena/R/maxlikec.r
pkg/RSiena/R/phase3.r
pkg/RSiena/R/print07Report.r
pkg/RSiena/R/siena01.r
pkg/RSiena/R/siena07.r
pkg/RSiena/R/sienaeffects.r
pkg/RSiena/R/sienaprint.r
pkg/RSiena/R/simstatsc.r
pkg/RSiena/changeLog
pkg/RSiena/man/RSiena-package.Rd
pkg/RSiena/man/sienaNet.Rd
pkg/RSiena/man/simstats0c.Rd
pkg/RSiena/src/model/EpochSimulation.cpp
pkg/RSiena/src/model/EpochSimulation.h
pkg/RSiena/src/model/Model.cpp
pkg/RSiena/src/model/Model.h
pkg/RSiena/src/model/StatisticCalculator.cpp
pkg/RSiena/src/model/StatisticCalculator.h
pkg/RSiena/src/model/effects/BehaviorEffect.cpp
pkg/RSiena/src/model/effects/BehaviorEffect.h
pkg/RSiena/src/model/effects/NetworkEffect.cpp
pkg/RSiena/src/model/effects/NetworkEffect.h
pkg/RSiena/src/model/effects/OutTruncEffect.cpp
pkg/RSiena/src/model/effects/StructuralRateEffect.cpp
pkg/RSiena/src/model/effects/StructuralRateEffect.h
pkg/RSiena/src/model/ml/BehaviorChange.cpp
pkg/RSiena/src/model/ml/BehaviorChange.h
pkg/RSiena/src/model/ml/Chain.cpp
pkg/RSiena/src/model/ml/Chain.h
pkg/RSiena/src/model/ml/MLSimulation.cpp
pkg/RSiena/src/model/ml/MLSimulation.h
pkg/RSiena/src/model/ml/MiniStep.cpp
pkg/RSiena/src/model/ml/MiniStep.h
pkg/RSiena/src/model/ml/NetworkChange.cpp
pkg/RSiena/src/model/ml/NetworkChange.h
pkg/RSiena/src/model/variables/BehaviorVariable.cpp
pkg/RSiena/src/model/variables/BehaviorVariable.h
pkg/RSiena/src/model/variables/DependentVariable.cpp
pkg/RSiena/src/model/variables/DependentVariable.h
pkg/RSiena/src/model/variables/NetworkVariable.cpp
pkg/RSiena/src/model/variables/NetworkVariable.h
pkg/RSiena/src/siena07internals.cpp
pkg/RSiena/src/siena07models.cpp
pkg/RSiena/src/siena07setup.cpp
pkg/RSiena/src/siena07setup.h
pkg/RSiena/src/siena07utilities.cpp
pkg/RSiena/src/siena07utilities.h
pkg/RSienaTest/DESCRIPTION
pkg/RSienaTest/R/Sienatest.r
pkg/RSienaTest/R/bayes.r
pkg/RSienaTest/R/effects.r
pkg/RSienaTest/R/effectsMethods.r
pkg/RSienaTest/R/initializeFRAN.r
pkg/RSienaTest/R/maxlikec.r
pkg/RSienaTest/R/phase3.r
pkg/RSienaTest/R/print07Report.r
pkg/RSienaTest/R/siena01.r
pkg/RSienaTest/R/siena07.r
pkg/RSienaTest/R/sienaeffects.r
pkg/RSienaTest/R/sienaprint.r
pkg/RSienaTest/R/simstatsc.r
pkg/RSienaTest/changeLog
pkg/RSienaTest/doc/RSienaDeveloper.tex
pkg/RSienaTest/inst/examples/algorithms.r
pkg/RSienaTest/inst/examples/runalg.r
pkg/RSienaTest/man/RSiena-package.Rd
pkg/RSienaTest/man/sienaNet.Rd
pkg/RSienaTest/man/simstats0c.Rd
pkg/RSienaTest/src/model/EpochSimulation.cpp
pkg/RSienaTest/src/model/EpochSimulation.h
pkg/RSienaTest/src/model/Model.cpp
pkg/RSienaTest/src/model/Model.h
pkg/RSienaTest/src/model/StatisticCalculator.cpp
pkg/RSienaTest/src/model/StatisticCalculator.h
pkg/RSienaTest/src/model/effects/AllEffects.h
pkg/RSienaTest/src/model/effects/BehaviorEffect.cpp
pkg/RSienaTest/src/model/effects/BehaviorEffect.h
pkg/RSienaTest/src/model/effects/EffectFactory.cpp
pkg/RSienaTest/src/model/effects/NetworkEffect.cpp
pkg/RSienaTest/src/model/effects/NetworkEffect.h
pkg/RSienaTest/src/model/effects/OutTruncEffect.cpp
pkg/RSienaTest/src/model/effects/StructuralRateEffect.cpp
pkg/RSienaTest/src/model/effects/StructuralRateEffect.h
pkg/RSienaTest/src/model/ml/BehaviorChange.cpp
pkg/RSienaTest/src/model/ml/BehaviorChange.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/ml/MiniStep.cpp
pkg/RSienaTest/src/model/ml/MiniStep.h
pkg/RSienaTest/src/model/ml/NetworkChange.cpp
pkg/RSienaTest/src/model/ml/NetworkChange.h
pkg/RSienaTest/src/model/variables/BehaviorVariable.cpp
pkg/RSienaTest/src/model/variables/BehaviorVariable.h
pkg/RSienaTest/src/model/variables/DependentVariable.cpp
pkg/RSienaTest/src/model/variables/DependentVariable.h
pkg/RSienaTest/src/model/variables/NetworkVariable.cpp
pkg/RSienaTest/src/model/variables/NetworkVariable.h
pkg/RSienaTest/src/siena07internals.cpp
pkg/RSienaTest/src/siena07models.cpp
pkg/RSienaTest/src/siena07setup.cpp
pkg/RSienaTest/src/siena07setup.h
pkg/RSienaTest/src/siena07utilities.cpp
pkg/RSienaTest/src/siena07utilities.h
Log:
Creation effects, ML more work in progress: fixes for bipartite and constraints.
Modified: pkg/RSiena/DESCRIPTION
===================================================================
--- pkg/RSiena/DESCRIPTION 2011-05-29 18:42:29 UTC (rev 152)
+++ pkg/RSiena/DESCRIPTION 2011-06-04 13:03:10 UTC (rev 153)
@@ -1,8 +1,8 @@
Package: RSiena
Type: Package
Title: Siena - Simulation Investigation for Empirical Network Analysis
-Version: 1.0.12.152
-Date: 2011-05-29
+Version: 1.0.12.153
+Date: 2011-06-04
Author: Various
Depends: R (>= 2.10.0), xtable
Imports: Matrix
Modified: pkg/RSiena/R/Sienatest.r
===================================================================
--- pkg/RSiena/R/Sienatest.r 2011-05-29 18:42:29 UTC (rev 152)
+++ pkg/RSiena/R/Sienatest.r 2011-06-04 13:03:10 UTC (rev 153)
@@ -16,53 +16,59 @@
## can be obtained via svd(data) (which is stored in z$sf). Square of ratio
## of smallest to largest singular value.
Report('Instability Analysis\n')
- # pp<- length(z$diver)
- constant<- z$diver
- test<- z$test
- covtheta<- z$covtheta
- covZ<- z$msf
- covth<- covtheta[!(test|constant),!(test|constant)]
- covth<- MatrixNorm(covth)
- eigenv<- eigen(covth,symmetric=TRUE)$values
- ma<- max(eigenv)
- mi<- min(eigenv)
+ constant <- z$diver
+ test <- z$test
+ covtheta <- z$covtheta
+ covZ <- z$msf
+ covth <- covtheta[!(test|constant), !(test|constant)]
+ covth <- MatrixNorm(covth)
+ eigenv <- eigen(covth,symmetric=TRUE)$values
+ ma <- max(eigenv)
+ mi <- min(eigenv)
if (mi!=0)
- cond.n <- ma/mi
- Report('Instability analysis\n',lf)
- Report('--------------------\n\n',lf)
- Report('Variance-covariance matrix of parameter estimates',lf)
+ {
+ cond.n <- ma / mi
+ }
+ Report('Instability analysis\n', lf)
+ Report('--------------------\n\n', lf)
+ Report('Variance-covariance matrix of parameter estimates', lf)
##if (global boolean1 )
## Report(' (without coordinates that are kept constant):\n',lf)
##else
- Report(c(':\n\nCondition number = ',format(cond.n,width=4,nsmall=4,digits=1),
+ Report(c(':\n\nCondition number = ',format(cond.n, width=4, nsmall=4,
+ digits=1),
' \n\n'),sep='',lf)
- Report(c('Eigen Values ',format(eigenv,width=6,nsmall=6,digits=1)),lf)
+ Report(c('Eigen Values ', format(eigenv, width=6, nsmall=6, digits=1)), lf)
Report('\n\n',lf)
- covZ<- MatrixNorm(covZ)
- eigenvZ<-eigen(covZ,symmetric=TRUE)$values
- ma<- max(eigenvZ)
- mi<- min(eigenvZ)
+ covZ <- MatrixNorm(covZ)
+ eigenvZ <-eigen(covZ, symmetric=TRUE)$values
+ ma <- max(eigenvZ)
+ mi <- min(eigenvZ)
if (mi!=0)
- cond.n <- ma/mi
+ {
+ cond.n <- ma / mi
+ }
Report('Variance-covariance matrix of X',lf)
- Report(c(':\n\nCondition number = ',format(cond.n,width=4,nsmall=4,digits=1),
- ' \n\n'),sep='',lf)
- Report(c('Eigen Values ',format(eigenvZ,width=6,nsmall=6,digits=1)),lf)
+ Report(c(':\n\nCondition number = ', format(cond.n, width=4, nsmall=4,
+ digits=1),
+ ' \n\n'), sep='', lf)
+ Report(c('Eigen Values ', format(eigenvZ, width=6, nsmall=6, digits=1)), lf)
Report(c('\n\n',date(),'\n'),sep='',lf)
- mysvd<- svd(z$sf)$d
- ma<- max(mysvd)
- mi<- min(mysvd)
- cond.n<- (ma/mi)^2
- Report(c(':\n\nCondition number2 = ',format(cond.n,width=4,nsmall=4,digits=1),
- ' \n\n'),sep='',lf)
- Report(c('Singular Values ',format(mysvd,width=6,nsmall=6,digits=1)),lf)
- Report(c('\n\n',date(),'\n'),sep='',lf)
+ mysvd <- svd(z$sf)$d
+ ma <- max(mysvd)
+ mi <- min(mysvd)
+ cond.n <- (ma/mi)^2
+ Report(c(':\n\nCondition number2 = ',
+ format(cond.n, width=4, nsmall=4, digits=1),
+ ' \n\n'), sep='', lf)
+ Report(c('Singular Values ', format(mysvd, width=6, nsmall=6, digits=1)), lf)
+ Report(c('\n\n', date(), '\n'), sep='', lf)
}
##@MatrixNorm siena07 Not currently used. May be incorrect.
-MatrixNorm<- function(mat)
+MatrixNorm <- function(mat)
{
- tmp<- apply(mat,2,function(x)x/sqrt(crossprod(x)))
+ tmp <- apply(mat, 2, function(x) x / sqrt(crossprod(x)))
##or sweep(mat,2,apply(mat,2,function(x)x/sqrt(crossprod(x))
tmp
}
@@ -70,27 +76,33 @@
TestOutput <- function(z, x)
{
testn <- sum(z$test)
- # browser()
+ ## browser()
if (testn)
{
if (x$maxlike)
+ {
Heading(2, outf,'Score test <c>')
+ }
else
+ {
Heading(2, outf, 'Generalised score test <c>')
- Report('Testing the goodness-of-fit of the model restricted by\n',outf)
- j<- 0
+ }
+ Report('Testing the goodness-of-fit of the model restricted by\n', outf)
+ j <- 0
for (k in 1:z$pp)
+ {
if (z$test[k])
{
- j<- j+1
- Report(c(" (",j,") ",
+ j <- j + 1
+ Report(c(" (", j, ") ",
format(paste(z$requestedEffects$type[k], ": ",
z$requestedEffects$effectName[k],
- sep=""),
- width=50), " = ",
- sprintf("%8.4f",z$theta[k]),"\n"),
+ sep=""),
+ width=50), " = ",
+ sprintf("%8.4f", z$theta[k]), "\n"),
sep = "", outf)
}
+ }
Report('_________________________________________________\n',outf)
Report(' ',outf)
Report(' \n',outf)
@@ -118,22 +130,29 @@
' d.f. = 1 p-value '), sep = '', outf)
pvalue<- 1-pchisq(z$testresult[k],1)
if (pvalue < 0.0001)
+ {
Report('< 0.0001\n',outf)
+ }
else
+ {
Report(c('= ', sprintf("%8.4f", pvalue), '\n'), sep = '',
outf)
+ }
Report(c(' - one-sided (normal variate): ',
sprintf("%8.4f", z$testresulto[k])), sep = '', outf)
- if (k<j)
+ if (k < j)
+ {
Report('\n\n',outf)
+ }
}
}
- Report(' \n_________________________________________________\n\n',outf)
- Report('One-step estimates: \n\n',outf)
+ Report(' \n_________________________________________________\n\n',
+ outf)
+ Report('One-step estimates: \n\n', outf)
for (i in 1 : z$pp)
{
- onestepest<- z$oneStep[i]+z$theta[i]
- Report(c(format(paste(z$requestedEffects$type[i],': ',
+ onestepest <- z$oneStep[i] + z$theta[i]
+ Report(c(format(paste(z$requestedEffects$type[i], ': ',
z$requestedEffects$effectName[i], sep = ''),
width=50),
sprintf("%8.4f", onestepest), '\n'), sep = '', outf)
@@ -141,10 +160,11 @@
Report('\n',outf)
}
}
+
##@ScoreTest siena07 Do score tests
ScoreTest<- function(pp, dfra, msf, fra, test, maxlike)
{
- testresult<- rep(NA, pp) ##for chisq per parameter
+ testresult <- rep(NA, pp) ##for chisq per parameter
testresulto <- rep(NA, pp) ##for one-sided tests per parameter
##first the general one
ans <- EvaluateTestStatistic(maxlike, test, dfra, msf, fra)
@@ -163,7 +183,7 @@
{
if (test[i])
{
- k <- k+1
+ k <- k + 1
use[i] <- TRUE
ans <- EvaluateTestStatistic(maxlike, test[use], dfra[use, use],
msf[use, use], fra[use])
Modified: pkg/RSiena/R/bayes.r
===================================================================
--- pkg/RSiena/R/bayes.r 2011-05-29 18:42:29 UTC (rev 152)
+++ pkg/RSiena/R/bayes.r 2011-06-04 13:03:10 UTC (rev 153)
@@ -10,29 +10,32 @@
## ****************************************************************************/
##@bayes algorithm fit a Bayesian model
bayes <- function(data, effects, model, nwarm=100, nmain=100, nrunMHBatches=20,
- nrunMH=100, save=TRUE)
+ nrunMH=100, save=TRUE, nbrNodes=1, seed=1)
{
- createStores <- function()
+ 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=nmain * nrunMHBatches,
- ncol=sum(basicRate))
- z$lambdas <- matrix(NA, nrow=nmain * nrunMHBatches,
- ncol=sum(basicRate))
- z$betas <- matrix(NA, nrow=nmain * nrunMHBatches, ncol=sum(!basicRate))
- z$candidates <- matrix(NA, nrow=nmain * nrunMHBatches,
- ncol=sum(!basicRate))
- z$acceptances <- rep(NA, nmain * nrunMHBatches)
- z$MHacceptances <- matrix(NA, nrow=nmain * nrunMHBatches, ncol=7)
- z$MHrejections <- matrix(NA, nrow=nmain * nrunMHBatches , ncol=7)
- z$MHproportions <- matrix(NA, nrow=nmain * nrunMHBatches, ncol=7)
+ 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)
z
}
storeData <- function()
{
start <- z$sub + 1
nrun <- nrow(z$parameters)
+ basicRate <- z$basicRate
+ npar <- length(z$theta)
end <- start + nrun - 1
z$accepts <- as.logical(z$accepts)
z$acceptances[start:end] <- z$accepts
@@ -43,7 +46,6 @@
z$parameters[, !basicRate] <-
carryForward(z$parameters[, !basicRate], z$betas[1:end,])
z$betas[start:end, ] <- z$parameters[, !basicRate]
-
z$posteriorTot <- z$posteriorTot + colSums(z$parameters)
for (i in npar)
{
@@ -56,7 +58,9 @@
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)
@@ -65,6 +69,10 @@
{
parameters <- rbind(betas[(nbeta - npar), ], parameters)
}
+ else
+ {
+ parameters <- rbind(z$theta[!z$basicRate], parameters)
+ }
parameters <-
apply(parameters, 2, function(x)
{
@@ -77,16 +85,8 @@
length(x) - x.pos[length(x.pos)] + 1))]
}
)
- if (npar < nbeta)
- {
- parameters[-1, ]
- }
- else
- {
- parameters
- }
+ parameters[-1, ]
}
-
improveMH <- function(z, x, tiny=1.0e-15, desired=40, maxiter=100,
tolerance=15)
{
@@ -119,7 +119,7 @@
iter <- iter + 1
z <- MCMCcycle(z, nrunMH=1, nrunMHBatches=100)
actual <- z$BayesAcceptances ## acceptances
- ans <- rescaleCGD(100)
+ ans <- rescaleCGD(100 * (z$observations - 1))
z$scaleFactor <- z$scaleFactor * ans
if (success | iter == maxiter)
{
@@ -127,7 +127,7 @@
}
if (z$scaleFactor < tiny)
{
- cat('calefactor < tiny\n')
+ cat('scalefactor < tiny\n')
browser()
}
}
@@ -135,9 +135,19 @@
z$scaleFactor, '\n')
z
}
+ ## ################################
+ ## start of function proper
+ ## ################################
- ## initialise
- Report(openfiles=TRUE, type="n") #initialise with no file
+ ## 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 <- createStores(z)
+
+ z$sub <- 0
+
+ z <- improveMH(z)
+
if (!save)
{
require(lattice)
@@ -146,39 +156,7 @@
dev.new()
acceptsplot = dev.cur()
}
- z <- NULL
- z$FinDiff.method <- FALSE
- z$maxlike <- TRUE
- model$maxlike <- TRUE
- model$FRANname <- "maxlikec"
- z$print <- FALSE
- z$int <- 1
- z$int2 <- 1
- model$cconditional <- FALSE
- if (!is.null(model$randomSeed))
- {
- 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,
- returnDeps=TRUE)
- ##if (useCluster)
- ## cl <- makeCluster(rep("localhost", 2), type = "SOCK")
- is.batch(TRUE)
-
- WriteOutTheta(z)
-
- npar <- length(z$theta)
- basicRate <- z$effects$basicRate
- ##iter <- 0
- z$numm <- 20
- z$scaleFactor <- 1
- z <- createStores()
- z$sub <- 0
- z <- improveMH(z)
-
for (ii in 1:nwarm)
{
z <- MCMCcycle(z, nrunMH=4, nrunMHBatches=20)
@@ -187,16 +165,20 @@
for (ii in 1:nmain)
{
- z <- MCMCcycle(z, nrunMH=100, nrunMHBatches=20)
+ 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')
dev.set(thetaplot)
- thetadf <- data.frame(z$lambdas, z$betas)
- acceptsdf <- data.frame(z$MHproportions[, 1:5],
+ basicRate <- z$basicRate
+ lambdas <- z$lambdas
+ lambdas[lambdas == 0] <- NA
+ thetadf <- data.frame(lambdas, z$betas)
+ acceptsdf <- data.frame(z$MHproportions,
z$acceptances)
lambdaNames <- paste(z$effects$name[basicRate],
z$effects$shortName[basicRate],
@@ -204,20 +186,32 @@
z$effects$group[basicRate], sep=".")
betaNames <- paste(z$effects$name[!basicRate],
z$effects$shortName[!basicRate], sep=".")
- names(thetadf) <- c(lambdaNames, betaNames)
+ names(thetadf) <- make.names(c(lambdaNames, betaNames),
+ unique=TRUE)
names(acceptsdf) <- c("InsDiag", "CancDiag", "Permute", "InsPerm",
- "CancPerm", #"Missing",
+ "DelPerm", "InsMissing", "DelMissing",
"BayesAccepts")
varnames <- paste(names(thetadf), sep="", collapse= " + ")
varcall <- paste("~ ", varnames, sep="", collapse="")
print(histogram(as.formula(varcall), data=thetadf, scales="free",
- outer=TRUE, breaks=NULL))
+ outer=TRUE, breaks=NULL, type="density",
+ panel=function(x, ...)
+ {
+ panel.histogram(x, ...)
+ 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"),
- outer=TRUE, breaks=NULL, type="density"))
+ outer=TRUE, breaks=NULL, type="density",
+ panel=function(x, ...)
+ {
+ panel.histogram(x, ...)
+ panel.densityplot(x, darg=list(na.rm=TRUE), ...)
+ }))
}
}
z
@@ -225,17 +219,45 @@
MCMCcycle <- function(z, nrunMH, nrunMHBatches)
{
- ## this function assumes only one period at the moment, but would deal with
- ## multiple periods locally, calling MCMCcycle for each. Not yet available
- ## in C, due to need to keep multiple chains in existence. I have not
- ## decided how to do it yet!
- period <- 1
- group <- 1
- f <- FRANstore()
- ans <- .Call("MCMCcycle", PACKAGE=pkgname, f$pData, f$pModel,
- f$myeffects, as.integer(period),
- as.integer(group),
- z$scaleFactor, 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)
+ {
+ 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)
+ {
+ 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]]
@@ -243,8 +265,66 @@
z$shapes <- ans[[3]]
z$MHaccepts <- ans[[4]]
z$MHrejects <- ans[[5]]
+ z$chains <- ans[[6]]
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)
+{
+ ## initialise
+ set.seed(seed)
+ Report(openfiles=TRUE, type="n") #initialise with no file
+ z <- NULL
+ z$FinDiff.method <- FALSE
+ z$maxlike <- TRUE
+ model$maxlike <- TRUE
+ model$FRANname <- "maxlikec"
+ z$print <- FALSE
+ z$int <- 1
+ z$int2 <- nbrNodes
+ model$cconditional <- FALSE
+ if (!is.null(model$randomSeed))
+ {
+ 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)
+ is.batch(TRUE)
+
+ WriteOutTheta(z)
+
+ if (nbrNodes > 1)
+ {
+ require(snow)
+ require(rlecuyer)
+ clusterString <- rep("localhost", nbrNodes)
+ z$cl <- makeCluster(clusterString, type = "SOCK",
+ outfile = "cluster.out")
+ clusterCall(z$cl, library, pkgname, character.only = TRUE)
+ clusterCall(z$cl, storeinFRANstore, FRANstore())
+ clusterCall(z$cl, FRANstore)
+ clusterCall(z$cl, initializeFRAN, z, model,
+ 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
+}
Modified: pkg/RSiena/R/effects.r
===================================================================
--- pkg/RSiena/R/effects.r 2011-05-29 18:42:29 UTC (rev 152)
+++ pkg/RSiena/R/effects.r 2011-06-04 13:03:10 UTC (rev 153)
@@ -51,9 +51,10 @@
if (!all(is.na(effects$endowment)))
{
neweffects <- effects[rep(1:nn,
- times=(1 + as.numeric(effects$endowment))), ]
- neweffects$type <- unlist(lapply(effects$endowment, function(x) if (x)
- c('eval', 'endow') else 'eval'))
+ times=(1 + 2 * as.numeric(effects$endowment))), ]
+ neweffects$type <-
+ unlist(lapply(effects$endowment, function(x) if (x)
+ c('eval', 'endow', 'creation') else 'eval'))
effects <- neweffects
nn <- nrow(effects)
}
@@ -306,10 +307,13 @@
#objEffects <- createObjEffectList(objEffects, varname)
#rateEffects <- createRateEffectList(rateEffects, varname)
- ## replace the text for endowment effects
+ ## replace the text for endowment and creation effects
tmp <- objEffects$functionName[objEffects$type =='endow']
tmp <- paste('Lost ties:', tmp)
objEffects$functionName[objEffects$type == 'endow'] <- tmp
+ tmp <- objEffects$functionName[objEffects$type =='creation']
+ tmp <- paste('New ties:', tmp)
+ objEffects$functionName[objEffects$type == 'creation'] <- tmp
## get starting values
starts <- getNetworkStartingVals(depvar)
@@ -504,10 +508,13 @@
rateEffects[1:noPeriods, 'initialValue'] <- starts$startRate
rateEffects$basicRate[1:observations] <- TRUE
- ## alter the text for endowment names
+ ## alter the text for endowment and creation names
objEffects$effectName[objEffects$type == 'endow'] <-
sub('behavior', 'dec. beh.',
objEffects$effectName[objEffects$type == 'endow'])
+ objEffects$effectName[objEffects$type == 'creation'] <-
+ sub('behavior', 'inc. beh.',
+ objEffects$effectName[objEffects$type == 'creation'])
list(effects = rbind(rateEffects = rateEffects,
objEffects = objEffects), starts=starts)
@@ -666,6 +673,9 @@
objEffects$functionName[objEffects$type == 'endow'] <-
paste('Lost ties:',
objEffects$functionName[objEffects$type =='endow'])
+ objEffects$functionName[objEffects$type == 'creation'] <-
+ paste('New ties:',
+ objEffects$functionName[objEffects$type =='creation'])
## get starting values
starts <- getBipartiteStartingVals(depvar)
Modified: pkg/RSiena/R/effectsMethods.r
===================================================================
--- pkg/RSiena/R/effectsMethods.r 2011-05-29 18:42:29 UTC (rev 152)
+++ pkg/RSiena/R/effectsMethods.r 2011-06-04 13:03:10 UTC (rev 153)
@@ -44,7 +44,7 @@
nDependents <- length(unique(x$name))
userSpecifieds <- x$shortName[x$include] %in% c("unspInt", "behUnspInt")
endowments <- !x$type[x$include] %in% c("rate", "eval")
- timeDummies <- !x$timeDummies[x$include] == ","
+ timeDummies <- !x$timeDummy[x$include] == ","
specs <- x[, c("name", "effectName", "include", "fix", "test",
"initialValue", "parm")]
if (includeOnly)
@@ -61,7 +61,7 @@
}
if (any(timeDummies))
{
- specs <- cbind(specs, type=x[x$include, "timeDummies"])
+ specs <- cbind(specs, timeDummy=x[x$include, "timeDummy"])
}
if (any(userSpecifieds))
{
Modified: pkg/RSiena/R/initializeFRAN.r
===================================================================
--- pkg/RSiena/R/initializeFRAN.r 2011-05-29 18:42:29 UTC (rev 152)
+++ pkg/RSiena/R/initializeFRAN.r 2011-06-04 13:03:10 UTC (rev 153)
@@ -167,7 +167,6 @@
}
}
types <- sapply(data[[1]]$depvars, function(x) attr(x, 'type'))
- ##nets <- sum(types != "behavior")
## now check if conditional estimation is OK and copy to z if so
z$cconditional <- FALSE
if (x$cconditional)
@@ -489,43 +488,67 @@
if (x$maxlike)
{
- ## set up chains and do initial steps
- simpleRates <- TRUE
+ simpleRates <- TRUE ## get a sensible test for this soon!
+ if (any(!z$effects$basicRate & z$effects$type =="rate"))
+ {
+ # browser()
+ simpleRates <- FALSE
+ }
+ z$simpleRates <- simpleRates
+ if (!initC)
+ {
+ ## set up chains and do initial steps
- types <- attr(f, "types")
- nbrNonMissNet <- attr(f, "numberNonMissingNetwork")
- nbrMissNet <- attr(f, "numberMissingNetwork")
- nbrNonMissBeh <- attr(f, "numberNonMissingBehavior")
- nbrMissBeh <- attr(f, "numberMissingBehavior")
+ types <- attr(f, "types")
+ nbrNonMissNet <- attr(f, "numberNonMissingNetwork")
+ nbrMissNet <- attr(f, "numberMissingNetwork")
+ nbrNonMissBeh <- attr(f, "numberNonMissingBehavior")
+ nbrMissBeh <- attr(f, "numberMissingBehavior")
- z$prmin <- nbrMissNet/ (nbrMissNet + nbrNonMissNet)
- if (sum(nbrMissBeh + nbrNonMissBeh) > 0)
- {
- z$prmib <- nbrMissBeh/ (nbrMissBeh + nbrNonMissBeh)
+ if (sum(nbrMissNet + nbrNonMissNet) > 0)
+ {
+ z$prmin <- nbrMissNet/ (nbrMissNet + nbrNonMissNet)
+ }
+ else
+ {
+ z$prmin <- rep(0, length(nbrMissNet))
+ }
+ if (sum(nbrMissBeh + nbrNonMissBeh) > 0)
+ {
+ z$prmib <- nbrMissBeh/ (nbrMissBeh + nbrNonMissBeh)
+ }
+ else
+ {
+ z$prmib <- rep(0, length(nbrMissBeh))
+ }
+ ## cat (z$prmin, z$prmib, '\n')
+ z$probs <- c(x$pridg, x$prcdg, x$prper, x$pripr, x$prdpr, x$prirms,
+ x$prdrms)
+ cat(z$probs,'\n')
+ ans <- .Call("mlMakeChains", PACKAGE=pkgname, pData, pModel,
+ simpleRates, z$probs, z$prmin, z$prmib,
+ x$minimumPermutationLength,
+ x$maximumPermutationLength,
+ x$initialPermutationLength)
+
+ f$minimalChain <- ans[[1]]
+ f$chain <- ans[[2]]
+ f$simpleRates <- simpleRates
+ ## print(nrow(ans[[1]][[1]]))
+ ## print(nrow(ans[[2]][[1]]))
+ ## browser()
}
- else
+ else ## set up the initial chains in the sub processes
{
- z$prmib <- rep(0, length(nbrMissBeh))
- }
- ## cat (z$prmin, z$prmib, '\n')
- z$probs <- c(x$pridg, x$prcdg, x$prper, x$pripr, x$prdpr, x$prirms,
- x$prdrms)
- cat(z$probs,'\n')
- ans <- .Call("mlMakeChains", PACKAGE=pkgname, pData, pModel,
- simpleRates, z$probs, z$prmin, z$prmib,
- x$minimumPermutationLength,
- x$maximumPermutationLength,
- x$initialPermutationLength)
-
- f$minimalChain <- ans[[1]]
- f$chain <- ans[[2]]
- # print(nrow(ans[[1]][[1]]))
- # print(nrow(ans[[2]][[1]]))
- # browser()
- ##store address of simulation object
- # f$pMLSimulation <- ans[[1]][[1]]
- #ans1 <- reg.finalizer(f$pMLSimulation, clearMLSimulation,
- # onexit = FALSE)
+ ans <- .Call("mlInitializeSubProcesses",
+ PACKAGE=pkgname, pData, pModel,
+ simpleRates, z$probs, z$prmin, z$prmib,
+ x$minimumPermutationLength,
+ x$maximumPermutationLength,
+ x$initialPermutationLength, ff$chain, ff$missingChain)
+ f$chain <- ff$chain
+ f$missingChain <- ff$missingChain
+ }
}
f$myeffects <- myeffects
f$myCompleteEffects <- myCompleteEffects
@@ -551,7 +574,9 @@
z$f <- f
z <- initForAlgorithms(z)
z$periodNos <- attr(data, "periodNos")
- z$f[1:nGroup] <- NULL
+ if (! returnDeps) {
+ z$f[1:nGroup] <- NULL
+ }
}
if (initC || (z$int == 1 && z$int2 == 1 &&
(is.null(z$nbrNodes) || z$nbrNodes == 1)))
@@ -606,15 +631,15 @@
## add attribute of size
if (bipartite)
{
- attr(mat1,'nActors') <- c(nrow(mat), ncol(mat))
- attr(mat2,'nActors') <- c(nrow(mat), ncol(mat))
- attr(mat3,'nActors') <- c(nrow(mat), ncol(mat))
+ attr(mat1, 'nActors') <- c(nrow(mat), ncol(mat))
+ attr(mat2, 'nActors') <- c(nrow(mat), ncol(mat))
+ attr(mat3, 'nActors') <- c(nrow(mat), ncol(mat))
}
else
{
- attr(mat1,'nActors') <- nrow(mat)
- attr(mat2,'nActors') <- nrow(mat)
- attr(mat3,'nActors') <- nrow(mat)
+ attr(mat1, 'nActors') <- nrow(mat)
+ attr(mat2, 'nActors') <- nrow(mat)
+ attr(mat3, 'nActors') <- nrow(mat)
}
list(mat1 = t(mat1), mat2 = t(mat2), mat3 = t(mat3))
@@ -642,8 +667,8 @@
}, y = matorig)
mat2 <- do.call(rbind, tmp)
## add attribute of size
- attr(mat1,'nActors1') <- nrow(mat)
- attr(mat1,'nActors2') <- ncol(mat)
+ attr(mat1, 'nActors1') <- nrow(mat)
+ attr(mat1, 'nActors2') <- ncol(mat)
list(mat1=t(mat1), mat2=t(mat2))
}
##@unpackOneMode siena07 Reformat data for C++
@@ -985,14 +1010,14 @@
if (i < observations)
{
## recreate distances, as we have none in c++. (no longer true)
- mymat1 <- depvar[,,i, drop=FALSE]
- mymat2 <- depvar[,,i + 1,drop=FALSE]
+ mymat1 <- depvar[, ,i, drop=FALSE]
+ mymat2 <- depvar[, ,i + 1, drop=FALSE]
##remove structural values
- mymat1[mymat1 %in% c(10,11)] <- NA
- mymat2[mymat2 %in% c(10,11)] <- NA
+ mymat1[mymat1 %in% c(10, 11)] <- NA
+ mymat2[mymat2 %in% c(10, 11)] <- NA
## and the diagonal
- diag(mymat1[,,1]) <- 0
- diag(mymat2[,,1]) <- 0
+ diag(mymat1[, ,1]) <- 0
+ diag(mymat2[, ,1]) <- 0
mydiff <- mymat2 - mymat1
attr(depvar, 'distance')[i] <- sum(mydiff != 0,
na.rm = TRUE)
@@ -1606,6 +1631,11 @@
gsub("#", y$parm, y$functionName)
}, y=effects)
+ if (any(effects$shortName == "behUnspInt" & effects$include &
+ effects$effect1 > 0))
+ {
+ stop("User specified behavior interactions are not yet implemented")
+ }
##validate user-specified network interactions
interactions <- effects[effects$shortName == "unspInt" & effects$include &
effects$effect1 > 0, ]
@@ -1651,8 +1681,8 @@
# "with different specifications eval/endow/rate ",
# "trying with experimental code. Remove these ",
# "Interactions if this does not work.")
- stop("invalid interaction specification: ",
- "must be same type: evaluation or endowment")
+ stop("invalid interaction specification: must",
+ "be same type: evaluation, endowment or creation")
}
}
else
@@ -1666,8 +1696,8 @@
if (inter1$type != inter2$type ||
inter1$type != inter3$type)
{
- stop("invalid interaction specification:",
- "must all be same type: evaluation or endowment")
+ stop("invalid interaction specification: must all be",
+ "same type: evaluation, endowment or creation")
}
}
## check types
Modified: pkg/RSiena/R/maxlikec.r
===================================================================
--- pkg/RSiena/R/maxlikec.r 2011-05-29 18:42:29 UTC (rev 152)
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/rsiena -r 153
More information about the Rsiena-commits
mailing list