Nothing
#/******************************************************************************
# * SIENA: Simulation Investigation for Empirical Network Analysis
# *
# * Web: http://www.stats.ox.ac.uk/~snijders/siena
# *
# * File: maxlikec.r
# *
# * Description: This module contains the code for simulating the process,
# * communicating with C++. For use with maximum likelihood method, so
# * never conditional or from finite differences, or parallel testing.
# *****************************************************************************/
##@maxlikec siena07 ML Simulation Module
maxlikec <- function(z, x, data=NULL, effects=NULL,
returnChains=FALSE, byGroup=FALSE, byWave=FALSE,
returnDataFrame=FALSE,
returnLoglik=FALSE, onlyLoglik=FALSE)
{
# note: parameter x is not used. Just for consistency with other possibilities for FRAN.
## retrieve stored information
f <- FRANstore()
callGrid <- z$callGrid
if (nrow(callGrid) == 1)
{
if (byGroup)
{
theta <- z$thetaMat[1,]
}
else
{
theta <- z$theta
}
ans <- .Call(C_mlPeriod, PACKAGE=pkgname, z$Deriv, f$pData,
f$pModel, f$myeffects, theta,
1, 1, z$nrunMH, z$addChainToStore,
z$returnDataFrame, z$returnDeps,
z$returnChains, returnLoglik, onlyLoglik)
if (!onlyLoglik)
{
ans[[6]] <- list(ans[[6]])
ans[[7]] <- list(ans[[7]])
if (byGroup)
{
ans[[8]] <- list(ans[[8]])
ans[[9]] <- list(ans[[9]])
ans[[10]] <- list(ans[[10]])
ans[[11]] <- list(ans[[11]])
}
}
else
{
ans[[2]] <- list(ans[[2]])
ans[[3]] <- list(ans[[3]])
ans[[4]] <- list(ans[[4]])
}
if (z$returnDeps)
{
sims <- ans[[12]]
}
else
{
sims <- 'there are no simulated dependent variables'
}
}
else
{
## z$int2 is the number of processors if iterating by period, so 1 means
## we are not. Can only parallelize by period withmaxlike.
if (z$int2 == 1)
{
anss <- apply(cbind(callGrid, 1:nrow(callGrid)),
1, doMLModel, z$Deriv, z$thetaMat,
z$nrunMH, z$addChainToStore,
z$returnDataFrame, z$returnDeps,
z$returnChains, byGroup, z$theta, returnLoglik,
onlyLoglik)
}
else
{
use <- 1:(min(nrow(callGrid), z$int2))
anss <- parRapply(z$cl[use], cbind(callGrid, 1:nrow(callGrid)),
doMLModel, z$Deriv, z$thetaMat,
z$nrunMH, z$addChainToStore,
z$returnDataFrame, z$returnDeps, z$returnChains, byGroup,
z$theta, returnLoglik, onlyLoglik)
}
## reorganize the anss so it looks like the normal one
ans <- list()
if (!onlyLoglik)
{
ans[[1]] <- sapply(anss, "[[", 1) ## statistics
ans[[2]] <- NULL ## scores
ans[[3]] <- NULL ## seeds
ans[[4]] <- NULL ## ntim
ans[[5]] <- NULL # randomseed
if (z$returnChains)
{
fff <- lapply(anss, function(x) x[[6]][[1]])
fff <- split(fff, callGrid[, 1 ]) ## split by group
ans[[6]] <- fff
}
ans[[7]] <- lapply(anss, "[[", 7) ## derivative
ans[[8]] <- lapply(anss, "[[", 8)
ans[[9]] <- lapply(anss, "[[", 9)
ans[[10]] <- lapply(anss, "[[", 10)
if (!byGroup)
{
ans[[8]] <- Reduce("+", ans[[8]]) ## accepts
ans[[9]] <- Reduce("+", ans[[9]]) ## rejects
ans[[10]] <- Reduce("+", ans[[10]]) ## aborts
}
ans[[11]] <- sapply(anss, "[[", 11)
if (z$returnDeps)
{
fff <- lapply(anss, function(x) x[[12]])
sims <- split(fff, callGrid[, 1 ]) ## split by group
}
else
{
sims <- 'no simulated dependent variables'
}
}
else ##onlyLoglik is always byGroup (sienaBayes)
{
ans[[1]] <- sum(sapply(anss, '[[', 1))## loglik
ans[[2]] <- lapply(anss, "[[", 2)
ans[[3]] <- lapply(anss, "[[", 3)
ans[[4]] <- lapply(anss, "[[", 4)
}
}
FRANstore(f)
if (z$Deriv && !onlyLoglik) ## need to reformat the derivatives
{
resp <- reformatDerivs(z, f, ans[[7]])
dff <- resp[[1]]
dff2 <- resp[[2]]
}
else
{
dff <- NULL
dff2 <- NULL
}
if (!onlyLoglik)
{
fra <- -t(ans[[1]]) ##note sign change
list(fra=fra, ntim0=NULL, feasible=TRUE, OK=TRUE,
sims=sims, dff=dff, dff2=dff2,
chain=ans[[6]], accepts=ans[[8]],
rejects=ans[[9]], aborts=ans[[10]], loglik=ans[[11]])
}
else
{
list(loglik=ans[[1]], accepts=ans[[2]],
rejects=ans[[3]], aborts=ans[[4]])
}
}
##@doMLModel Maximum likelihood
doMLModel <- function(x, Deriv, thetaMat, nrunMH, addChainToStore,
returnDataFrame, returnDeps, returnChains,
byGroup, theta, returnLoglik, onlyLoglik)
{
f <- FRANstore()
if (byGroup)
{
theta <- thetaMat[x[1], ]
}
else
{
theta <- theta
}
.Call(C_mlPeriod, PACKAGE=pkgname, Deriv, f$pData,
f$pModel, f$myeffects, theta,
as.integer(x[1]), as.integer(x[2]), nrunMH[x[3]], addChainToStore,
returnDataFrame, returnDeps, returnChains,
returnLoglik, onlyLoglik)
}
##@reformatDerivs Maximum likelihood move this back to inline or internal function
reformatDerivs <- function(z, f, derivList)
{
## current format is a list of vectors of the lower? triangle
## of the matrix. Need to be put into a symmetric matrix.
## Tricky part is getting the rates in the right place!
## Note that we do not yet deal with rate effects other than the basic
dff <- as(matrix(0, z$pp, z$pp), "dsyMatrix")
nPeriods <- length(derivList)
dff2 <- vector("list", nPeriods)
if (z$byWave)
{
tmp <- as(matrix(0, z$pp, z$pp), "dsyMatrix")
dff2[] <- tmp
}
for (period in 1:nPeriods)
{
dffraw <- derivList[[period]]
## start indexes row/col for first effect for the variable
start <- 1
## rawsub is subscript in the vector
rawsub <- 1
## f$myeffects is a list of data frames per dependent variable
## of selected effects.
for (i in 1:length(f$myeffects))
{
dffPeriod <- matrix(0, nrow=z$pp, ncol=z$pp)
## rows is the number of effects for this variable
rows <- nrow(f$myeffects[[i]])
## nRates is the number of rate effects for this variable
## at present nRates will be the number of periods.
nRates <- sum(f$myeffects[[i]]$type == 'rate')
## nonRates is the number of rows/cols in the objective function
## part of the derivative matrix dff.
nonRates <- rows - nRates
## first put the basic rate for this variable in the right place
dffPeriod[cbind(start:(start + nRates -1),
start:(start + nRates - 1))] <-
dffraw[rawsub:(rawsub + nRates - 1)]
##
## now the matrix of objective function effects
##
start <- start + nRates
##
rawsub <- rawsub + nRates
##
dffPeriod[start : (start + nonRates - 1 ),
start : (start + nonRates - 1)] <-
dffraw[rawsub:(rawsub + nonRates * nonRates - 1)]
start <- start + nonRates
rawsub <- rawsub + nonRates * nonRates
dffPeriod <- dffPeriod + t(dffPeriod)
diag(dffPeriod) <- diag(dffPeriod) / 2
if (z$byWave)
{
dff2[[period]] <- dff2[[period]] - dffPeriod
}
dff <- dff - dffPeriod
}
}
list(dff, dff2)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.