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.