# pfMLLik1: Create a function for computing the log of an unbiased... In smfsb: Stochastic Modelling for Systems Biology

## Description

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 "log-sum-exp trick" for avoiding numerical underflow. See `pfMLLik` for a version which does.

## Usage

 `1` ```pfMLLik1(n,simx0,t0,stepFun,dataLik,data) ```

## Arguments

 `n` An integer representing the number of particles to use in the particle filter. `simx0` A function with interface `simx0(n,t0,...)`, where `n` is the number of rows of a matrix and `t0` is a time at which to simulate from an initial distribution for the state of the particle filter. The return value should be a matrix whose rows are random samples from this distribution. The function therefore represents a prior distribution on the initial state of the Markov process. `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 `StepGillespie`. `dataLik` A function with interface `dataLik(x,t,y,log=TRUE,...)`, where `x` and `t` represent the true state and time of the process, and `y` is the observed data. The return value should be the (log of the) likelihood of the observation. The function therefore represents the observation model. `data` A timed data matrix representing the observations, such as produced by `simTimes` or `as.timedData`.

## Value

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))) ```