95% probability intervals (PI)
I typicallly work with models and MCMC algorithms with many unknown parameters. At the end of the algorithm, all I really want is the mean and 95% PI. Here is some code that only saves the extreme 5% (2.5%,97.5%) of the posterior draws for a parameter. This cuts down on memory space if there are a lot of unknown parameters.
find95=function(y,List,q,Q) #y is the new draw or updated parameter
{
v=round(.025*Q,0) #length of 2.5 percent
if(q < 2*v) #first 5 percent are used as intializers
{ init=cbind(List$init,y)
return(list("init"=init))
}
if(q==2*v)
{ List$init=cbind(List$init,y)
lb =sort(List$init)[1:v]
ub =sort(List$init)[(v+1):(2*v)]
}
if(q > 2*v)
{ if(y < max(List$lb)){ lb=sort(c(y,List$lb))[1:v]}else{ lb=List$lb }
if(y > min(List$ub)){ ub=sort(c(y,List$ub))[2:(v+1)]}else{ ub=List$ub }
}
return(list("lb"=lb, "ub"=ub))
}
###### Example
List=list("init"=numeric(0))
Q=1001 #number of iterations
for(q in 1:Q)
{ y=rnorm(1,0,1)
List=find95(y, List, q,Q)
}
List$lb[v] #2.5 percent lower bound
List$ub[1] #97.5 percent upper bound