[Rsiena-commits] r150 - in pkg: RSiena RSiena/R RSiena/man RSiena/src RSiena/src/model/ml RSienaTest RSienaTest/R RSienaTest/doc RSienaTest/inst/doc RSienaTest/man RSienaTest/src RSienaTest/src/model/ml
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri May 27 17:20:18 CEST 2011
Author: ripleyrm
Date: 2011-05-27 17:20:17 +0200 (Fri, 27 May 2011)
New Revision: 150
Modified:
pkg/RSiena/DESCRIPTION
pkg/RSiena/R/RSienaRDocumentation.r
pkg/RSiena/R/Sienatest.r
pkg/RSiena/R/bayes.r
pkg/RSiena/R/effects.r
pkg/RSiena/R/initializeFRAN.r
pkg/RSiena/R/maxlike.r
pkg/RSiena/R/maxlikecalc.r
pkg/RSiena/R/observationErrors.r
pkg/RSiena/R/phase1.r
pkg/RSiena/R/phase2.r
pkg/RSiena/R/phase3.r
pkg/RSiena/R/print01Report.r
pkg/RSiena/R/printDataReport.r
pkg/RSiena/R/printInitialDescription.r
pkg/RSiena/R/robmon.r
pkg/RSiena/R/siena01.r
pkg/RSiena/R/siena07.r
pkg/RSiena/R/siena08.r
pkg/RSiena/R/sienaDataCreate.r
pkg/RSiena/R/sienaDataCreateFromSession.r
pkg/RSiena/R/simstatsc.r
pkg/RSiena/changeLog
pkg/RSiena/cleanup
pkg/RSiena/cleanup.win
pkg/RSiena/man/RSiena-package.Rd
pkg/RSiena/src/Makevars
pkg/RSiena/src/Makevars.win
pkg/RSiena/src/model/ml/Chain.cpp
pkg/RSiena/src/siena07models.cpp
pkg/RSiena/src/siena07utilities.cpp
pkg/RSienaTest/DESCRIPTION
pkg/RSienaTest/R/RSienaRDocumentation.r
pkg/RSienaTest/R/Sienatest.r
pkg/RSienaTest/R/bayes.r
pkg/RSienaTest/R/initializeFRAN.r
pkg/RSienaTest/R/maxlike.r
pkg/RSienaTest/R/maxlikecalc.r
pkg/RSienaTest/R/observationErrors.r
pkg/RSienaTest/R/phase1.r
pkg/RSienaTest/R/phase2.r
pkg/RSienaTest/R/phase3.r
pkg/RSienaTest/R/print01Report.r
pkg/RSienaTest/R/printDataReport.r
pkg/RSienaTest/R/printInitialDescription.r
pkg/RSienaTest/R/robmon.r
pkg/RSienaTest/R/siena01.r
pkg/RSienaTest/R/siena07.r
pkg/RSienaTest/R/siena08.r
pkg/RSienaTest/R/sienaDataCreate.r
pkg/RSienaTest/R/sienaDataCreateFromSession.r
pkg/RSienaTest/R/sienaGOF.r
pkg/RSienaTest/R/simstatsc.r
pkg/RSienaTest/changeLog
pkg/RSienaTest/cleanup
pkg/RSienaTest/cleanup.win
pkg/RSienaTest/doc/RSiena.bib
pkg/RSienaTest/doc/s_man400.tex
pkg/RSienaTest/inst/doc/s_man400.pdf
pkg/RSienaTest/man/RSiena-package.Rd
pkg/RSienaTest/src/Makevars
pkg/RSienaTest/src/Makevars.win
pkg/RSienaTest/src/model/ml/Chain.cpp
pkg/RSienaTest/src/siena07models.cpp
pkg/RSienaTest/src/siena07utilities.cpp
Log:
Removing check warnings, fixing manual and bibtex errors
Modified: pkg/RSiena/DESCRIPTION
===================================================================
--- pkg/RSiena/DESCRIPTION 2011-05-26 16:14:47 UTC (rev 149)
+++ pkg/RSiena/DESCRIPTION 2011-05-27 15:20:17 UTC (rev 150)
@@ -1,8 +1,8 @@
Package: RSiena
Type: Package
Title: Siena - Simulation Investigation for Empirical Network Analysis
-Version: 1.0.12.148
-Date: 2011-05-25
+Version: 1.0.12.150
+Date: 2011-05-27
Author: Various
Depends: R (>= 2.10.0), xtable
Imports: Matrix
Modified: pkg/RSiena/R/RSienaRDocumentation.r
===================================================================
--- pkg/RSiena/R/RSienaRDocumentation.r 2011-05-26 16:14:47 UTC (rev 149)
+++ pkg/RSiena/R/RSienaRDocumentation.r 2011-05-27 15:20:17 UTC (rev 150)
@@ -58,7 +58,7 @@
mystr <- paste("##", "@", sep="")
comms1 <- strsplit(comms, mystr)
## join up rest
- comms2 <- do.call(rbind, comms1)
+ ## comms2 <- do.call(rbind, comms1)
## turn into dataframe
comms3 <- sapply(comms1, function(x)
{
@@ -107,7 +107,6 @@
tt <- lapply(1:length(ttmp), function(x, y, z)
{
yy <- y[x]
- zz <- z[[x]]
yy <- getFromNamespace(yy, pkgname)
targs <- formals(yy)
n <- length(targs)
@@ -126,7 +125,7 @@
yy <- y[[x]]
n <- length(y[[x]])
bb <- names(yy)
- t1<- lapply(1:length(yy), function(x, b, a)
+ t1<- lapply(1:n, function(x, b, a)
{
y <- a[[x]]
bb <- b[[x]]
@@ -136,7 +135,7 @@
else
c( bb, " ")
}, a=yy, b=bb)
- t2 <- do.call(rbind,t1)
+ do.call(rbind,t1)
}, y=tt
)
@@ -200,7 +199,7 @@
names(tttt7) <- c("Called from", "Function")
tttt7 <- tttt7[order(tttt7[,2],tttt7[,1]), ]
- tmp7bit <- tmp7[tmp7$Function %in% tttt7$Function, ]
+ ## tmp7bit <- tmp7[tmp7$Function %in% tttt7$Function, ]
tmp7new <- merge(tmp7, tttt7, by=c("Function", "Called from"), all=TRUE)
Modified: pkg/RSiena/R/Sienatest.r
===================================================================
--- pkg/RSiena/R/Sienatest.r 2011-05-26 16:14:47 UTC (rev 149)
+++ pkg/RSiena/R/Sienatest.r 2011-05-27 15:20:17 UTC (rev 150)
@@ -16,7 +16,7 @@
## 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)
+ # pp<- length(z$diver)
constant<- z$diver
test<- z$test
covtheta<- z$covtheta
Modified: pkg/RSiena/R/bayes.r
===================================================================
--- pkg/RSiena/R/bayes.r 2011-05-26 16:14:47 UTC (rev 149)
+++ pkg/RSiena/R/bayes.r 2011-05-27 15:20:17 UTC (rev 150)
@@ -100,7 +100,7 @@
{
u <- 1 / (2 - (actual / desired))
}
- if (tol <- abs(actual - desired) <= tolerance)
+ if (abs(actual - desired) <= tolerance)
{
number <<- number + 1
if (number == 2) success <<- TRUE
@@ -172,7 +172,7 @@
npar <- length(z$theta)
basicRate <- z$effects$basicRate
- iter <- 0
+ ##iter <- 0
z$numm <- 20
z$scaleFactor <- 1
z <- createStores()
Modified: pkg/RSiena/R/effects.r
===================================================================
--- pkg/RSiena/R/effects.r 2011-05-26 16:14:47 UTC (rev 149)
+++ pkg/RSiena/R/effects.r 2011-05-27 15:20:17 UTC (rev 150)
@@ -264,26 +264,20 @@
{
objEffects <-
rbind(objEffects,
- createEffects("covarNetNetObjective",
- otherName,
- names(xx$cCovars)[k],
- name=varname,
- groupName=groupName, group=group,
- netType=netType))
+ covarNetNetEff(otherName, names(xx$cCovars)[k],
+ attr(xx$cCovars[[k]], 'poszvar'),
+ name=varname))
}
}
- for (k in seq(along=xx$vCovars))
+ for (k in seq(along=xx$vCovars))
{
if (attr(xx$vCovars[[k]], 'nodeSet') == nodeSet)
{
objEffects <-
rbind(objEffects,
- createEffects("covarNetNetObjective",
- otherName,
- names(xx$vCovars)[k],
- name=varname,
- groupName=groupName, group=group,
- netType=netType))
+ covarNetNetEff(otherName, names(xx$vCovars)[k],
+ attr(xx$vCovars[[k]], 'poszvar'),
+ name=varname))
}
}
for (k in seq(along=xx$depvars))
@@ -293,12 +287,9 @@
{
objEffects <-
rbind(objEffects,
- createEffects("covarNetNetObjective",
- otherName,
- names(xx$depvars)[k],
- name=varname,
- groupName=groupName, group=group,
- netType=netType))
+ covarNetNetEff(otherName, names(xx$depvars)[k],
+ poszvar=TRUE,
+ name=varname))
}
}
@@ -862,13 +853,13 @@
}
##@covarNetNetEff internal getEffects
covarNetNetEff<- function(othernetname,
- covarname, poszvar, moreThan2, name)
+ covarname, poszvar, name)
{
objEffects <- createEffects("covarNetNetObjective", othernetname,
covarname, name=name,
groupName=groupName, group=group,
netType=netType)
- if (!(poszvar && (!moreThan2)))
+ if (!poszvar)
{
objEffects <- objEffects[objEffects$shortName != "covNetNet", ]
}
Modified: pkg/RSiena/R/initializeFRAN.r
===================================================================
--- pkg/RSiena/R/initializeFRAN.r 2011-05-26 16:14:47 UTC (rev 149)
+++ pkg/RSiena/R/initializeFRAN.r 2011-05-27 15:20:17 UTC (rev 150)
@@ -119,7 +119,6 @@
interactionNos <- interactionNos[interactionNos > 0]
interactions <- effects$effectNumber %in%
interactionNos
- interactionMainEffects <- effects[interactions, ]
effects$requested <- effects$include
requestedEffects <- effects[effects$include, ]
@@ -168,7 +167,7 @@
}
}
types <- sapply(data[[1]]$depvars, function(x) attr(x, 'type'))
- nets <- sum(types != "behavior")
+ ##nets <- sum(types != "behavior")
## now check if conditional estimation is OK and copy to z if so
z$cconditional <- FALSE
if (x$cconditional)
@@ -1354,7 +1353,6 @@
{
beh <- depvar[, 1, ]
behmiss <- is.na(beh)
- origbeh <- beh
allna <- apply(beh, 1, function(x)all(is.na(x)))
modes <- attr(depvar, "modes")
## carry forward missings ### nb otherwise use the mode
@@ -1463,7 +1461,6 @@
unpackVDyad<- function(dyvCovar, observations)
{
edgeLists <- vector('list', observations)
- varmats <- vector('list', observations)
sparse <- attr(dyvCovar, 'sparse')
means <- attr(dyvCovar, "meanp")
nodeSets <- attr(dyvCovar, "nodeSet")
@@ -1565,7 +1562,6 @@
atts <- attributes(compositionChange)
events <- atts$events
activeStart <- atts$activeStart
- nActors <- nrow(activeStart)
observations <- ncol(activeStart)
## check that there is someone there always
for (i in 1:(observations - 1))
Modified: pkg/RSiena/R/maxlike.r
===================================================================
--- pkg/RSiena/R/maxlike.r 2011-05-26 16:14:47 UTC (rev 149)
+++ pkg/RSiena/R/maxlike.r 2011-05-27 15:20:17 UTC (rev 150)
@@ -1,33 +1,40 @@
-maxlikefn<- function(z,x,INIT=FALSE,TERM=FALSE, data, effects=NULL,nstart=1000,
- pinsdel=0.6,pperm=0.3,prelins=0.1,multfactor=2.0,
- promul=0.1,promul0=0.5,pdiaginsdel=0.1,
+maxlikefn<- function(z, x, INIT=FALSE, TERM=FALSE, data, effects=NULL,
+ nstart=1000,
+ pinsdel=0.6, pperm=0.3, prelins=0.1, multfactor=2.0,
+ promul=0.1, promul0=0.5, pdiaginsdel=0.1,
fromFiniteDiff=FALSE, noSamples=1, sampInterval=50, int=1)
{
- mlInit<- function(z,x,data,effects)
+ mlInit <- function(z, x, data, effects)
{
f <- NULL
- if (!inherits(data, 'siena'))
- stop('not valid siena data object')
+ if (!inherits(data, "siena"))
+ {
+ stop("not valid siena data object")
+ }
if (is.null(effects))
+ {
effects <- getEffects(data)
+ }
if (!is.data.frame(effects))
+ {
stop('effects is not a data.frame')
- effects <- effects[effects$include,]
+ }
+ effects <- effects[effects$include, ]
z$theta <- effects$initialValue
z$fixed <- effects$fix
z$test <- effects$test
z$pp <- length(z$test)
- z$posj <- rep(FALSE,z$pp)
+ z$posj <- rep(FALSE, z$pp)
z$targets <- rep(0, z$pp)
## effectsNames<- getEffectNames(effects)
- z$posj[grep('basic', effects$effectName)] <- TRUE
- z$posj[grep('constant', effects$effectName)] <- TRUE
+ z$posj[grep("basic", effects$effectName)] <- TRUE
+ z$posj[grep("constant", effects$effectName)] <- TRUE
z$BasicRateFunction <- z$posj
observations <- data$observations
- mats <- vector('list', observations)
- f$mynets <- vector('list', observations)
- types <- sapply(data$depvars, function(x)attr(x,'type'))
- netsubs <- which(types=='oneMode')
+ mats <- vector("list", observations)
+ f$mynets <- vector("list", observations)
+ types <- sapply(data$depvars, function(x) attr(x, "type"))
+ netsubs <- which(types=="oneMode")
netsub <- min(netsubs) ### only one for now
actsubs <- which(types=='behavior')
for (i in 1:observations)
@@ -35,10 +42,14 @@
mats[[i]] <- data$depvars[[netsub]][, , i]
f$mynets[[i]] <- mats[[i]]
if (i==1)
- f$mynets[[i]][is.na(mats[[i]])] <-0
+ {
+ f$mynets[[i]][is.na(mats[[i]])] <- 0
+ }
else ##carry missing forward!
+ {
f$mynets[[i]][is.na(mats[[i]])] <-
f$mynets[[i - 1]][is.na(mats[[i]])]
+ }
f$mynets[[i]][mats[[i]]==10] <- 0
f$mynets[[i]][mats[[i]]==11] <- 1
}
@@ -56,28 +67,28 @@
f$mats[[i]][mats[[i]]==11] <- 1
f$mats[[i]][mats[[i]]==10] <- 0
}
- if (length(actsubs)>0)
+ if (length(actsubs) > 0)
{
acts <- matrix(data$depvars[[actsubs[1]]],
ncol=observations)
f$acts <- acts
f$myacts <- acts
f$myacts[is.na(acts)] <- 0
- f$meanact <- round(mean(acts,na.rm=TRUE))
+ f$meanact <- round(mean(acts, na.rm=TRUE))
}
f$observations <- observations
## browser()
if (any(z$targets!=0))
{
- Report(c('Targets should be zero for maximum likelihood:',
- 'they have been zeroed\n'))
+ Report(c("Targets should be zero for maximum likelihood:',
+ 'they have been zeroed\n"))
z$targets <- rep(0, z$pp)
}
mat1 <- data$depvars[[netsub]][, , 1]
mat2 <- data$depvars[[netsub]][, , 2]
- # f$mat1<- mat1
- # f$mat2<- mat2
+ ## f$mat1<- mat1
+ ## f$mat2<- mat2
startmat <- mat1
startmat[is.na(startmat)] <- 0
endmat <- mat2
@@ -146,7 +157,7 @@
{
f <- FRANstore()
niter <- f$niter
- nactors <- nrow(f$startmat)
+ ## nactors <- nrow(f$startmat)
promul <- promul0
# int <- x$int
if (z$Phase==2)
@@ -262,9 +273,9 @@
tmpnet[ii,jj] <- 1- tmpnet[ii,jj]
}
diag(tmpnet)<- 0
- ps<- calcprobs(i,tmpnet,betapar,nactors)
- pr<- 1/sum(ps)
- ndiag<-nrow(chain[chain[,1]==chain[,2],,drop=FALSE])
+ ps <- calcprobs(i,tmpnet,betapar,nactors)
+ pr <- 1/sum(ps)
+ ndiag <- nrow(chain[chain[, 1] == chain[, 2], , drop=FALSE])
## ans <- kappasigmamu(nactors,nrow(chain),lambda,add1=TRUE)
## mu<- ans$mu + 1/lambda/nactors
## sigma2 <- ans$sigma2 + 1/lambda/lambda/nactors/nactors
@@ -328,51 +339,65 @@
list(ccps=ccps, nc=nc)
}
##fix up numm
- numm<- f$numm
+ numm <- f$numm
if (numm > nrow(chain))
+ {
numm <- nrow(chain)
+ }
if (numm > 40)
+ {
numm <- 40
+ }
if (numm < 2)
+ {
numm <- 2
+ }
num <- trunc(numm)
if (!f$madechain)
+ {
num <- 0
+ }
## else
## cat('num=', num, 'numm=', numm, '\n')
if (insdel)
{
- mults<-as.matrix(unique(chain[duplicated(chain[,1:2])&
- chain[,1]!=chain[,2],,drop=FALSE]))
- nmul<- nrow(mults)
+ mults <- as.matrix(unique(chain[duplicated(chain[, 1:2]) &
+ chain[, 1]!=chain[, 2], , drop=FALSE]))
+ nmul <- nrow(mults)
if (is.null(nmul))
nmul <- 0
- if (nmul>0 && runif(1)<promul)
+ if (nmul > 0 && runif(1) < promul)
{
##choose one of unique duplicates
- mypair <- mults[sample(1:nmul,1),]
+ mypair <- mults[sample(1:nmul, 1), ]
}
else
{
- mypair <- sample(1:nactors,2)
+ mypair <- sample(1:nactors, 2)
}
from <- mypair[1]
to <- mypair[2]
- similar<- chain[chain[,1]==from&chain[,2]==to,,drop=FALSE]
- nk<- nrow(similar)
+ similar <- chain[chain[, 1] == from & chain[, 2]==to, , drop=FALSE]
+ nk <- nrow(similar)
##nk is number of this connection
tmp <- findCCPs(similar)
nc <- tmp$nc
ccps <- tmp$ccps
## nc<- nrow(similar[similar[,3]>0,,drop=FALSE])/2
if (nc < 1)
- ins<- TRUE
+ {
+ ins <- TRUE
+ }
else
{
if (runif(1) < prelins)
+ {
ins <- TRUE
+ }
else
+ {
ins <- FALSE
+ }
}
del <- !ins
if (ins)
@@ -389,39 +414,43 @@
else
{
## if (nc>0)
-#### {
- ## ccps<- chain[chain[,1]==from&chain[,2]==to&
- ## chain[,3]>0,]
- ## ccpsubs<- unique(ccps[,3])
- ## subslist<- 1:(max(ccpsubs)+1)
- ## newsub<- min(subslist[!subslist %in% ccpsubs])
+ ## {
+ ## ccps <- chain[chain[, 1]==from & chain[, 2]==to &
+ ## chain[, 3] > 0, ]
+ ## ccpsubs <- unique(ccps[, 3])
+ ## subslist <- 1:(max(ccpsubs) + 1)
+ ## newsub <- min(subslist[!subslist %in% ccpsubs])
## }
## else
- ## newsub<- 1
- samplist<- 1:nrow(chain)
- samplist<- samplist[!(from==chain[,1]&to==chain[,2])]
- samplist<- c(samplist,nrow(chain)+1)
- numav1<- length(samplist)
+ ## {
+ ## newsub <- 1
+ ## }
+ samplist <- 1:nrow(chain)
+ samplist <- samplist[!(from==chain[, 1] & to==chain[, 2])]
+ samplist <- c(samplist, nrow(chain) + 1)
+ numav1 <- length(samplist)
- k1<- sample(samplist,1)
+ k1 <- sample(samplist, 1)
## find last previous and next of this kind
- thiskind<-(1:nrow(chain))[from==chain[,1]&to==chain[,2]]
- k<- max(c(0,thiskind[thiskind<k1]))+1
- kk<- min(c(thiskind[thiskind>k1],nrow(chain)+1))
+ thiskind <- (1:nrow(chain))[from==chain[, 1] & to==chain[, 2]]
+ k <- max(c(0, thiskind[thiskind<k1])) + 1
+ kk <- min(c(thiskind[thiskind>k1], nrow(chain) + 1))
if (k>nrow(chain)) ##not possible to proceed
- return(list(chain=chain,accept=FALSE))
- numav2<- kk-k
+ {
+ return(list(chain=chain, accept=FALSE))
+ }
+ numav2 <- kk - k
- k2<- sample(1:numav2,1)
+ k2 <- sample(1:numav2, 1)
if (k2==k1)
{
k2 <- kk
}
- if (k2<k1)
+ if (k2 < k1)
{
- kk<- k2
- k2<- k1
- k1<- kk
+ kk <- k2
+ k2 <- k1
+ k1 <- kk
}
}
}
@@ -649,34 +678,52 @@
##cat('probrat ',probrat1,' ')
## probrat <- prdelphi * pp
## cat(probrat,' ',probrat-probrat1,'\n')
- if (insdel&ins)
- nn<- 3
- else if (insdel&del)
- nn<- 4
+ ## if (insdel & ins)
+ ## {
+ ## nn <- 3
+ ## }
+ ##else if (insdel & del)
+ ##{
+ ## nn <- 4
+ ##}
+ ##else
+ ##{
+ ## nn <- 5
+ ##}
+ ## cat(probrat, '\n')
+ accept <- FALSE
+ if (probrat > 1)
+ {
+ accept <- TRUE
+ }
else
- nn<- 5
- ## cat(probrat, '\n')
- accept<- FALSE
- if (probrat>1)
{
- accept<- TRUE
+ if (runif(1) < probrat)
+ {
+ accept <- TRUE
+ }
}
- else
- if (runif(1)< probrat)
- accept<- TRUE
- if (sum(tempchain[,3]==1)%%2!=0)
+ if (sum(tempchain[, 3]==1) %% 2!=0)
+ {
browser()
+ }
if (accept)
- chain<- tempchain
+ {
+ chain <- tempchain
+ }
## cat(num,f$numm,accept , insdel,truncated,accept,'\n')
## if (!insdel)
## browser()
if (accept && !insdel && !truncated)
- numm<- numm+0.5
+ {
+ numm <- numm + 0.5
+ }
if (!accept && !insdel && !truncated)
- numm <- numm-0.5
+ {
+ numm <- numm - 0.5
+ }
#browser()
- list(chain=chain,accept=accept,numm=numm)
+ list(chain=chain, accept=accept, numm=numm)
}##end of procedure
#########################################################################
##start of mhstep
Modified: pkg/RSiena/R/maxlikecalc.r
===================================================================
--- pkg/RSiena/R/maxlikecalc.r 2011-05-26 16:14:47 UTC (rev 149)
+++ pkg/RSiena/R/maxlikecalc.r 2011-05-27 15:20:17 UTC (rev 150)
@@ -127,21 +127,25 @@
grad <- derivs(pars, Z, varmat)
grad
}
-kappa1<- function(n,nc,lambda)
+kappa1 <- function(n, nc, lambda)
{
- nc<- nc+1
- mu<- nc/n/lambda
- sigma2<- nc/n/n/lambda/lambda
- logp<- -(1-mu)^2/2/sigma2 - log(sigma2)/2 -log(2*pi)/2 -log(n*lambda)
- pp<-dpois(nc,n*lambda,log=TRUE)
- # cat(logp,pp,'\n')
- -logp+nc*log(n)
+ nc <- nc + 1
+ mu <- nc/ n / lambda
+ sigma2 <- nc/ n / n / lambda / lambda
+ logp <- -(1 - mu)^2 / 2 / sigma2 - log(sigma2) / 2 - log(2 * pi) / 2 -
+ log(n * lambda)
+ ## pp <- dpois(nc, n * lambda, log=TRUE)
+ ## cat(logp, pp, '\n')
+ -logp + nc * log(n)
}
-kappasigmamu<- function(n,nc,lambda,add1=FALSE)
+kappasigmamu<- function(n, nc, lambda, add1=FALSE)
{
- if (add1) nc<- nc+1
- mu<- nc/n/lambda
- sigma2<- nc/n/n/lambda/lambda
- kappa1<- -(1-mu)^2/2/sigma2 - log(sigma2)/2
- list(kappa= kappa1,mu=mu,sigma2=sigma2)
+ if (add1)
+ {
+ nc <- nc + 1
+ }
+ mu <- nc / n / lambda
+ sigma2 <- nc / n / n / lambda / lambda
+ kappa1 <- -(1 - mu)^2 / 2 / sigma2 - log(sigma2) / 2
+ list(kappa=kappa1, mu=mu, sigma2=sigma2)
}
Modified: pkg/RSiena/R/observationErrors.r
===================================================================
--- pkg/RSiena/R/observationErrors.r 2011-05-26 16:14:47 UTC (rev 149)
+++ pkg/RSiena/R/observationErrors.r 2011-05-27 15:20:17 UTC (rev 150)
@@ -17,9 +17,9 @@
##@deviance.observationErrors siena08: deviance under normality assumptions
deviance.observationErrors <- function(mu, sig, x, se)
{
-## deviance.observationErrors: given vector observations x, assuming model
-## x ~ normal(expected = mu, variance = sig^2 + se^2)
-## required is length(x) = length(se)
+ ## deviance.observationErrors: given vector observations x, assuming model
+ ## x ~ normal(expected = mu, variance = sig^2 + se^2)
+ ## required is length(x) = length(se)
sig2 <- sig^2
sum((x - mu)^2/(sig2 + se^2)) + sum(log(sig2 + se^2))
}
@@ -27,36 +27,46 @@
##@profdev.mu siena08 (profile deviance =) minus twice profile loglik for mu
profdev.mu <- function(mu,x,se)
{
-## deviance.observationErrors maximized over sigma
+ ## deviance.observationErrors maximized over sigma
ra <- max(x) - min(x) # easy upper bound for sigma
optimize(function(sig)deviance.observationErrors(mu, sig, x, se),
- c(0, ra))
+ c(0, ra))
}
##@profdev.sig siena08 (profile deviance =) minus twice prof. loglik for sigma
profdev.sig <- function(sig,x,se)
{
-## deviance.observationErrors maximized over mu
+ ## deviance.observationErrors maximized over mu
optimize(function(mu)deviance.observationErrors(mu, sig, x, se),
- range(x))
+ range(x))
}
##@maxlik siena08 maximum likelihood estimator
maxlik <- function(x,se)
{
- ## deviance.observationErrors minimized over mu and sigma
- ## Minimize profile deviance for mu:
- opmu <- optimize(function(mu)profdev.mu(mu, x, se)$objective,
- range(x))
- ## MLE for mu:
- mu <- opmu$minimum
- ## minimized deviance:
- dev <- opmu$objective
- ## Location for minimum of deviance for this value of mu:
- sig <- profdev.mu(mu, x ,se)$minimum
- ## Standard error of MLE(mu):
- se.mu <- sqrt(1 / sum(1 / (se^2 + sig^2)))
+ if (length(x) > 1)
+ {
+ ## deviance.observationErrors minimized over mu and sigma
+ ## Minimize profile deviance for mu:
+ opmu <- optimize(function(mu)profdev.mu(mu, x, se)$objective,
+ range(x))
+ ## MLE for mu:
+ mu <- opmu$minimum
+ ## minimized deviance:
+ dev <- opmu$objective
+ ## Location for minimum of deviance for this value of mu:
+ sig <- profdev.mu(mu, x , se)$minimum
+ ## Standard error of MLE(mu):
+ se.mu <- sqrt(1 / sum(1 / (se^2 + sig^2)))
+ }
+ else
+ {
+ mu <- x
+ se.mu <- NA
+ sig <- NA
+ dev <- NA
+ }
return(list(mu = mu, se.mu = se.mu, sigma = sig, deviance = dev))
}
@@ -88,7 +98,7 @@
{
x2 <- x2 + ra
}
- if (f(x1)*f(x2) < 0)
+ if (f(x1) * f(x2) < 0)
{
break
}
@@ -104,7 +114,7 @@
}
}
-##@ confint.mu siena08 confidence interval for mu
+##@confint.mu siena08 confidence interval for mu
confint.mu <- function(x, se, alpha=0.05)
{
## confidence interval for mu in model
@@ -114,33 +124,40 @@
## confidence level (1-alpha).
ma <- max(x)
mi <- min(x)
- ra <- max(x) - min(x) # easy upper bound for sigma
maxli <- maxlik(x, se) # ML estimators
mindev <- maxli$deviance # minimized deviance
mlemu <- maxli$mu # MLE for mu
chidev <- qchisq(1 - alpha, 1) # critical value chi-squared distribution
- ## The end points of the confidence interval are the values
- ## where the profile deviance for mu is equal to mindev + chidev.
- ## Now first the solution for the left side of the confidence interval:
- tmp1 <- unisroot(function(mu)
- {
- tmp <- profdev.mu(mu, x, se)
- tmp$objective - mindev - chidev
- },
- c(mi, mlemu))
- mu1 <- tmp1$root
- ## and then the solution for the right side of the confidence interval:
- tmp1 <- unisroot(function(mu)
- {
- tmp <- profdev.mu(mu, x, se)
- tmp$objective - mindev - chidev
- },
- c(mlemu, ma), left=FALSE)
- mu2 <- tmp1$root
+ if (length(x) > 1)
+ {
+ ## The end points of the confidence interval are the values
+ ## where the profile deviance for mu is equal to mindev + chidev.
+ ## Now first the solution for the left side of the confidence interval:
+ tmp1 <- unisroot(function(mu)
+ {
+ tmp <- profdev.mu(mu, x, se)
+ tmp$objective - mindev - chidev
+ },
+ c(mi, mlemu))
+ mu1 <- tmp1$root
+ ## and then the solution for the right side of the confidence interval:
+ tmp1 <- unisroot(function(mu)
+ {
+ tmp <- profdev.mu(mu, x, se)
+ tmp$objective - mindev - chidev
+ },
+ c(mlemu, ma), left=FALSE)
+ mu2 <- tmp1$root
+ }
+ else
+ {
+ mu1 <- NA
+ mu2 <- NA
+ }
c(mu1, mu2, 1 - alpha)
}
-##@ confint.sig siena08 confidence interval for sigma
+##@confint.sig siena08 confidence interval for sigma
confint.sig <- function(x, se, alpha=0.05)
{
## confidence interval for sigma in model
@@ -148,7 +165,6 @@
## with length(x) = length(se)
## returns list consisting of bounds confidence interval and
## confidence level (1 - alpha). ma <- max(x)
- mi <- min(x)
ra <- max(x) - min(x) # easy upper bound for sigma
maxli <- maxlik(x, se) # ML estimators
mindev <- maxli$deviance # minimized deviance
@@ -159,20 +175,28 @@
## unless the profile deviance in 0 is smaller than this.
## Now first the solution for the left side of the confidence interval;
## this may be 0.
- if (profdev.sig(0, x, se)$objective <= mindev + chidev)
+ if (length(x) > 1)
{
- sig1 <- 0
+ if (profdev.sig(0, x, se)$objective <= mindev + chidev)
+ {
+ sig1 <- 0
+ }
+ else
+ {
+ sig1 <- uniroot(function(sig)
+ profdev.sig(sig, x, se)$objective - mindev - chidev,
+ c(0, mlesig))$root
+ }
+ ## and then the solution for the right side of the confidence interval:
+ sig2 <- unisroot(function(sig)
+ profdev.sig(sig, x, se)$objective - mindev - chidev,
+ c(mlesig, ra), left=FALSE)$root
}
else
{
- sig1 <- uniroot(function(sig)
- profdev.sig(sig, x, se)$objective - mindev - chidev,
- c(0, mlesig))$root
+ sig1 <- NA
+ sig2 <- NA
}
- ## and then the solution for the right side of the confidence interval:
- sig2 <- unisroot(function(sig)
- profdev.sig(sig,x,se)$objective - mindev - chidev,
- c(mlesig, ra), left=FALSE)$root
return(c(sig1, sig2, 1 - alpha))
}
Modified: pkg/RSiena/R/phase1.r
===================================================================
--- pkg/RSiena/R/phase1.r 2011-05-26 16:14:47 UTC (rev 149)
+++ pkg/RSiena/R/phase1.r 2011-05-27 15:20:17 UTC (rev 150)
@@ -408,7 +408,6 @@
##@CalculateDerivative siena07 Calculates derivative in Phase 1
CalculateDerivative <- function(z, x)
{
- f <- FRANstore()
if (z$FinDiff.method || x$maxlike)
{
dfra <- t(apply(z$sdf, c(2, 3), mean))
Modified: pkg/RSiena/R/phase2.r
===================================================================
--- pkg/RSiena/R/phase2.r 2011-05-26 16:14:47 UTC (rev 149)
+++ pkg/RSiena/R/phase2.r 2011-05-27 15:20:17 UTC (rev 150)
@@ -31,8 +31,8 @@
z$phase2fras <- array(0, dim=c(4, z$pp, 1000))
z$rejectprops <- array(0, dim=c(4, 4, 1000))
}
- int <- 1
- f <- FRANstore()
+ ## int <- 1
+ ## f <- FRANstore()
z$Phase <- 2
z$writefreq <- 1
if (!is.batch())
@@ -107,7 +107,6 @@
## report truncations and positivizations
if (any(z$truncated))
{
- truncations<- (1 : length(z$truncated))[z$truncated]
msg<- paste('Intervention 2.', subphase,
'.1: changes truncated, iterations: ',sep='')
Report(msg, cf)
@@ -600,9 +599,9 @@
nc[names(ncvals)] <- ncvals
logChoiceProb <- sapply(chain, function(x)x[[9]])
logOptionSetProb <- sapply(chain, function(x)x[[8]])
- loglik <- sum(logChoiceProb) # + sum(logOptionSetProb)
+ loglik <- sum(logChoiceProb) + sum(logOptionSetProb)
#print(sum(logOptionSetProb))
- loglik <- loglik - sum(nactors * lambda) + sum(nc * log(lambda))
+ #loglik <- loglik - sum(nactors * lambda) + sum(nc * log(lambda))
loglik
}
Modified: pkg/RSiena/R/phase3.r
===================================================================
--- pkg/RSiena/R/phase3.r 2011-05-26 16:14:47 UTC (rev 149)
+++ pkg/RSiena/R/phase3.r 2011-05-27 15:20:17 UTC (rev 150)
@@ -454,7 +454,6 @@
##@CalulateDerivative3 siena07 Calculates derivative at end of phase 3
CalculateDerivative3<- function(z,x)
{
- f <- FRANstore()
z$mnfra <- colMeans(z$sf)
if (z$FinDiff.method || x$maxlike)
{
@@ -479,7 +478,7 @@
z$msf <- cov(z$sf)
if (z$Phase3nits > 2)
{
- z$sfl <- apply(z$sf, 2, function(x)acf(x, plot=FALSE, lag=1)[[1]][[2]])
+ z$sfl <- apply(z$sf, 2, function(x)acf(x, plot=FALSE, lag.max=1)[[1]][[2]])
}
z$dfra1 <- z$dfra
z$dfra <- dfra
Modified: pkg/RSiena/R/print01Report.r
===================================================================
--- pkg/RSiena/R/print01Report.r 2011-05-26 16:14:47 UTC (rev 149)
+++ pkg/RSiena/R/print01Report.r 2011-05-27 15:20:17 UTC (rev 150)
@@ -334,7 +334,6 @@
reportBehaviors <- function()
{
Heading(2, outf, "Reading dependent actor variables.")
- anymissings <- FALSE
iBehav <- 0
for (i in 1:length(x$depvars))
{
@@ -343,10 +342,8 @@
depvar <- x$depvars[[i]]
atts <- attributes(depvar)
netname <- atts$name
- type <- atts$type
iBehav <- iBehav + 1
- mymat <- depvar[, 1, ]
mystr <- paste(iBehav, switch(as.character(iBehav),
"1"=, "21"=, "31"= "st",
"2"=, "22"=, "32"= "nd",
@@ -718,7 +715,6 @@
reportCompositionChange <- function()
{
comps <- x$compositionChange
- nComps <- length(comps)
Heading(2, outf, "Reading files with times of composition change.")
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/rsiena -r 150
More information about the Rsiena-commits
mailing list