[R] Preserving output of MCMC iterations

[Ben Bolker]

Francesco Checchi wrote:
> 
> Dear colleagues,
>  
> I'm an epidemiologist with no background in programming and just
> started using R a few weeks ago. I am working on the epidemiology of
> African sleeping sickness, and am trying to use R to perform a Monte
> Carlo Markov Chain analysis to estimate three unknown parameters within
> a model of African sleeping sickness case detection that is mainly
> informed by actual field programme data.
> Basically, I have a set of four differential equations governing my
> model, which I am implementing stochastically (small populations, low
> transmission). I iterate the model a large number of times (10 000) for
> each unique combination of plausible values of the three unknown
> parameters a, b and c. After each iteration, I ask R to check whether
> the output 'fits' the observed data: if it does, I ask R to add 1 to the
> (i,j,k)th element of a 3-dimensional 'scores' array (3 unknown
> parameters = 3 dimensions), where i,j, and k are the plausible values of
> a,b and c that are currently being tested.
> My intended output is in fact this 3-dimensional array, which is
> essentially a likelihood surface wherein each (i,j,k)th element is the
> likelihood of the corresponding combination of plausible parameter
> values a, b and c.
>  
> The code I have written looks more or less like this:
>  
> FUN1(args) {
> specify parameters based on arguments
>  
> replicate (10 000 times, {
>   {implement differential equations once}
>   if output fits data, scores[i,j,k]<-scores[i,j,k]+1
>   })
> }
>  
> FUN2(args) {
> create 3-dim 'scores' array
> apply FUN1 to each combination of a,b and c values, i.e. to each row of
> a dataframe containing each combination
> return(scores)
> }
>  
> I've tested the various sub-routines within the code repeatedly and
> everything seems to work fine with the exception of one problem: the
> scores array that R returns when I call FUN2 is unchanged from its
> initial specification within FUN2 (e.g. if I tell R the initial value of
> the array elements should be 0, R returns an array with all 0s). The
> scores array IS being updated with +1 after each data-fitting iteration
> of the model (I checked this), but somehow these changes are not being
> preserved outside of the replicate {} environment.
>  
> I've tried to look at the help archives and other help sources but
> haven't been able to find or perhaps interpret solutions to my problem.
> I've also looked at the contributed MCMC packages already available, but
> for a variety of reasons these do not seem suitable for my purposes, and
> I would prefer writing the full code myself.
>  
> Apologies if the explanation is not very clear - I'm still in the very
> steep learning curve! I will be happy to clarify further.
>  
> Thanks in advance for your consideration.
>  
> Francesco
>  
>  
>  
>  
>  
>  
> Francesco Checchi
> Honorary Lecturer, Part-time PhD Student
> Disease Control and Vector Biology Unit
> Department of Infectious and Tropical Diseases
> London School of Hygiene and Tropical Medicine
> Room 51G6, 49-51 Bedford Square
> London WC1B 3DP, United Kingdom
> office +44 (0)20 7927 2336
> mobile +44 (0)79 0671 9637
> fax +44 (0)20 7927 2918
> e-mail  francesco.checchi at lshtm.ac.uk
> 
> 	[[alternative HTML version deleted]]
> 
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
> 
> 

The quickest solution to your problem will probably be to use
the dreaded "global assignment" operator, <<-  (double-arrow)
to alter the global copy of "scores".  (What's happening is that
a copy of scores is getting created in the environment of
FUN1, and then discarded when FUN1 completes.)

 The more R-idiomatic way of doing this would be:


FUN1(args) {
specify parameters based on arguments
  score <- 0
replicate (10 000 times, {
  {implement differential equations once}
  if output fits data, score <-score+1
  })
  score  ## or return(score)
}
 
FUN2(args) {
create 3-dim 'scores' array
apply FUN1 to each combination of a,b and c values, i.e. to each row of
a dataframe containing each combination
return(scores)
   ## what do you mean by apply?  
  ## literal: parameters stored in a 3-dimensional array ...
    scores <- apply(param_array,c(1,2),FUN1)
  ## not so literal:  

  create scores array

   param1vec = seq(0,5,by=0.1)
  param2vec =  seq(0,5,by=0.1)
  param3vec = seq(0,5,by=0.1)
  for (i in seq(along=param1vec)) {
     for (j in seq(along=param2vec)) {
       for (k in seq(along=param3vec)) {
           scores[i,j,k]  <- FUN1(param1vec[i],param2vec[j],param3vec[k])
      }
   }
}
 
  ... or something like that.


-- 
View this message in context: http://www.nabble.com/Preserving-output-of-MCMC-iterations-tf4877082.html#a13960804
Sent from the R help mailing list archive at Nabble.com.