https://github.com/cran/HLSM
Tip revision: 31ec11bc3c0ac504839d98ce0e6e8feb8951fbbe authored by Samrachana Adhikari on 18 June 2014, 00:00:00 UTC
version 0.1
version 0.1
Tip revision: 31ec11b
HLSM_run.R
####Function to run the sampler#####
##SPECIFYING PRIORS######
#########################
##if priors = NULL => uses randomly generated priors
##else: priors is a list of following objects:
##MuBeta:= prior mean for betas & intercepts;
##SigmaBeta:= prior variance for betas & intercept;
##MuAlpha:= prior mean for alpha;
##SigmaAlpha:= prior variance for alpha;
##MuZ; VarZ;
##PriorA; PriorB
##TUNING PARAMETERS#######
##########################
##if tune = NULL => uses auto tuning
##else: tune is a list of following objects:
##tuneAlpha = 0.9
##tuneBeta = array, dim=c(PP,KK)
##tuneInt = vec, len = KK
##tuneZ = list( vec(len = nn[x]])) length of list = KK
##############################################################
#library(MASS)
HLSM=function(X, Y, initialVals = NULL, priors = NULL, tune = NULL,
tuneIn = TRUE, TT = NULL,dd, niter,intervention)
{
#X and Y are provided as list.
nn = sapply(1:length(X),function(x) nrow(X[[x]]))
KK = length(X)
PP = dim(X[[1]])[3]
XX = unlist(X)
YY = unlist(Y)
YY[which(is.na(YY))] = 0
XX[which(is.na(XX))] = 0
#Priors
if(is.null(priors)){
MuBeta= rep(0,(PP+1))
SigmaBeta = rep(1,(PP+1))
MuAlpha=0
SigmaAlpha = 100
MuZ = c(0,0)
VarZ = c(100,100)
PriorA = 4.5
PriorB = 17.5
}else{
if(class(priors) != 'list')(stop("priors must be of class list, if not NULL"))
MuBeta = priors$MuBeta
SigmaBeta = priors$SigmaBeta
MuAlpha = priors$MuAlpha
SigmaAlpha = priors$SigmaAlpha
MuZ = priors$MuZ
VarZ = priors$VarZ
PriorA = priors$PriorA
PriorB = priors$PriorB
}
##starting values
if(is.null(initialVals)){
Z0 = array(0, dim = c(sum(nn),dd))
cc = 1
for(i in 1:KK){
cc1 = (cc-1)+nn[i]
ZZ = t(replicate(nn[i],rnorm(dd,0,1)))
ZZ[1,]=c(0,0)
ZZ[2,2]=0
if(ZZ[2,1] < ZZ[1,1]){
ZZ[2,1] = -1*ZZ[2,1]}
ZZ[3,2] = abs(ZZ[3,2])
Z0[cc:cc1,] = ZZ
cc = cc+nn[i]
}
Z0 = unlist(Z0)
beta0 = replicate(KK,rnorm(PP,0,1))
intercept0 = rnorm(KK, 0,1)
if(intervention == 1){ alpha0=rnorm(1, 0, 1) }
print("Starting Values Set")
}else{
if(class(initialVals)!= 'list')(stop("initialVals must be of class list, if not NULL"))
Z0 = initialVals$ZZ
beta0 = initialVals$beta
intercept0 = initialVals$intercept
if(intervention == 1){ alpha0 = initialVals$alpha}
}
if(intervention == 0){
alpha0 = 0
TT = rep(0, KK)
}
###tuning parameters#####
if(is.null(tune)){
a.number = 5
tuneAlpha = 0.9
tuneBeta = array(1,dim=c(PP,KK))
tuneInt = rep(0.2,KK)
tuneZ = lapply(1:KK,function(x) rep(1.2,nn[x]))
} else{
if(class(tune) != 'list')(stop("tune must be of class list, if not NULL"))
a.number = 1
tuneAlpha = tune$tuneAlpha
tuneBeta = tune$tuneBeta
tuneInt = tune$tuneInt
tuneZ = tune$tuneZ
}
###Tuning the Sampler####
do.again = 1
tuneX = 1
if(tuneIn == TRUE){
while(do.again ==1){
print('Tuning the Sampler')
for(counter in 1:a.number){
rslt = MCMCfunction(nn=nn,PP=PP,KK=KK,dd=dd,XX = XX,YY = YY,ZZ = Z0,TT = TT,beta = beta0 ,intercept = intercept0,
alpha = alpha0,MuAlpha = MuAlpha,SigmaAlpha = SigmaAlpha,MuBeta = MuBeta,SigmaBeta = SigmaBeta,MuZ = MuZ,
VarZ = VarZ,tuneBetaAll = tuneBeta, tuneInt = tuneInt, tuneAlpha = tuneAlpha,tuneZAll = unlist(tuneZ),
niter = 200,PriorA = PriorA, PriorB = PriorB, intervention = intervention)
tuneAlpha = adjust.my.tune(tuneAlpha, rslt$acc$alpha,1)
tuneZ = lapply(1:KK,function(x)adjust.my.tune(tuneZ[[x]], rslt$acc$Z[[x]], 2))
tuneBeta = array(sapply(1:KK,function(x)adjust.my.tune(tuneBeta[,x], rslt$acc$beta[,x],1)),dim = c(PP,KK))
tuneInt = sapply(1:KK,function(x)adjust.my.tune(tuneInt[x],rslt$acc$intercept[x], 1))
print(paste('TuneDone = ',tuneX))
tuneX = tuneX+1
}
extreme = lapply(1:KK,function(x)which.suck(rslt$acc$Z[[x]],2))
do.again = max(sapply(extreme, length)) > 5
}
print("Tuning is finished")
}
rslt = MCMCfunction(nn=nn,PP=PP,KK=KK,dd= dd,XX = XX,YY = YY,ZZ = Z0,TT = TT,beta = beta0 ,intercept = intercept0,
alpha = alpha0,MuAlpha = MuAlpha,SigmaAlpha = SigmaAlpha,MuBeta = MuBeta,SigmaBeta = SigmaBeta,MuZ = MuZ,
VarZ = VarZ,tuneBetaAll = tuneBeta, tuneInt = tuneInt, tuneAlpha = tuneAlpha,tuneZAll = unlist(tuneZ),
niter = niter,PriorA = PriorA, PriorB = PriorB,intervention = intervention)
rslt$call = match.call()
class(rslt) = 'HLSM'
rslt
}