Nothing
library(glmm)
data(BoothHobert)
clust <- makeCluster(2)
set.seed(1234)
out<-glmm(y~0+x1,list(y~0+z1),varcomps.names=c("z1"),data=BoothHobert,
family.glmm=bernoulli.glmm,m=50,doPQL=FALSE,debug=TRUE, cluster=clust)
vars <- new.env(parent = emptyenv())
debug<-out$debug
vars$m1 <- debug$m1
m2 <- debug$m2
m3 <- debug$m3
vars$zeta <- 5
vars$cl <- clust
registerDoParallel(vars$cl) #making cluster usable with foreach
vars$no_cores <- length(vars$cl)
vars$umat<-debug$umat
vars$newm <- nrow(vars$umat)
vars$u.star<-debug$u.star
vars$mod.mcml<-out$mod.mcml
vars$nu.pql <- debug$nu.pql
D.star.inv <- Dstarnotsparse <- vars$D.star <- as.matrix(debug$D.star)
getEk<-glmm:::getEk
addVecs<-glmm:::addVecs
genRand<-glmm:::genRand
vars$family.glmm<-out$family.glmm
vars$ntrials<- rep(1, length(out$y))
beta.pql <- debug$beta.pql
if(is.null(out$weights)){
wts <- rep(1, length(out$y))
} else{
wts <- out$weights
}
vars$wts<-wts
simulate <- function(vars, Dstarnotsparse, m2, m3, beta.pql, D.star.inv){
#generate m1 from t(0,D*)
if(vars$m1>0) genData<-rmvt(ceiling(vars$m1/vars$no_cores),sigma=Dstarnotsparse,df=vars$zeta,type=c("shifted"))
if(vars$m1==0) genData<-NULL
#generate m2 from N(u*,D*)
if(m2>0) genData2<-genRand(vars$u.star,vars$D.star,ceiling(m2/vars$no_cores))
if(m2==0) genData2<-NULL
#generate m3 from N(u*,(Z'c''(Xbeta*+zu*)Z+D*^{-1})^-1)
if(m3>0){
Z=do.call(cbind,vars$mod.mcml$z)
eta.star<-as.vector(vars$mod.mcml$x%*%beta.pql+Z%*%vars$u.star)
if(vars$family.glmm$family.glmm=="bernoulli.glmm") {cdouble<-vars$family.glmm$cpp(eta.star)}
if(vars$family.glmm$family.glmm=="poisson.glmm"){cdouble<-vars$family.glmm$cpp(eta.star)}
if(vars$family.glmm$family.glmm=="binomial.glmm"){cdouble<-vars$family.glmm$cpp(eta.star, vars$ntrials)}
#still a vector
cdouble<-Diagonal(length(cdouble),cdouble)
wtsmat <- diag(vars$wts)
Sigmuh.inv<- t(Z)%*%cdouble%*%wtsmat%*%Z+D.star.inv
Sigmuh<-solve(Sigmuh.inv)
genData3<-genRand(vars$u.star,Sigmuh,ceiling(m3/vars$no_cores))
}
if(m3==0) genData3<-NULL
# #these are from distribution based on data
# if(distrib=="tee")genData<-genRand(sigma.gen,s.pql,mod.mcml$z,m1,distrib="tee",gamm)
# if(distrib=="normal")genData<-genRand(sigma.pql,s.pql,mod.mcml$z,m1,distrib="normal",gamm)
# #these are from standard normal
# ones<-rep(1,length(sigma.pql))
# zeros<-rep(0,length(s.pql))
# genData2<-genRand(ones,zeros,mod.mcml$z,m2,distrib="normal",gamm)
umat<-rbind(genData,genData2,genData3)
m <- nrow(umat)
list(umat=umat, m=m, Sigmuh.inv=Sigmuh.inv)
}
clusterSetRNGStream(vars$cl, 1234)
clusterExport(vars$cl, c("vars", "Dstarnotsparse", "m2", "m3", "beta.pql", "D.star.inv", "simulate", "genRand"), envir = environment()) #installing variables on each core
noprint <- clusterEvalQ(vars$cl, umatparams <- simulate(vars=vars, Dstarnotsparse=Dstarnotsparse, m2=m2, m3=m3, beta.pql=beta.pql, D.star.inv=D.star.inv))
vars$nbeta <- 1
vars$p1=vars$p2=vars$p3=1/3
par<-c(6,1.5)
del<-rep(10^-8,2)
objfun<-glmm:::objfun
core2<-objfun(par=par, vars=vars)
umats <- clusterEvalQ(vars$cl, umatparams$umat)
umat <- Reduce(rbind, umats)
Sigmuh.invs <- clusterEvalQ(vars$cl, umatparams$Sigmuh.inv)
Sigmuh.inv <- Sigmuh.invs[[1]]
stopCluster(clust)
vars$cl <- makeCluster(1)
registerDoParallel(vars$cl) #making cluster usable with foreach
vars$no_cores <- length(vars$cl)
clusterExport(vars$cl, c("umat", "vars", "Sigmuh.inv"), envir = environment())
noprint <- clusterEvalQ(vars$cl, umatparams <- list(umat=umat, m=vars$newm, Sigmuh.inv=Sigmuh.inv))
core1<-objfun(par=par, vars=vars)
all.equal(core1, core2)
stopCluster(vars$cl)
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.