Nothing
library(glmm)
data(BoothHobert)
clust <- makeCluster(2)
set.seed(1234)
mod.mcml1<-glmm(y~0+x1,list(y~0+z1),varcomps.names=c("z1"), data=BoothHobert, family.glmm=bernoulli.glmm, m=21, doPQL=TRUE, debug=TRUE, cluster=clust)
mod.mcml<-mod.mcml1$mod.mcml
z<-mod.mcml$z[[1]]
x<-mod.mcml$x
y<-mod.mcml$y
stuff<-mod.mcml1$debug
beta.pql<-stuff$beta.pql
nu.pql<-stuff$nu.pql
u.pql<-u.star<-stuff$u.star
umat<-stuff$umat
m1<-stuff$m1
family.glmm<-bernoulli.glmm
objfun<-glmm:::objfun
getEk<-glmm:::getEk
addVecs<-glmm:::addVecs
if(is.null(mod.mcml$weights)){
wts <- rep(1, length(y))
} else{
wts <- mod.mcml$weights
}
############################################
#this should be the same as elc
logfyuk<-function(eta,x,y){
value<-sum(y*eta)-sum(log(1+exp(eta)))
Pi<-exp(eta)/(1+exp(eta))
gradient<-sum(y*x)-sum(x*Pi)
hessian<-sum(x^2*(-Pi+Pi^2) )
list(value=value,gradient=gradient,hessian=hessian)
}
#compare elc and logfyuk for a value of eta
eta<-rep(2,150)
ntrials <- rep(1, 150)
this<-.C(glmm:::C_elc,as.double(mod.mcml$y), as.double(mod.mcml$x), as.integer(nrow(mod.mcml$x)), as.integer(ncol(mod.mcml$x)), as.double(eta), as.integer(1), ntrials=as.integer(ntrials), wts=as.double(wts),value=double(1), gradient=double(ncol(mod.mcml$x)), hessian=double((ncol(mod.mcml$x)^2)))
that<-logfyuk(eta,mod.mcml$x,mod.mcml$y)
all.equal(as.numeric(this$value),as.numeric(that$value))
all.equal(as.numeric(this$gradient),as.numeric(that$gradient))
all.equal(as.numeric(this$hessian),as.numeric(that$hessian))
#compare elval to logfyuk
this<-.C(glmm:::C_elval, as.double(mod.mcml$y), as.integer(nrow(mod.mcml$x)), as.integer(ncol(mod.mcml$x)), as.double(eta), as.integer(1), ntrials=as.integer(ntrials), wts=as.double(wts), value=double(1))
all.equal(as.numeric(this$value),as.numeric(that$value))
#compare elGH to logfyuk
this<-.C(glmm:::C_elGH, as.double(mod.mcml$y), as.double(mod.mcml$x), as.integer(nrow(mod.mcml$x)), as.integer(ncol(mod.mcml$x)), as.double(eta), as.integer(1), ntrials=as.integer(ntrials), wts=as.double(wts), gradient=double(ncol(mod.mcml$x)), hessian=double((ncol(mod.mcml$x)^2)))
all.equal(as.numeric(this$gradient),as.numeric(that$gradient))
all.equal(as.numeric(this$hessian),as.numeric(that$hessian))
############################################
#want to check distRand when we use a normal distribution to get our random effects
#check written for BH
distRandCheck<-function(nu,uvec,muvec){
ukmuk<-sum((uvec-muvec)^2)
value<- -length(uvec)*.5*log(2*pi)-5*log(nu)-ukmuk/(2*nu)
gradient<- -5/nu +ukmuk/(2*nu^2)
hessian<- 5/(nu^2)-ukmuk/(nu^3)
hessian<-as.matrix(hessian)
list(value=value,gradient=gradient,hessian=hessian)
}
#function written originally in R
distRand <-
function(nu,U,z.list,mu){
# T=number variance components
T<-length(z.list)
#nrandom is q_t
nrand<-lapply(z.list,ncol)
nrandom<-unlist(nrand)
totnrandom<-sum(nrandom)
mu.list<-U.list<-NULL
if(T==1) {
U.list[[1]]<-U
mu.list[[1]]<-mu
}
if(T>1){
U.list[[1]]<-U[1:nrandom[1]]
mu.list[[1]]<-mu[1:nrandom[1]]
for(t in 2:T){
thing1<-sum(nrandom[1:t-1])+1
thing2<-sum(nrandom[1:t])
U.list[[t]]<-U[thing1:thing2]
mu.list[[t]]<-mu[thing1:thing2]
}
}
val<-gradient<-Hessian<-rep(0,T)
#for each variance component
for(t in 1:T){
you<-as.vector(U.list[[t]])
mew<-as.vector(mu.list[[t]])
Umu<-(you-mew)%*%(you-mew)
val[t]<--length(you)*.5*log(2*pi)+ as.numeric(-.5*nrandom[t]*log(nu[t])-Umu/(2*nu[t]))
gradient[t]<- -nrandom[t]/(2*nu[t])+Umu/(2*(nu[t])^2)
Hessian[t]<- nrandom[t]/(2*(nu[t])^2)- Umu/((nu[t])^3)
}
value<-sum(val)
if(T>1) hessian<-diag(Hessian)
if(T==1) hessian<-matrix(Hessian,nrow=1,ncol=1)
list(value=value,gradient=gradient,hessian=hessian)
}
you<-umat[1,]
this<-distRandCheck(2,you,u.pql)
that<-distRand(2,you,mod.mcml$z,u.pql)
all.equal(this,that)
#use finite diffs to make sure distRandCheck (and distRand) have correct derivs
del<-10^(-9)
thisdel<-distRandCheck(2+del,you,u.pql)
firstthing<-thisdel$value-this$value
secondthing<-as.vector(this$gradient%*%del)
all.equal(firstthing,secondthing)
#firstthing
#secondthing
#firstthing-secondthing
#compare the gradient and hessian of the C functions by using these functions
#(the value is checked in distRandGeneral)
mynu<-2
mymu<-rep(0,10)
T<-1
nrandom<-10
meow<-c(0,10)
set.seed(1234)
myyou<-rnorm(10)
hohum<-.C(glmm:::C_distRand3C,as.double(mynu), as.double(mymu), as.integer(T), as.integer(nrandom), as.integer(meow), as.double(myyou), double(T), double(T^2))
drcheck<-distRandCheck(mynu,myyou,mymu)
all.equal(drcheck$gradient,hohum[[7]])
all.equal(drcheck$hessian,matrix(hohum[[8]],nrow=T,byrow=F))
###############################################
#distRandGeneral in R first
distRandGeneral<-function(uvec,mu,Sigma.inv){
logDetSigmaInv<-sum(log(eigen(Sigma.inv,symmetric=TRUE)$values))
umu<-uvec-mu
piece2<-t(umu)%*%Sigma.inv%*%umu
out<-as.vector(.5*(logDetSigmaInv-piece2))
const<-length(uvec)*.5*log(2*pi)
out<-out-const
out
}
D.star<-2*diag(10)
D.star.inv<-.5*diag(10)
A.star<-sqrt(2)*diag(10)
this<-distRandGeneral(you,u.pql,D.star.inv)
all.equal(this,that$value)
#check distRandGenC
logdet<-sum(log(eigen(D.star.inv)$values))
stuff<-.C(glmm:::C_distRandGenC,as.double(D.star.inv),as.double(logdet), as.integer(length(you)), as.double(you), as.double(u.pql), double(1))[[6]]
all.equal(that$value,stuff)
############################################
#want to check that the value of objfun is the same for a value of nu and beta
vars <- new.env(parent = emptyenv())
debug<-mod.mcml1$debug
vars$m1 <- debug$m1
m2 <- debug$m2
m3 <- debug$m3
vars$zeta <- 5
vars$cl <- mod.mcml1$cluster
registerDoParallel(vars$cl) #making cluster usable with foreach
vars$no_cores <- length(vars$cl)
vars$mod.mcml<-mod.mcml1$mod.mcml
vars$nu.pql <- debug$nu.pql
vars$umat<-debug$umat
vars$newm <- nrow(vars$umat)
vars$u.star<-debug$u.star
D <- vars$D.star <- Dstarnotsparse<-2*diag(10)
D.inv <- D.star.inv <-.5*diag(10)
getEk<-glmm:::getEk
addVecs<-glmm:::addVecs
genRand<-glmm:::genRand
vars$family.glmm<-mod.mcml1$family.glmm
vars$ntrials<- rep(1, length(mod.mcml1$y) )
beta.pql <- debug$beta.pql
vars$wts<-wts
length(wts) == length(vars$ntrials)
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)
Sigmuh.inv<- t(Z)%*%cdouble%*%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
objfun<-glmm:::objfun
umats <- clusterEvalQ(vars$cl, umatparams$umat)
umat <- Reduce(rbind, umats)
Sigmuh.invs <- clusterEvalQ(vars$cl, umatparams$Sigmuh.inv)
Sigmuh.inv <- Sigmuh.invs[[1]]
Sigmuh <- solve(Sigmuh.inv)
dbb<-db<-b<-rep(0,vars$newm)
sigsq<-nu<-2
beta<-6
Z<-vars$mod.mcml$z[[1]]
A<-sqrt(2)*diag(10)
eta.star<-x*beta.pql+as.vector(Z%*%u.star)
cdouble<-as.vector(bernoulli.glmm()$cpp(eta.star)) #still a vector
cdouble<-diag(cdouble)
piece3<-rep(0,3)
#calculate objfun's value for comparison
cache<-new.env(parent = emptyenv())
that<-objfun(c(beta,nu), cache=cache,vars=vars)
#get t stuff ready
tconstant<-glmm:::tconstant
zeta<-5
tconst<-tconstant(zeta,10,diag(D.star.inv))
tdist2<-function(tconst,u, Dstarinv,zeta,myq){
inside<-1+t(u)%*%Dstarinv%*%u/zeta
logft<-tconst - ((zeta+myq)/2)*log(inside)
as.vector(logft)
}
#now go through row by row of umat
#ie go through each vector of gen rand eff
for(k in 1:vars$newm){
uvec<-umat[k,]
eta<-x*beta+as.vector(Z%*%uvec)
piece1<- logfyuk(eta,x,y)$value
piece2<- distRandCheck(nu,uvec,rep(0,10))$value
piece3[1]<-tdist2(tconst,uvec,D.star.inv,zeta,10)
piece3[2]<- distRandGeneral(uvec, vars$u.star, D.star.inv)
piece3[3]<-distRandGeneral(uvec,vars$u.star,Sigmuh.inv)
damax<-max(piece3)
blah<-sum(exp(piece3-damax)/3)
lefoo<-damax+log(blah)
b[k]<-piece1+piece2-lefoo
}
a<-max(b)
top<-exp(b-a)
value<-a+log(mean(top))
all.equal(value,that$value)
#Given generated random effects, the value of the objective function is correct.
#This plus the test of finite diffs for objfun should be enough.
stopCluster(clust)
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.