Error message

  • Deprecated function: The each() function is deprecated. This message will be suppressed on further calls in book_prev() (line 775 of /home3/gardeoi3/public_html/iamrandom/modules/book/book.module).
  • Deprecated function: implode(): Passing glue string after array is deprecated. Swap the parameters in drupal_get_feeds() (line 394 of /home3/gardeoi3/public_html/iamrandom/includes/common.inc).

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

 

Tags: