##############################
Hook and regal data
R<-c(0,0,0,1,1,1,1)
B<-c(0,1,1,0,0,1,1)
D<-c(1,0,1,0,1,0,1)
y.HR<-c(49,247,142,60,4,112,12)

regstuff<-glm(y.HR~R+B+D+D:R,family=poisson)
modmat<-model.matrix(regstuff)
k<-3
X.list<-make.X.list(k)
X<-make.X.all(k)
beta<-coef(regstuff)
repcount<-2^c(2,2,4,2,4,4,6)
lambda<-as.vector(exp(predict(regstuff)))
piA<-log(0.05/(1-0.05))
stuff<-make.log.p(k,piA)
theta<-c(beta,piA)
Nsim<-1000

result.mat<-matrix(0,length(theta),Nsim)
for(i in 1:Nsim){
print(i)
M<-rpois(5^k-1, rep(lambda,repcount)*exp(unlist(stuff[[1]])))
y<-as.vector(X%*%M)
result.mat[,i]<-solver(y,k, X.list, modmat, theta)$theta
}

###########################################################
#

############################################################
get.var<-function(k, X.list, modmat, theta){
# calculates covariance matrix of thetahat, 
# given value of theta
	mu<-0
        V <- 0
        D1 <- 0
        D2 <- 0        
	r <- length(theta)
        piA <- theta[r]
        beta <- theta[1:(r - 1)]
        stuff <- make.log.p(k, piA)
        log.p <- stuff[[1]]
        d.log.p <- stuff[[2]]
        lambda <- as.vector(exp(modmat %*% beta))
        K <- 2^k - 1
        for(j in (1:K)) {
	 Xj <- X.list[[j]] 
	 	pj <- exp(log.p[[j]])
	 	muj <- as.vector(lambda[j] * Xj %*% pj)
                mu <- mu + muj
                V <- V + lambda[j] * (Xj %*% (pj * t(Xj)))
                D1 <- D1 + outer(muj, modmat[j,  ])
                D2 <- D2 + lambda[j] * Xj %*% (pj * d.log.p[[j]])
                
        }
        D <- cbind(D1, D2)
        A <- t(D) %*% moore.penrose(V)
        solve(A%*%D)
}

  
######################################################


get.var.b<-function(theta,k,X.list,modmat)get.var(k,X.list,modmat,theta)[1,1]
vars<-apply(result.mat,2,get.var.b,k,X.list,modmat)
est.se<-sqrt(vars)*exp(result.mat[1,])


# results

# piA

theta[6] 
-2.944439
> mean(result.mat[6,])
[1] -2.957468

> get.var(k,X.list,modmat,theta)[6,6]
[1] 0.02359298

> var(result.mat[6,])
[1] 0.02475154



# a_0^T\beta
> theta[1]
 (Intercept) 
     4.65328
> mean(result.mat[1,])
[1] 4.648438

> get.var(k,X.list,modmat,theta)[1,1]
[1] 0.017499

> var(result.mat[1,])
[1] 0.01743578



# lambda_0

> exp(theta[1])
 (Intercept) 
    104.9286

> mean(exp(result.mat[1,]))
[1] 105.3287


true.se<-sqrt(get.var.b(theta,k,X.list,modmat))*exp(theta[1])
> true.se
 (Intercept) 
    13.88035

> sqrt(var(exp(result.mat[1,])))
[1] 13.83014


# estimated standard errors

> mean(est.se)
[1] 13.93948

# conf int

> mean(abs(theta[1]-result.mat[1,])<1.96*sqrt(vars))
[1] 0.954