Description Usage Arguments Value See Also Examples
Create a function for computing the log of an unbiased estimate of marginal likelihood of a time course data set using a simple bootstrap particle filter. This version does not use the "logsumexp trick" for avoiding numerical underflow.
See pfMLLik
for a version which does.
1 
n 
An integer representing the number of particles to use in the particle filter. 
simx0 
A function with interface 
t0 
The time corresponding to the starting point of the Markov process. Can be no bigger than the smallest observation time. 
stepFun 
A function for advancing the state of the Markov process, such as returned by 
dataLik 
A function with interface

data 
A timed data matrix representing the observations, such as produced by 
An R function with interface (...)
which evaluates to the log of the particle filter's unbiased estimate of the marginal likelihood of the data.
pfMLLik
, StepGillespie
, as.timedData
,
simTimes
, stepLVc
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25  noiseSD=5
# first simulate some data
truth=simTs(c(x1=50,x2=100),0,20,2,stepLVc)
data=truth+rnorm(prod(dim(truth)),0,noiseSD)
data=as.timedData(data)
# measurement error model
dataLik < function(x,t,y,log=TRUE,...)
{
ll=sum(dnorm(y,x,noiseSD,log=TRUE))
if (log)
return(ll)
else
return(exp(ll))
}
# now define a sampler for the prior on the initial state
simx0 < function(N,t0,...)
{
mat=cbind(rpois(N,50),rpois(N,100))
colnames(mat)=c("x1","x2")
mat
}
mLLik=pfMLLik1(1000,simx0,0,stepLVc,dataLik,data)
print(mLLik())
print(mLLik(th=c(th1 = 1, th2 = 0.005, th3 = 0.6)))
print(mLLik(th=c(th1 = 1, th2 = 0.005, th3 = 0.5)))

Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.