User:Carl Boettiger/Notebook/Stochastic Population Dynamics/2010/05/06
From OpenWetWare
Stochastic Population Dynamics  Main project page Previous entry Next entry 
Goals
Test CasedX_{t} = α_{t}(θ − X_{t})dt + σdt dα_{t} = − β
Approach
Conditional probabilitySolution to the time dependent OU process, see Gardiner 4th ed. pg 111, eq. (4.5.109). Solve using and apply Ito formula, giving And the moments are:
numeric_V < function(Dt, pars){ vint < function(x){ exp(pars$beta*x^2 2* pars$alpha_0*x) } int < integrate(vint, 0, Dt) int$value*pars$sigma^2 } analytic_V < function(Dt, pars){ if( erf(pars$alpha/sqrt(pars$beta) ) != 1){ Log_Vx < log(pars$sigma^2 * sqrt(pi) ) + pars$alpha_0^2/(pars$beta) + log( ( erf( (pars$alpha_0+pars$beta*Dt)/(sqrt(pars$beta) ) )  erf(pars$alpha_0/(sqrt(pars$beta) ) ))/(2*sqrt(pars$beta) )) class(Log_Vx) = "numeric" Vx < exp(Log_Vx) } else { warning("beta near zero, using approximation") Vx < pars$sigma^2*(1exp(2*pars$alpha*Dt))/(2*pars$alpha) } Vx } # Parametrization dXt = alpha(theta  Xt)dt + sigma*dWt, alpha(t) = beta*t+alpha_0 warning_model < function(Dt, Xo, pars){ int < pars$beta*Dt^2/2 + pars$alpha_0*Dt Ex < Xo * exp(int) + pars$theta * (1  exp(int) ) Vx < analytic_V(Dt, pars) return(list(Ex=Ex,Vx=Vx)) } dcWarning < function(x, Dt, x0, pars, log = FALSE){ P < warning_model(Dt, x0, pars) dnorm(x, mean=P$Ex, sd=sqrt(P$Vx), log=log) } warning.lik < function(alpha_0, theta, sigma, beta){ n < length(X) dt < deltat(X) pars = list(alpha_0=alpha_0, theta=theta, sigma=sigma, beta=beta) sum( dcWarning(X[2:n], dt, X[1:(n1)], pars, log=TRUE) ) } require(stats4) require(sde) require(NORMT3) set.seed(123) theta < 3 alpha < 1 sigma < 2 alpha_0 < 1 beta < .01 pars = list(alpha_0=alpha, theta=theta, sigma=sigma, beta=beta) warning_model(1, 1, pars) setOU(1, 1, theta=c(theta*alpha, alpha, sigma) ) numeric_V(Dt, pars) analytic_V(Dt, pars)
Note on Numerics
Misc Notes
#!/bin/bash xrandr newmode "1920x1200_60.00" 193.16 1920 2048 2256 2592 1200 1201 1204 1242 HSync +Vsync xrandr newmode "1280x1040_60.00" 110.63 1280 1360 1496 1712 1040 1041 1044 1077 HSync +Vsync xrandr addmode VGA1 1920x1200_60.00 xrandr addmode VGA1 1280x1040_60.00
