# Implementing Likelihoods for quadratic models

Working from the quadratic model for the saddle node bifurcation rather than the linear (OU) model. Recall that the dynamics are:

$dX = \left( r - (x-\theta)^2\right)dt + \beta \sqrt{X} dB_t$

and the linear noise approximation gives a linear Fokker Planck equation such that the probability distribution is Gaussian with moments ((corrected from yesterday to scale noise, basically as pure-birth Poisson process))

$\frac{d \hat{x}}{dt} &= r - (\hat{x}-\theta)^2$ $\frac{d \sigma^2}{dt} &= -2 (\hat x - \theta )\sigma^2 + \beta \left(r - (\hat{x}-\theta)^2 \right)$

So we can evaluate the probability density of being at $x_t$ at time t given that we started at x_0 at time 0 by integrating the above equations (numerically, in this case, as we lack a closed form solution for the above).

Being able to calculate this probability we can easily calculate likelihood of the timeseries under the model. We code this in R (likelihood_bifur_models.R) and write a quick example script (also in repo as bifur_likelihood_ex.R):


pars = c(r=10, theta=3, beta=1)
m <- init_sdemodel(pars =pars, Xo = 6.2, model="SN", N=200)
X <- simulate.SN(m)
# Set initial guess
m\$pars <- c(r=11,theta=4,beta=1)
out <- update.SN(m, X, method="SANN")

Exploring different fitting, seems to be a challenge. For instance, finds wrong optimum sometimes – perhaps due to a bug somewhere (see more examples in flickr stream):

As likelihood surface (shown here for theta) appears challenging:

and numerical errors make it worse:

Never the less, the fit does okay when the initial condition for theta starts nearer the better well (right). See flickr library for more examples and git for code. Still not obvious why the unstable equilibrium is quite such a likelihood peak, but anyway…

Success of fits varies quite a bit, depending on data, initial guesses for parameters, etc. Note that for many of these the optimization fails to converge (more reasons discussed below).


system(paste('hpc-autotweets "', gitcom, ' done"'))