[Rsiena-commits] r76 - in pkg/RSienaTest: . R data doc inst/doc inst/examples man src src/data src/model src/model/effects src/model/ml src/model/variables src/network src/utils tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Apr 12 18:01:04 CEST 2010
Author: ripleyrm
Date: 2010-04-12 18:01:03 +0200 (Mon, 12 Apr 2010)
New Revision: 76
Added:
pkg/RSienaTest/R/maxlikec.r
Modified:
pkg/RSienaTest/R/globals.r
pkg/RSienaTest/R/iwlsm.R
pkg/RSienaTest/R/maxlike.r
pkg/RSienaTest/R/maxlikecalc.r
pkg/RSienaTest/R/phase1.r
pkg/RSienaTest/R/phase2.r
pkg/RSienaTest/R/phase3.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/sienaModelCreate.r
pkg/RSienaTest/R/sienaprint.r
pkg/RSienaTest/R/simstatsc.r
pkg/RSienaTest/changeLog
pkg/RSienaTest/data/allEffects.csv
pkg/RSienaTest/doc/RSiena.bib
pkg/RSienaTest/doc/s_man400.tex
pkg/RSienaTest/inst/doc/s_man400.pdf
pkg/RSienaTest/inst/examples/EIESATT.dat
pkg/RSienaTest/inst/examples/FR0.DAT
pkg/RSienaTest/inst/examples/FR1.DAT
pkg/RSienaTest/inst/examples/FR2.DAT
pkg/RSienaTest/inst/examples/FR3.DAT
pkg/RSienaTest/inst/examples/FR4.DAT
pkg/RSienaTest/inst/examples/FR5.DAT
pkg/RSienaTest/inst/examples/FR6.DAT
pkg/RSienaTest/inst/examples/VARS.DAT
pkg/RSienaTest/inst/examples/VRND32T0.DAT
pkg/RSienaTest/inst/examples/VRND32T1.DAT
pkg/RSienaTest/inst/examples/VRND32T2.DAT
pkg/RSienaTest/inst/examples/VRND32T5.DAT
pkg/RSienaTest/inst/examples/cbc3.dat
pkg/RSienaTest/inst/examples/klas12b-advice.dat
pkg/RSienaTest/inst/examples/klas12b-alcohol.dat
pkg/RSienaTest/inst/examples/klas12b-delinquency.dat
pkg/RSienaTest/inst/examples/klas12b-demographics.dat
pkg/RSienaTest/inst/examples/klas12b-net-1.dat
pkg/RSienaTest/inst/examples/klas12b-net-2.dat
pkg/RSienaTest/inst/examples/klas12b-net-3.dat
pkg/RSienaTest/inst/examples/klas12b-net-3m.dat
pkg/RSienaTest/inst/examples/klas12b-net-4.dat
pkg/RSienaTest/inst/examples/klas12b-net-4m.dat
pkg/RSienaTest/inst/examples/klas12b-present.dat
pkg/RSienaTest/inst/examples/klas12b-primary.dat
pkg/RSienaTest/inst/examples/klas12b_netbeh.csv
pkg/RSienaTest/inst/examples/klas12b_netbeha.csv
pkg/RSienaTest/man/coDyadCovar.Rd
pkg/RSienaTest/man/includeEffects.Rd
pkg/RSienaTest/man/siena07.Rd
pkg/RSienaTest/man/sienaModelCreate.Rd
pkg/RSienaTest/src/data/ActorSet.cpp
pkg/RSienaTest/src/data/NetworkConstraint.cpp
pkg/RSienaTest/src/data/NetworkConstraint.h
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/effects/AverageReciprocatedAlterEffect.cpp
pkg/RSienaTest/src/model/effects/Effect.cpp
pkg/RSienaTest/src/model/effects/Effect.h
pkg/RSienaTest/src/model/effects/ReciprocalDegreeBehaviorEffect.cpp
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/EffectValueTable.cpp
pkg/RSienaTest/src/model/variables/NetworkVariable.cpp
pkg/RSienaTest/src/model/variables/NetworkVariable.h
pkg/RSienaTest/src/network/OneModeNetwork.cpp
pkg/RSienaTest/src/siena07.cpp
pkg/RSienaTest/src/utils/SqrtTable.cpp
pkg/RSienaTest/tests/parallel.Rout.save
Log:
many changes. bug fixes, compiler problems, maxlike
Modified: pkg/RSienaTest/R/globals.r
===================================================================
--- pkg/RSienaTest/R/globals.r 2010-03-31 10:15:52 UTC (rev 75)
+++ pkg/RSienaTest/R/globals.r 2010-04-12 16:01:03 UTC (rev 76)
@@ -139,9 +139,13 @@
Report(c("@", level, "\n", text, "\n"), hdest=dest, sep="", fill=fill)
Report(rep(ch, sum(nchar(text))), hdest=dest, sep="", fill=fill)
if (level < 3)
+ {
Report("\n\n", hdest = dest)
+ }
else
+ {
Report("\n", hdest = dest)
+ }
}
}
@@ -149,10 +153,14 @@
PrtOutMat<- function(mat, dest)
{
if (missing(dest))
- Report(format(t(mat)), sep=c(rep.int(" ", ncol(mat) - 1), "\n"))
+ {
+ Report(format(t(mat), scientific=FALSE),
+ sep=c(rep.int(" ", ncol(mat) - 1), "\n"))
+ }
else
{
- Report(format(t(mat)), sep=c(rep.int(" ", ncol(mat) - 1), "\n"),
+ Report(format(t(mat), scientific=FALSE),
+ sep=c(rep.int(" ", ncol(mat) - 1), "\n"),
hdest=deparse(substitute(dest)))
Report("\n", hdest=deparse(substitute(dest)))
}
Property changes on: pkg/RSienaTest/R/iwlsm.R
___________________________________________________________________
Name: svn:eol-style
+ native
Modified: pkg/RSienaTest/R/maxlike.r
===================================================================
--- pkg/RSienaTest/R/maxlike.r 2010-03-31 10:15:52 UTC (rev 75)
+++ pkg/RSienaTest/R/maxlike.r 2010-04-12 16:01:03 UTC (rev 76)
@@ -1,797 +1,797 @@
-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)
- {
- f <- NULL
- 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,]
- z$theta <- effects$initialValue
- z$fixed <- effects$fix
- z$test <- effects$test
- z$pp <- length(z$test)
- 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$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')
- netsub <- min(netsubs) ### only one for now
- actsubs <- which(types=='behavior')
- for (i in 1:observations)
- {
- mats[[i]] <- data$depvars[[netsub]][, , i]
- f$mynets[[i]] <- mats[[i]]
- if (i==1)
- 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
- }
- f$mystructs <- vector('list',observations)
- for (i in 1:observations)
- {
- f$mystructs[[i]] <- mats[[i]]
- f$mystructs[[i]][, ] <- 0
- f$mystructs[[i]][mats[[i]]==11] <- 1
- f$mystructs[[i]][mats[[i]]==10] <- 1
- }
- f$mats <- mats
- for (i in 1:observations)
- {
- f$mats[[i]][mats[[i]]==11] <- 1
- f$mats[[i]][mats[[i]]==10] <- 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$observations <- observations
- ## browser()
-
- if (any(z$targets!=0))
- {
- 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
- startmat <- mat1
- startmat[is.na(startmat)] <- 0
- endmat <- mat2
- endmat[is.na(endmat)] <- startmat[is.na(endmat)]
- diffmat <- startmat != endmat
- if (is.null(x$multfactor))
- f$niter <- multfactor * sum(diffmat)
- else
- f$niter <- x$multfactor * sum(diffmat)
-### create initial chain
- chain <- matrix(0, nrow=sum(diffmat), ncol=4)
- chain[,1] <- row(diffmat)[diffmat]
- chain[,2] <- col(diffmat)[diffmat]
- chain <- chain[sample(1:nrow(chain)),]
- chain[, 4] <- 1:nrow(chain)
- ##chain<- chain ##(here you can put a known chain in (eg from
- ##delphi!)
- cat(nrow(chain), '\n')
-### initialise
- pinsdel <- pinsdel/(1 - pperm)
- pdiaginsdel <- pdiaginsdel/(1 - pperm)
- iter <- 0
- ##burnin
-###construct a max like object to be passed to FRAN
- f$startmat <- startmat
- f$endmat <- endmat
- f$chain <- chain
- f$accepts <- rep(0,4)
- f$rejects <- rep(0,4)
- f$probs <- c(pinsdel, 0, pdiaginsdel)#
- f$madechain <- FALSE
- f$numm <- 20
- for (i in 1:nstart)
- {
- iter <- iter+1
- ## cat(iter,'\n')
- f <- mhstep(z$theta, f, promul, prelins)
- }
- f$madechain <- TRUE
- pinsdel <- pinsdel * (1-pperm)
- pdiaginsdel <- pdiaginsdel * ( 1-pperm)
- f$probs <- c(pinsdel, pperm, pdiaginsdel)
- f$mats <- f$mystructs <- f$mynets <- NULL
- FRANstore(f)
- z
- }
-
- ##start of function
- scores <- NULL
- dff <- NULL
- theta <- z$theta
- ## f<- z$f
- if (INIT)
- {
- z <- mlInit(z, x, data, effects)
- ## f <<-f
- return(z)
- }
- else if (TERM)
- {
- f <- FRANstore()
- z$f <- f
- return(z)
- }
- else
- {
- f <- FRANstore()
- niter <- f$niter
- nactors <- nrow(f$startmat)
- promul <- promul0
- # int <- x$int
- if (z$Phase==2)
- {
- f$accepts <- rep(0, 4)
- f$rejects <- rep(0, 4)
- varmat <- FALSE
- ## browser()
- if (z$nit == 1)## beginning of a subphase
- {
- i <- 1
- while (i <= 10 * niter)
- {
- f <- mhIntStep(theta, f, promul, prelins, int)
- i <- i + f$n
- }
- }
- Z <- vector("list", noSamples)
- i <- 1
- while (i <= niter)
- {
- f <- mhIntStep(theta, f, promul, prelins, int)
- i <- i + f$n
- }
- Z[[1]] <- f$chain
- if (noSamples > 1)
- {
- for (sample in 2:noSamples)
- {
- i <- 1
- while (i < sampInterval)
- {
- f <- mhIntStep(theta, f, promul, prelins, int)
- i <- i + f$n
- }
- Z[[sample]] <- f$chain
- }
- }
- ans <- calcgrad(theta, Z, f$startmat, varmat)
- # browser()
- f$Z <- Z
- f$chain <- f$Z[[noSamples]]
- }
- else
- {
- varmat <- TRUE
- i <- 1
- while (i <= niter)
- {
- f <- mhIntStep(theta, f, promul, prelins, int)
-
- i <- i + f$n
- }
- ans <- calcgrad(theta, f$chain, f$startmat, varmat)
- }
-
- ## browser()
- scores <- ans$sc
- dff <- ans$dff
- ##cat(object.size(f), object.size(f$chain), object.size(f$startmat), '\n')
- FRANstore(f)
- # cat(scores,'\n')
- ##browser()
- list(fra=matrix(scores, nrow=1), sc=NULL, dff=dff, OK=TRUE,
- rejectprop=f$rejects / (f$accepts + f$rejects))
- }
-}
-mhIntStep <- function(theta, f, promul, prelins, int)
-{
- ff <- lapply(1:int, function(x, theta, f, promul, prelins)
- mhstep(theta, f, promul, prelins),
- theta=theta, f=f, promul=promul,
- prelins=prelins)
- acts <- sapply(ff, function(x)x$act)
- steptypes <- sapply(ff, function(x)x$steptype)
- if (any(acts))
- {
- firstAct <- min(which(acts))
- n <- firstAct
- f$chain <- ff[[n]]$chain
- f$act <- c(rep(FALSE, n-1), TRUE)
- f$steptype <- steptypes[1:n]
- f$accepts[steptypes[n]] <- f$accepts[steptypes[n]]+1
- f$numm <- ff[[n]]$numm
- }
- else
- {
- n <- int
- f$chain <- ff[[n]]$chain
- f$act <- rep(FALSE, n)
- f$steptype <- steptypes
- f$numm <- ff[[n]]$numm
- }
- rejects <- f$steptype[!f$act]
- rejNos <- table(rejects)
- rejs <- as.numeric(names(rejNos))
- f$rejects[rejs] <- f$rejects[rejs] + as.vector(rejNos)
- f$n <- n
- f
-}
-mhstep <- function(theta, f, promul, prelins)
-{
- mhtryinsertdiag<- function(i,kpos)
- {
- tmpnet<- startmat
- mychain<- chain[1:(kpos-1),,drop=FALSE]
- ## tmpnet[mychain[,1:2]]<-
- ## 1-tmpnet[mychain[,1:2]]
- for (mysub in 1:nrow(mychain))
- {
- ii <- mychain[mysub,1]
- jj <- mychain[mysub,2]
- 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])
- ## ans <- kappasigmamu(nactors,nrow(chain),lambda,add1=TRUE)
- ## mu<- ans$mu + 1/lambda/nactors
- ## sigma2 <- ans$sigma2 + 1/lambda/lambda/nactors/nactors
- prr<- pr*lambda*nactors/(ndiag+1 )
- ## prdelphi <- log(pr) - log(nactors) - ans$kappa - 0.5*log(sigma2) -
- ## (1-mu)^2/2/sigma2
- ## prr <- exp(prdelphi) * (nrow(chain)+1)/(ndiag+1) * nactors
-##cat(prr,'\n')
- if (runif(1)<prr)
- accept<-TRUE
- else
- accept<-FALSE
- accept
- }
- mhtrycanceldiag<- function(i,dpos)
- {
- tmpnet<- startmat
- mychain<- chain[1:(dpos-1),,drop=FALSE]
- ## tmpnet[mychain[,1:2]]<-
- ## 1-tmpnet[mychain[,1:2]]
- for (mysub in 1:nrow(mychain))
- {
- ii <- mychain[mysub,1]
- jj <- mychain[mysub,2]
- tmpnet[ii,jj] <- 1- tmpnet[ii,jj]
- }
- diag(tmpnet)<-0
- ps<-calcprobs(i,tmpnet,betapar,nactors)
- pr<- 1/sum(ps)
- ## browser()
- ndiag<-nrow(chain[chain[,1]==chain[,2],,drop=FALSE])
- prr<- ndiag/pr/lambda/nactors
- ## ans <- kappasigmamu(nactors,nrow(chain+1),lambda,add1=TRUE)
- ## mu<- ans$mu - 1/lambda/nactors
- ## sigma2 <- ans$sigma2 - 1/lambda/lambda/nactors/nactors
- ## prdelphi <- -log(pr) + log(nactors) - ans$kappa - 0.5*log(sigma2) -
- ## (1-mu)^2/2/sigma2
- ## prr <- exp(prdelphi) / (nrow(chain))*(ndiag) / nactors
-##cat(prr,'\n')
- if (runif(1)<prr)
- accept<-TRUE
- else
- accept<-FALSE
- accept
- }
- mhinsdelperm<- function(insdel)
- {
- findCCPs <- function(mat)
- {
- if (nrow(mat) > 0)
- {
- dontuse <- c(diff(mat[, 4]), 0) == 1
- ccps <- mat[!dontuse, , drop=FALSE]
- nc <- nrow(ccps) - 1
- }
- else
- {
- ccps <- NULL
- nc <- 0
- }
- list(ccps=ccps, nc=nc)
- }
- ##fix up 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)
- if (is.null(nmul))
- nmul <- 0
- if (nmul>0 && runif(1)<promul)
- {
- ##choose one of unique duplicates
- mypair <- mults[sample(1:nmul,1),]
- }
- else
- {
- mypair <- sample(1:nactors,2)
- }
- from <- mypair[1]
- to <- mypair[2]
- 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
- else
- {
- if (runif(1) < prelins)
- ins <- TRUE
- else
- ins <- FALSE
- }
- del <- !ins
- if (ins)
- {
- if (nk==0)
- {
- kmypair<- sample(1:(nrow(chain)+1),2)
- numav1<- nrow(chain)+1
- numav2<- nrow(chain)
- k1<- min(kmypair)
- k2<- max(kmypair)
- ## newsub<- 1
- }
- 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])
- ## }
- ## else
- ## 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)
- ## 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))
- if (k>nrow(chain)) ##not possible to proceed
- return(list(chain=chain,accept=FALSE))
- numav2<- kk-k
-
- k2<- sample(1:numav2,1)
- if (k2==k1)
- {
- k2 <- kk
- }
- if (k2<k1)
- {
- kk<- k2
- k2<- k1
- k1<- kk
- }
- }
- }
- else ##del
- {
- ## possible<- chain[chain[,1]==from&chain[,2]==to&chain[,3]>0,]
- ## subs<- unique(possible[,3])
- subs <- 1:nc
- if (length(subs)>1)
- ## kmypair<- sample(subs,1)
- kccp <- sample(subs, 1)
- else
- kccp <- subs
- ## kmypair<- subs
- pair <- ccps[kccp : (kccp + 1), 4]
- ## pair<- chain[chain[,1]==from&chain[,2]==to&
- ## chain[,3]==kmypair,4]
- k1<- min(pair)
- k2<- max(pair)
- numav1<-nrow(chain)+1-nrow(similar)
- tempchain <- chain[-pair, ]
- thiskind <- (1:nrow(tempchain))[from == tempchain[, 1] &
- to == tempchain[, 2]]
- ## thiskind<- (1:nrow(chain))[from==chain[,1]&to==chain[,2]]
- k<- max(c(0,thiskind[thiskind<k1]))
- kk<- min(c(thiskind[thiskind>k1],nrow(chain)+1))
- if (k==0 &kk>nrow(chain))
- numav2<- numav1-1
- else
- numav2 <- kk-k-3
- ## numav2<-kk-k-2
- }
-
- } #if insdel
- else
- {
- ins<- FALSE
- del<- FALSE
- k2<- nrow(chain)+1
- k1<- sample(1:nrow(chain),1)
- }
- ##set up network to just before k1
- tmpnet<- startmat
- if (k1>2)
- {
- mychain<- chain[1:(k1-1),,drop=FALSE]
- ## tmpnet[mychain[,1:2]]<-
- ## 1-tmpnet[mychain[,1:2]]
- for (mysub in 1:nrow(mychain))
- {
- ii <- mychain[mysub,1]
- jj <- mychain[mysub,2]
- tmpnet[ii,jj] <- 1- tmpnet[ii,jj]
- }
- diag(tmpnet)<- 0
- }
- k1mat<- tmpnet
- ##decide on permutation
- ## permute num steps but not including k2, stop before end
- ##if deleting k1, remove that one first
- if (del)
- {
- num<- min(num,k2-k1-1)
- end<- k1+num
- start<- k1+1
- }
- else
- {
- num <- min(num,k2-k1)
- end <- k1+num-1
- start <- k1
- }
- truncated <- num!=trunc(numm)
- ##calculate probs of existing chain from k1 to end or k2
- if (ins)
- thisend<- k2-1
- else if (del)
- thisend<- k2
- else
- thisend<- end
- prbefore<- rep(0,thisend-k1+1)
- for (link in k1:thisend)
- {
- i<- chain[link,1]
- j<- chain[link,2]
- ps<- calcprobs(i,tmpnet,betapar,nactors)
- prbefore[link-k1+1]<- ps[j]/sum(ps)
- if (i!=j)
- tmpnet[i,j]<- 1-tmpnet[i,j]
- }
- ##reset up network to just before k1
- tmpnet<- k1mat
- if (num > 0)
- {
- permchain<- chain[start:end,,drop=FALSE]
-
- myorder<- sample(1:nrow(permchain))
- permchain<- permchain[myorder,]
- }
- else
- permchain <- NULL
- ##reconstruct new chain
- if (k1>1)
- tempchain<- chain[1:(k1-1),]
- else
- tempchain<- NULL
- if (ins)
- tempchain<- rbind(tempchain,c(from,to,0,0))
- ## tempchain<- rbind(tempchain,c(from,to,newsub,0))
- tempchain<- rbind(tempchain,permchain)
- if (end<(k2-1))
- tempchain<- rbind(tempchain,chain[(end+1):(k2-1),]) ##check this
- if (ins)
- tempchain<- rbind(tempchain,c(from,to,0,0))
- ## tempchain<- rbind(tempchain,c(from,to,newsub,0))
- if (!del & k2<=nrow(chain))
- tempchain<- rbind(tempchain, chain[k2:nrow(chain),])
- if (del &k2 < nrow(chain))
- tempchain<- rbind(tempchain, chain[(k2+1):nrow(chain),])
- ##calc all probs between k1 and k2 for new chain
- if (ins)
- {
- start<- k1
- end<- k2+1
- }
- else
- if (del)
- {
- start<- k1
- end<- k2-2
- }
- else
- {
- start<- k1
- end<- k1+num-1
- }
- if (end-start+1>0)
- {
- prafter<- rep(0,end-start+1)
- for (link in start:end)
- {
- i<- tempchain[link,1]
- j<- tempchain[link,2]
- ps<- calcprobs(i,tmpnet,betapar,nactors)
- prafter[link-start+1]<- ps[j]/sum(ps)
- if (i!=j)
- tmpnet[i,j]<- 1-tmpnet[i,j]
- }
- }
- else
- prafter<- 1
- ##now try to get the ratios
- ## if(del) browser()
- ## ans<- kappasigmamu(nactors,nrow(chain),lambda,add1=TRUE)
- ## if (ins)
- ## {
- ## mu<- ans$mu+2/nactors/lambda
- ## sigma2 <- ans$sigma2+2/nactors/nactors/lambda/lambda
- ## prdelphi <- sum(log(prafter)) - sum(log(prbefore)) -
- ## 2 * log(nactors) -ans$kappa - 0.5 * log(sigma2) -
- ## (1-mu)^2/2/sigma2
- ## prdelphi <- exp(prdelphi)
- ## }
- ## if (del)
- ## {
- ## mu<- ans$mu - 2/nactors/lambda
- ## sigma2 <- ans$sigma2 - 2/nactors/nactors/lambda/lambda
- ## prdelphi <- sum(log(prafter))-sum(log(prbefore)) +
- ## 2 * log(nactors) - ans$kappa - 0.5 * log(sigma2) -
- ## (1-mu)^2/2/sigma2
- ## prdelphi <- exp(prdelphi)
- ## }
- ## if (!insdel)
- ## {
- ## prdelphi<- sum(log(prafter))-sum(log(prbefore))
- ## }
- logprobrat<- 0
- if (ins)
- logprobrat<- 2*log(lambda)-log(nrow(chain)+2)-
- log(nrow(chain)+1)
- if (del)
- logprobrat<- -2*log(lambda)+log(nrow(chain))+
- log(nrow(chain)-1)
- logprobrat<- logprobrat - sum(log(prbefore))
- if (length(prafter)>0)
- logprobrat<- logprobrat +
- sum(log(prafter))
- probrat<- exp(logprobrat)
- pa<- (1-promul)/nactors/(nactors-1)
- ## tempchain[,4] <- 1:nrow(tempchain)
- if (insdel)
- {
- newsimilar <- tempchain[tempchain[, 1] == from &
- tempchain[, 2] == to, ,
- drop=FALSE]
- newnc <- findCCPs(newsimilar)$nc
- }
- if (ins)
- {
- if (nk<2)
- pp<- (promul/(nmul+1) +pa)/pa*(1-prelins)*
- ## numav1*numav2/2/(nc+1)
- numav1*numav2/2/(newnc)
- else
- pp<- (1-prelins)/numav1*numav2/2/(newnc)
- ## pp<- (1-prelins)/numav1*numav2/2/(nc+1)
- if (nc>=1) pp <- pp/prelins
- }
- else if (del)
- {
- if (nk<4)
- pp<- pa*2*nc/(promul/nmul +pa)/
- (1-prelins)/numav1/numav2
- else
- pp<- 1/(1-prelins)/numav1/numav2*2*nc
- ## if (nc>=2)
- if (newnc >= 1)
- pp <- pp*prelins
- }
- else
- {
- pp<- 1
- }
- probrat <- probrat*pp
- ##cat('probrat ',probrat1,' ')
- ## probrat <- prdelphi * pp
- ## cat(probrat,' ',probrat-probrat1,'\n')
- 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
- if (runif(1)< probrat)
- accept<- TRUE
- if (sum(tempchain[,3]==1)%%2!=0)
- browser()
- if (accept)
- chain<- tempchain
- ## cat(num,f$numm,accept , insdel,truncated,accept,'\n')
- ## if (!insdel)
- ## browser()
- if (accept && !insdel && !truncated)
- numm<- numm+0.5
- if (!accept && !insdel && !truncated)
- numm <- numm-0.5
- #browser()
- list(chain=chain,accept=accept,numm=numm)
- }##end of procedure
- #########################################################################
- ##start of mhstep
- #########################################################################
- #cat('start', f$numm,'\n')
- ## print(table(f$chain))
- startmat<- f$startmat
- nactors<- nrow(startmat)
- chain<- f$chain
- lambda<- theta[1]
- betapar<- theta[-1]
- steptype<- sample(1:3,1,prob=f$probs)
- # cat(steptype,'step type\n')
- if (steptype==1)
- {
- ans<- mhinsdelperm(TRUE)
- act<- ans$accept
- chain<- ans$chain
- }
- else if (steptype==2)
- {
- if (f$madechain)
- {
- ans<-mhinsdelperm(FALSE)
- act<- ans$accept
- chain<- ans$chain
- f$numm <- ans$numm
- }
- else
- {
- act<- NA
- }
- }
- else ##if (steptype==3)
- {
- actor <- sample(1:nactors,1)
- kpos<- sample(1:(nrow(chain)+1),1)
- act<- mhtryinsertdiag(actor,kpos)
-
- if (act)
- {
- if (kpos>1)
- tmp<- chain[1:(kpos-1),]
- else
- tmp<- NULL
- tmp<- rbind(tmp,c(actor,actor,0,0))
- if (kpos<=nrow(chain))
- tmp<- rbind(tmp, chain[kpos:nrow(chain),])
- ##if (sum(chain[,3]==1)%%2!=0)
- ## browser()
- chain <- tmp
- chain[, 4] <- 1:nrow(chain)
- }
- ## }
- ## else
- ## {
- diags<-chain[chain[,1]==chain[,2],,drop=FALSE]
- ndiag<- nrow(diags)
- act<- FALSE
- if (ndiag>0)
- {
- ## browser()
- kk<- sample(1:ndiag,1)
- actor<- diags[kk,1]
- dpos<- diags[kk,4]
- ## if(chain[dpos,1]!=actor) browser()
- act <- mhtrycanceldiag(chain[dpos,1],dpos)
- if (act)
- {
- if (dpos>1)
- tmp<- chain[1:(dpos-1),]
- else
- tmp<- NULL
- if (dpos<nrow(chain))
- tmp<- rbind(tmp, chain[(dpos+1):nrow(chain),])
- if (sum(chain[,3]==1)%%2!=0)
- browser()
- chain <- tmp
- }
- }
- }
- ## cat('act',act,'steptype',steptype,' ')
- chain[, 4] <- 1:nrow(chain)
- f$chain <- chain
- f$act <- act
- f$steptype <- steptype
- #cat(f$numm,'\n')
- f
-}
-
-calcprobs <- function(i, tmpnet, betapar, nactors)
-{
- ss <- matrix(0, nrow=nactors, ncol=2)
- myout <- tmpnet[i, ]
- myin <- tmpnet[, i]
- ss[, 1] <- 1 - 2 * myout
- ss[, 2] <- (1 - 2 * myout) * myin
- ss[i, ] <- 0
- fs <- drop(ss %*% betapar)
- exp(fs)
-}
-
-doMoreUpdates <- function(z, x, n)
-{
- f <- FRANstore()
- for (i in 1:n)
- {
- ans <- calcgrad(z$theta, f$Z, f$startmat, varmat=FALSE)
-
- fchange <- as.vector(z$gain * ans$sc %*% z$dinv)
- fchange <- ifelse(z$posj & (fchange > z$theta), z$theta * 0.5, fchange)
-
- z$theta <- z$theta - fchange
- z$thav <- z$thav + z$theta
- z$thavn <- z$thavn + 1
- }
- z
-}
+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)
+ {
+ f <- NULL
+ 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,]
+ z$theta <- effects$initialValue
+ z$fixed <- effects$fix
+ z$test <- effects$test
+ z$pp <- length(z$test)
+ 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$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')
+ netsub <- min(netsubs) ### only one for now
+ actsubs <- which(types=='behavior')
+ for (i in 1:observations)
+ {
+ mats[[i]] <- data$depvars[[netsub]][, , i]
+ f$mynets[[i]] <- mats[[i]]
+ if (i==1)
+ 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
+ }
+ f$mystructs <- vector('list',observations)
+ for (i in 1:observations)
+ {
+ f$mystructs[[i]] <- mats[[i]]
+ f$mystructs[[i]][, ] <- 0
+ f$mystructs[[i]][mats[[i]]==11] <- 1
+ f$mystructs[[i]][mats[[i]]==10] <- 1
+ }
+ f$mats <- mats
+ for (i in 1:observations)
+ {
+ f$mats[[i]][mats[[i]]==11] <- 1
+ f$mats[[i]][mats[[i]]==10] <- 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$observations <- observations
+ ## browser()
+
+ if (any(z$targets!=0))
+ {
+ 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
+ startmat <- mat1
+ startmat[is.na(startmat)] <- 0
+ endmat <- mat2
+ endmat[is.na(endmat)] <- startmat[is.na(endmat)]
+ diffmat <- startmat != endmat
+ if (is.null(x$multfactor))
+ f$niter <- multfactor * sum(diffmat)
+ else
+ f$niter <- x$multfactor * sum(diffmat)
+### create initial chain
+ chain <- matrix(0, nrow=sum(diffmat), ncol=4)
+ chain[,1] <- row(diffmat)[diffmat]
+ chain[,2] <- col(diffmat)[diffmat]
+ chain <- chain[sample(1:nrow(chain)),]
+ chain[, 4] <- 1:nrow(chain)
+ ##chain<- chain ##(here you can put a known chain in (eg from
+ ##delphi!)
+ cat(nrow(chain), '\n')
+### initialise
+ pinsdel <- pinsdel/(1 - pperm)
+ pdiaginsdel <- pdiaginsdel/(1 - pperm)
+ iter <- 0
+ ##burnin
+###construct a max like object to be passed to FRAN
+ f$startmat <- startmat
+ f$endmat <- endmat
+ f$chain <- chain
+ f$accepts <- rep(0,6)
+ f$rejects <- rep(0,6)
+ f$probs <- c(pinsdel, 0, pdiaginsdel)#
+ f$madechain <- FALSE
+ f$numm <- 20
+ for (i in 1:nstart)
+ {
+ iter <- iter+1
+ ## cat(iter,'\n')
+ f <- mhstep(z$theta, f, promul, prelins)
+ }
+ f$madechain <- TRUE
+ pinsdel <- pinsdel * (1-pperm)
+ pdiaginsdel <- pdiaginsdel * ( 1-pperm)
+ f$probs <- c(pinsdel, pperm, pdiaginsdel)
+ f$mats <- f$mystructs <- f$mynets <- NULL
+ FRANstore(f)
+ z
+ }
+
+ ##start of function
+ scores <- NULL
+ dff <- NULL
+ theta <- z$theta
+ ## f<- z$f
+ if (INIT)
+ {
+ z <- mlInit(z, x, data, effects)
+ ## f <<-f
+ return(z)
+ }
+ else if (TERM)
+ {
+ f <- FRANstore()
+ z$f <- f
+ return(z)
+ }
+ else
+ {
+ f <- FRANstore()
+ niter <- f$niter
+ nactors <- nrow(f$startmat)
+ promul <- promul0
+ # int <- x$int
+ if (z$Phase==2)
+ {
+ f$accepts <- rep(0, 6)
+ f$rejects <- rep(0, 6)
+ varmat <- FALSE
+ ## browser()
+ if (z$nit == 1)## beginning of a subphase
+ {
+ i <- 1
+ while (i <= 10 * niter)
+ {
+ f <- mhIntStep(theta, f, promul, prelins, int)
+ i <- i + f$n
+ }
+ }
+ Z <- vector("list", noSamples)
+ i <- 1
+ while (i <= niter)
+ {
+ f <- mhIntStep(theta, f, promul, prelins, int)
+ i <- i + f$n
+ }
+ Z[[1]] <- f$chain
+ if (noSamples > 1)
+ {
+ for (sample in 2:noSamples)
+ {
+ i <- 1
+ while (i < sampInterval)
+ {
+ f <- mhIntStep(theta, f, promul, prelins, int)
+ i <- i + f$n
+ }
+ Z[[sample]] <- f$chain
+ }
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/rsiena -r 76
More information about the Rsiena-commits
mailing list