Nothing
#ntrials is a vector with length equal to length(y). if Bern or Poisson, ntrials is a vec of 1s
objfun <-
function(par, cache, vars){
vars$par <- par
vars$beta<-beta<-vars$par[1:vars$nbeta]
vars$nu<-vars$par[-(1:vars$nbeta)]
if (!missing(cache)) stopifnot(is.environment(cache))
if(any(vars$nu<=0)){
out<-list(value=-Inf,gradient=rep(1,length(vars$par)),hessian=as.matrix(c(rep(1,length(vars$par)^2)),nrow=length(vars$par)))
return(out)
}
vars$Z=do.call(cbind,vars$mod.mcml$z)
vars$Tee<-length(vars$mod.mcml$z)
nrand<-lapply(vars$mod.mcml$z,ncol)
vars$nrandom<-unlist(nrand)
vars$family.glmm<-getFamily(vars$family.glmm)
if(vars$family.glmm$family.glmm=="bernoulli.glmm"){vars$family_glmm=1}
if(vars$family.glmm$family.glmm=="poisson.glmm"){vars$family_glmm=2}
if(vars$family.glmm$family.glmm=="binomial.glmm"){vars$family_glmm=3}
Dstarinvdiag<-1/diag(vars$D.star)
vars$D.star.inv<-diag(Dstarinvdiag)
vars$logdet.D.star.inv<- -sum(log(diag(vars$D.star)))
clusterExport(vars$cl, c("vars"), envir = environment())
clusterEvalQ(vars$cl, logdet.Sigmuh.inv<-sum(log(eigen(umatparams$Sigmuh.inv,symmetric=TRUE)$values)))
vars$myq<-nrow(vars$D.star.inv)
vars$tconst <- tconstant(vars$zeta, vars$myq, Dstarinvdiag)
#for the particular value of nu we're interested in, need to prep for distRandGenC
eek <- getEk(vars$mod.mcml$z)
preDinvfornu <- Map("*",eek,(1/vars$nu))
vars$Dinvfornu <- addVecs(preDinvfornu)
vars$logdetDinvfornu <- sum(log(vars$Dinvfornu))
vars$Dinvfornu <- diag(vars$Dinvfornu)
vars$meow <- rep(1,vars$Tee+1)
vars$meow[1] <- 0
throwaway <- vars$Tee+1
vars$meow[2:throwaway] <- cumsum(vars$nrandom)
vars$pea <- c(vars$p1,vars$p2,vars$p3)
vars$n <- nrow(vars$mod.mcml$x)
##need to scale first m1 vectors of generated random effects by multiplying by A
# preAfornu<-Map("*",eek,sqrt(nu))
# Afornu<-addVecs(preAfornu)
# for(k in 1:m1){
# u.swoop<-umat[k,]
# umat[k,]<-u.swoop*Afornu
# }
clusterEvalQ(vars$cl, C_valgrad <- glmm:::C_valgrad)
clusterEvalQ(vars$cl, C_hess <- glmm:::C_hess)
#parallelizing the calculations for the value of the log-likelihood approximation and gradient
clusterExport(vars$cl, c("vars"), envir = environment()) #installing variables on each core
out <- foreach(i=1:vars$no_cores) %dopar% {
# y is vector of length n
stopifnot(length(vars$mod.mcml$y) == vars$n)
# Umat is myq by m matrix, the R matrix transposed
stopifnot(nrow(t(umatparams$umat)) == vars$myq)
stopifnot(ncol(t(umatparams$umat)) == umatparams$m)
# x is n by nbeta matrix
stopifnot(nrow(vars$mod.mcml$x) == vars$n)
stopifnot(ncol(vars$mod.mcml$x) == vars$nbeta)
# beta is vector of length nbeta
stopifnot(length(vars$beta) == vars$nbeta)
# z is n by myq matrix
stopifnot(nrow(vars$Z) == vars$n)
stopifnot(ncol(vars$Z) == vars$myq)
# Dinvfornu is myq x myq
stopifnot(nrow(vars$Dinvfornu) == vars$myq)
stopifnot(ncol(vars$Dinvfornu) == vars$myq)
# logdetDinvfornu is scalar
stopifnot(length(vars$logdetDinvfornu) == 1)
# only YOU can prevent segfault errors
# prevent segfault errors by checking dims
stopifnot(length(vars$mod.mcml$y) == vars$n)
stopifnot(nrow(vars$mod.mcml$x) == vars$n)
stopifnot(ncol(vars$mod.mcml$x) == vars$nbeta)
stopifnot(length(vars$beta) == vars$nbeta)
stopifnot(nrow(vars$Z) == vars$n)
stopifnot(ncol(vars$Z) == vars$myq)
stopifnot(nrow(vars$Dinvfornu) == vars$myq)
stopifnot(ncol(vars$Dinvfornu) == vars$myq)
stopifnot(length(vars$logdetDinvfornu) == 1)
stopifnot(length(vars$family_glmm) == 1)
stopifnot(nrow(vars$D.star.inv) == vars$myq)
stopifnot(ncol(vars$D.star.inv) == vars$myq)
stopifnot(length(vars$logdet.D.star.inv) == 1)
stopifnot(length(vars$u.star) == vars$myq)
stopifnot(nrow(umatparams$Sigmuh.inv) == vars$myq)
stopifnot(ncol(umatparams$Sigmuh.inv) == vars$myq)
stopifnot(length(logdet.Sigmuh.inv) == 1)
stopifnot(length(vars$Tee) == 1)
stopifnot(length(vars$nrandom) == vars$Tee)
stopifnot(length(vars$nrandom) == length(vars$nu))
stopifnot(length(vars$meow) == (vars$Tee+1))
stopifnot(length(vars$nu) == vars$Tee)
stopifnot(length(vars$zeta) == 1)
stopifnot(length(vars$tconst) == 1)
stopifnot(length(vars$mod.mcml$ntrials) == vars$n)
stopifnot(length(vars$wts) == vars$n)
.C(C_valgrad,
y = as.double(vars$mod.mcml$y),
Umat = as.double(t(umatparams$umat)),
myq = as.integer(vars$myq),
m = as.integer(umatparams$m),
x = as.double(vars$mod.mcml$x),
n = as.integer(vars$n),
nbeta = as.integer(vars$nbeta),
beta = as.double(vars$beta),
z = as.double(vars$Z),
Dinvfornu = as.double(vars$Dinvfornu),
logdetDinvfornu = as.double(vars$logdetDinvfornu),
family_glmm = as.integer(vars$family_glmm),
Dstarinv = as.double(vars$D.star.inv),
logdetDstarinv = as.double(vars$logdet.D.star.inv),
ustar = as.double(vars$u.star),
Sigmuhinv = as.double(umatparams$Sigmuh.inv),
logdetSigmuhinv = as.double(logdet.Sigmuh.inv),
pee = as.double(vars$pea),
nps = as.integer(length(vars$pea)),
T = as.integer(vars$Tee),
nrandom = as.integer(vars$nrandom),
meow = as.integer(vars$meow),
nu = as.double(vars$nu),
zeta = as.integer(vars$zeta),
tconst = as.double(vars$tconst),
v = double(umatparams$m),
ntrials = as.integer(vars$ntrials),
value = double(1),
gradient = double(length(vars$par)),
b = double(umatparams$m),
wts = as.double(vars$wts))
}
#foreach runs loop in parallel, dopar operator sends each chunk of umat to seperate core and runs the .C function
#combining the b's from each core into one vector
vars$bs <- as.list(rep(0, vars$no_cores))
for(i in 1:vars$no_cores){
vars$bs[[i]] <- out[[i]]$b
}
vars$b <- Reduce(c, vars$bs)
#recalculating the approximation of the log-likelihood value
#finding a, the max of values
vals <- c(rep(0, vars$no_cores))
for(i in 1:vars$no_cores){
vals[i] <- out[[i]]$value
}
a <- max(vals)
#storing normalized b[k], collected from values
normalbk <- c(rep(0, vars$no_cores))
for(i in 1:vars$no_cores){
normalbk[i] <- out[[i]][[4]]*exp(out[[i]]$value - a)
}
sumnbk <- sum(normalbk)
#value <- log(sumnbk/length(out[[1]]$v)) + a
value <- log(sumnbk/vars$newm) + a
#recalculating the gradient
vars$gradient <- c(rep(0, length(out[[1]]$gradient)))
for(j in 1: length(out[[1]]$gradient)){
gradadd <- 0
for(i in 1:length(normalbk)){
gradadd <- gradadd + normalbk[i]*out[[i]]$gradient[j]
}
vars$gradient[j] <- gradadd/sumnbk
}
#parallelizing calculations for the hessian
clusterExport(vars$cl, c("vars"), envir = environment())
out2 <- foreach(i=1:vars$no_cores) %dopar% {
.C(C_hess,
as.double(vars$mod.mcml$y),
as.double(t(umatparams$umat)),
as.integer(vars$myq),
as.integer(umatparams$m),
as.double(vars$mod.mcml$x),
as.integer(vars$n),
as.integer(vars$nbeta),
as.double(vars$beta),
as.double(vars$Z),
as.double(vars$Dinvfornu),
as.double(vars$logdetDinvfornu),
as.integer(vars$family_glmm),
as.double(vars$D.star.inv),
as.double(vars$logdet.D.star.inv),
as.double(vars$u.star),
as.double(umatparams$Sigmuh.inv),
as.double(logdet.Sigmuh.inv),
pea=as.double(vars$pea),
nps=as.integer(length(vars$pea)),
T=as.integer(vars$Tee),
nrandom=as.integer(vars$nrandom),
meow=as.integer(vars$meow),
nu=as.double(vars$nu),
zeta=as.integer(vars$zeta),
tconst=as.double(vars$tconst),
v=double(umatparams$m),
ntrials=as.integer(vars$ntrials),
gradient=as.double(vars$gradient),
hessian=double((length(vars$par))^2),
b=as.double(vars$b),
length=as.integer(vars$newm),
q=as.double(vars$bs[[i]]),
wts=as.double(vars$wts))}
#adding hessian components
hessian <- c(rep(0, length(out2[[1]]$hessian)))
for(j in 1:length(out2[[1]]$hessian)){
hessadd <- 0
for(i in 1:vars$no_cores){
hessadd <- hessadd + out2[[i]]$hessian[j]
}
hessian[j] <- hessadd
}
#making weights accessible
weight <- as.list(rep(0, vars$no_cores))
for(i in 1:vars$no_cores){
weight[[i]] <- out2[[i]]$v
}
weights <- Reduce(c, weight)
if (!missing(cache)) cache$weights<-weights
list(value=value,gradient=vars$gradient,hessian=matrix(hessian,ncol=length(vars$par),byrow=FALSE))
}
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.