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 [EMAIL PROTECTED]
[[alternative HTML version deleted]]
______________________________________________
[email protected] 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.
