[Rsiena-commits] r81 - in pkg: RSiena RSiena/R RSiena/inst/doc RSiena/man RSiena/src RSiena/src/data RSiena/src/model RSiena/src/model/effects RSiena/src/model/ml RSiena/src/model/variables RSiena/src/network RSiena/src/utils RSienaTest RSienaTest/R RSienaTest/doc RSienaTest/inst/doc RSienaTest/src RSienaTest/src/model RSienaTest/src/model/effects RSienaTest/src/model/ml RSienaTest/src/model/variables RSienaTest/src/utils
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat Apr 24 17:09:31 CEST 2010
Author: ripleyrm
Date: 2010-04-24 17:09:30 +0200 (Sat, 24 Apr 2010)
New Revision: 81
Added:
pkg/RSiena/R/bayes.r
pkg/RSiena/R/maxlike.r
pkg/RSiena/R/maxlikec.r
pkg/RSiena/R/maxlikecalc.r
pkg/RSiena/man/maxlikefn.Rd
pkg/RSienaTest/R/bayes.r
Modified:
pkg/RSiena/R/phase2.r
pkg/RSiena/R/phase3.r
pkg/RSiena/R/simstatsc.r
pkg/RSiena/changeLog
pkg/RSiena/inst/doc/s_man400.pdf
pkg/RSiena/src/data/ActorSet.cpp
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/effects/AverageReciprocatedAlterEffect.cpp
pkg/RSiena/src/model/effects/Effect.cpp
pkg/RSiena/src/model/effects/Effect.h
pkg/RSiena/src/model/effects/ReciprocalDegreeBehaviorEffect.cpp
pkg/RSiena/src/model/effects/StructuralRateEffect.cpp
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/EffectValueTable.cpp
pkg/RSiena/src/model/variables/NetworkVariable.cpp
pkg/RSiena/src/model/variables/NetworkVariable.h
pkg/RSiena/src/network/OneModeNetwork.cpp
pkg/RSiena/src/siena07.cpp
pkg/RSiena/src/utils/Random.cpp
pkg/RSiena/src/utils/Random.h
pkg/RSiena/src/utils/SqrtTable.cpp
pkg/RSienaTest/R/phase3.r
pkg/RSienaTest/R/sienaTimeTest.r
pkg/RSienaTest/R/simstatsc.r
pkg/RSienaTest/changeLog
pkg/RSienaTest/doc/s_man400.tex
pkg/RSienaTest/inst/doc/s_man400.pdf
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/StructuralRateEffect.cpp
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/variables/BehaviorVariable.cpp
pkg/RSienaTest/src/model/variables/DependentVariable.cpp
pkg/RSienaTest/src/model/variables/DependentVariable.h
pkg/RSienaTest/src/model/variables/NetworkVariable.cpp
pkg/RSienaTest/src/siena07.cpp
pkg/RSienaTest/src/utils/Random.cpp
pkg/RSienaTest/src/utils/Random.h
Log:
Bug fixes, ML, Bayes
Added: pkg/RSiena/R/bayes.r
===================================================================
--- pkg/RSiena/R/bayes.r (rev 0)
+++ pkg/RSiena/R/bayes.r 2010-04-24 15:09:30 UTC (rev 81)
@@ -0,0 +1,254 @@
+##/*****************************************************************************
+## * SIENA: Simulation Investigation for Empirical Network Analysis
+## *
+## * Web: http://www.stats.ox.ac.uk/~snidjers/siena
+## *
+## * File: bayes.r
+## *
+## * Description: This file contains the code to run Bayesian simulation
+## *
+## ****************************************************************************/
+##@bayes algorithm fit a Bayesian model
+bayes <- function(data, effects, model, nwarm=100, nmain=100, nrunMHBatches=20,
+ nrunMH=100)
+{
+ createStores <- function()
+ {
+ 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=6)
+ z$MHrejections <- matrix(NA, nrow=nmain * nrunMHBatches , ncol=6)
+ z$MHproportions <- matrix(NA, nrow=nmain * nrunMHBatches, ncol=6)
+ z
+ }
+ storeData <- function()
+ {
+ start <- z$sub + 1
+ nrun <- nrow(z$parameters)
+ end <- start + nrun - 1
+ z$accepts <- as.logical(z$accepts)
+ z$acceptances[start:end] <- z$accepts
+ z$lambdas[start:end, ] <- z$parameters[, basicRate]
+ z$lambdadist[start:end, ] <- z$shapes[, basicRate]
+ z$candidates[start:end, ] <- z$parameters[, !basicRate]
+ z$parameters[!z$accepts, !basicRate] <- NA
+ z$parameters[, !basicRate] <-
+ carryForward(z$parameters[, !basicRate], z$betas[1:end,])
+ z$betas[start:end, ] <- z$parameters[, !basicRate]
+
+ z$posteriorTot <- z$posteriorTot + colSums(z$parameters)
+ for (i in npar)
+ {
+ z$posteriorMII <- z$posteriorMII +
+ outer(z$parameters[i, ], z$parameters[i, ])
+ }
+ z$MHacceptances[start:end, ] <- z$MHaccepts
+ z$MHrejections[start:end, ] <- z$MHrejects
+ z$MHproportions[start:end, ] <-z$MHaccepts/ (z$MHaccepts + z$MHrejects)
+ z$sub <- z$sub + nrun
+ z
+ }
+
+ carryForward <- function(parameters, betas)
+ {
+ npar <- nrow(parameters)
+ nbeta <- nrow(betas)
+ if (npar < nbeta)
+ {
+ parameters <- rbind(betas[(nbeta - npar), ], parameters)
+ }
+ parameters <-
+ apply(parameters, 2, function(x)
+ {
+ x.pos <- which(!is.na(x))
+ if (length(x.pos) == 0 || x.pos[1] != 1)
+ {
+ x.pos <- c(1, x.pos)
+ }
+ x[rep(x.pos, c(diff(x.pos),
+ length(x) - x.pos[length(x.pos)] + 1))]
+ }
+ )
+ if (npar < nbeta)
+ {
+ parameters[-1, ]
+ }
+ else
+ {
+ parameters
+ }
+ }
+
+ improveMH <- function(z, x, tiny=1.0e-15, desired=40, maxiter=100,
+ tolerance=15)
+ {
+ rescaleCGD <- function(iter)
+ {
+ if (actual > desired)
+ {
+ u <- 2 - ((iter - actual) / (iter - desired))
+ }
+ else
+ {
+ u <- 1 / (2 - (actual / desired))
+ }
+ if (tol <- abs(actual - desired) <= tolerance)
+ {
+ number <<- number + 1
+ if (number == 2) success <<- TRUE
+ }
+ else
+ {
+ number <<- 0
+ }
+ u
+ }
+ iter <- 0
+ number <- 0
+ success <- FALSE
+ repeat
+ {
+ iter <- iter + 1
+ z <- MCMCcycle(z, nrunMH=1, nrunMHBatches=100)
+ actual <- z$BayesAcceptances ## acceptances
+ ans <- rescaleCGD(100)
+ z$scaleFactor <- z$scaleFactor * ans
+ if (success | iter == maxiter)
+ {
+ break
+ }
+ if (z$scaleFactor < tiny)
+ {
+ cat('scalefactor < tiny\n')
+ browser()
+ }
+ }
+ cat('fine tuning took ', iter, ' iterations. Scalefactor:',
+ z$scaleFactor, '\n')
+ z
+ }
+
+ ## initialise
+ Report(openfiles=TRUE, type="n") #initialise with no file
+ require(lattice)
+ ## require(MASS)
+ ## dev.new()
+ ## histPlot = dev.cur()
+ dev.new()
+ thetaplot = dev.cur()
+ dev.new()
+ acceptsplot = dev.cur()
+ z <- NULL
+ z$FinDiff.method <- FALSE
+ z$maxlike <- TRUE
+ model$maxlike <- TRUE
+ model$FRANname <- "maxlikec"
+ 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)
+ numm <- z$numm
+ }
+
+ for (ii in 1:nmain)
+ {
+ z <- MCMCcycle(z, nrunMH=100, nrunMHBatches=20)
+ z <- storeData()
+
+ numm <- z$numm
+ if (ii %% 10 == 0) ## do some plots
+ {
+ cat('main after ii',ii,numm, '\n')
+ ## dev.set(histPlot)
+ ## par(mfrow=c(2,3))
+ ## truehist(z$lambdas)
+ ## truehist(z$betas[, 1])
+ ## truehist(z$candidates[, 1])
+ ## truehist(z$betas[, 2])
+ ## truehist(z$candidates[, 2])
+ dev.set(thetaplot)
+ thetadf <- data.frame(z$lambdas, z$betas)
+ acceptsdf <- data.frame(z$MHproportions, z$acceptances)
+ lambdaNames <- paste(z$effects$name[basicRate],
+ z$effects$shortName[basicRate],
+ z$effects$period[basicRate],
+ z$effects$group[basicRate], sep=".")
+ betaNames <- paste(z$effects$name[!basicRate],
+ z$effects$shortName[!basicRate], sep=".")
+ names(thetadf) <- c(lambdaNames, betaNames)
+ names(acceptsdf) <- c("InsDiag", "CancDiag", "Permute", "InsPerm",
+ "CancPerm", "Missing", "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))
+ 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"))
+ }
+ }
+ z
+}
+
+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$pMLSimulation, f$myeffects, as.integer(period),
+ as.integer(group),
+ z$scaleFactor, nrunMH, nrunMHBatches)
+ ## process the return values
+ z$BayesAcceptances <- sum(ans[[1]])
+ z$accepts <- ans[[1]]
+ z$parameters <- ans[[2]]
+ z$shapes <- ans[[3]]
+ z$MHaccepts <- ans[[4]]
+ z$MHrejects <- ans[[5]]
+ z
+}
+
+
+
Property changes on: pkg/RSiena/R/bayes.r
___________________________________________________________________
Name: svn:eol-style
+ native
Added: pkg/RSiena/R/maxlike.r
===================================================================
--- pkg/RSiena/R/maxlike.r (rev 0)
+++ pkg/RSiena/R/maxlike.r 2010-04-24 15:09:30 UTC (rev 81)
@@ -0,0 +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,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
+ }
+ }
+ 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), ntim0 = NULL, feasible = TRUE, OK = TRUE, dff=dff, accepts=f$accepts, rejects= 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
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/rsiena -r 81
More information about the Rsiena-commits
mailing list