Particle Filter
Particle Filter
The particle filter can be used for non-linear and non-Gaussian state space models, such as DGLMs. A general DGLM is given by:
\[\begin{align*} Y_t &\sim p(y_t | g(F_t^Tx_t), \theta) \\ X_t &= G_tX_{t-1} + w_t, \qquad w_t \sim \mathcal{N}(0, W), \\ X_0 &\sim \mathcal{N}(m_0, C_0). \end{align*}\]
Where \(p(y_t | F_t^Tx_t, \theta)\) is an Exponential Family distribution, with optional static parameter \(\theta\).
The particle filter can be used to find the filtering state \(p(x_{0:t}|y_{1:T}, \psi)\), where \(\psi\) represents the static parameters in the model. The particle filter algorithm proceeds as follows:
Initialisation 1. Sample \(k\) initial particles from the initial distribution, \(x_0^{(i)} \sim \mathcal{N}(m_0, C_0), i = 1,\dots,k\) 2. Initialise the particle weights \(w^{(i)} = \frac{1}{k}\)
Main Loop
- Assume at time \(t\) we have weighted sample from the posterior \(\{x_t^{(i)}, w_t^{(i)}, i = 1,\dots,k\}\)
- Advance the state using the state evolution kernel \(x_{t+1}^{(i)} \sim \mathcal{N}(G_tx_t^{(i)}, W), i = 1,\dots,k\)
- Calculate the weights using the conditional likelihood: \(w_{t+1}^{(i)} = p(y_{t+1}|x_{t+1}^{(i)}), i = 1,\dots,k\)
- Normalise the weights, \(\pi_{t+1}^{(i)} = \frac{w_{t+1}^{(i)}}{\sum_{j=1}^kw_{t+1}^{(j)}}, i = 1,\dots,k\)
- Resample the particles by sampling \(u_j \sim p(\pmb\pi_{t+1})\) then select \(x^{(u_j)}_{t+1}\) \(j = 1,\dots,k\)
Repeat the mean loop until \(t = T\).
To apply the particle filter to a model using this library, specify either a DGLM or a DLM model. Consider a poisson model
scala> import com.github.jonnylaw.dlm._
<console>:12: warning: Unused import
import com.github.jonnylaw.dlm._
^
import com.github.jonnylaw.dlm._
scala> import breeze.linalg.{DenseMatrix, DenseVector}
<console>:10: warning: Unused import
import com.github.jonnylaw.dlm._
^
<console>:15: warning: Unused import
import breeze.linalg.{DenseMatrix, DenseVector}
^
<console>:15: warning: Unused import
import breeze.linalg.{DenseMatrix, DenseVector}
^
import breeze.linalg.{DenseMatrix, DenseVector}
scala> val mod = Dglm.poisson(Dlm.polynomial(1))
<console>:12: warning: Unused import
import breeze.linalg.{DenseMatrix, DenseVector}
^
<console>:12: warning: Unused import
import breeze.linalg.{DenseMatrix, DenseVector}
^
mod: com.github.jonnylaw.dlm.Dglm = Dglm(com.github.jonnylaw.dlm.Dglm$$$Lambda$7171/2056819218@77b45cca,com.github.jonnylaw.dlm.Dlm$$$Lambda$6788/1540206298@362e0e9d,com.github.jonnylaw.dlm.Dlm$$$Lambda$6789/1989125347@68256bc7,com.github.jonnylaw.dlm.Dglm$$$Lambda$7172/1171156020@4982582c,com.github.jonnylaw.dlm.Dglm$$$Lambda$7173/490903945@7f9b32e4)
scala> val params = DlmParameters(DenseMatrix(2.0),
| DenseMatrix(0.05),
| DenseVector(0.0),
| DenseMatrix(1.0))
params: com.github.jonnylaw.dlm.DlmParameters = DlmParameters(2.0 ,0.05 ,DenseVector(0.0),1.0 )
scala> val sims = Dglm.simulateRegular(mod, params, 1.0).
| steps.take(100).map(_._1).toVector
<console>:12: warning: Unused import
import breeze.linalg.{DenseMatrix, DenseVector}
^
<console>:12: warning: Unused import
import breeze.linalg.{DenseMatrix, DenseVector}
^
sims: Vector[com.github.jonnylaw.dlm.Data] = Vector(Data(1.0,DenseVector(Some(3.0))), Data(2.0,DenseVector(Some(1.0))), Data(3.0,DenseVector(Some(1.0))), Data(4.0,DenseVector(Some(3.0))), Data(5.0,DenseVector(Some(1.0))), Data(6.0,DenseVector(Some(1.0))), Data(7.0,DenseVector(Some(0.0))), Data(8.0,DenseVector(Some(3.0))), Data(9.0,DenseVector(Some(3.0))), Data(10.0,DenseVector(Some(1.0))), Data(11.0,DenseVector(Some(1.0))), Data(12.0,DenseVector(Some(5.0))), Data(13.0,DenseVector(Some(2.0))), Data(14.0,DenseVector(Some(1.0))), Data(15.0,DenseVector(Some(4.0))), Data(16.0,DenseVector(Some(3.0))), Data(17.0,DenseVector(Some(2.0))), Data(18.0,DenseVector(Some(1.0))), Data(19.0,DenseVector(Some(2.0))), Data(20.0,DenseVector(Some(0.0))), Data(21.0,DenseVector(Some(1.0))), Data(22.0,DenseVect...A simulation from the Poisson model is presented in the figure below:
Then the bootstrap particle filter can be user to determine the latent-state of the partially observed inhomogeneous Poisson process:
scala> val n = 500
<console>:10: warning: Unused import
import com.github.jonnylaw.dlm._
^
<console>:12: warning: Unused import
import breeze.linalg.{DenseMatrix, DenseVector}
^
<console>:12: warning: Unused import
import breeze.linalg.{DenseMatrix, DenseVector}
^
n: Int = 500
scala> val filtered = ParticleFilter(n, n, ParticleFilter.metropolisResampling(10)).filter(mod, sims, params)
<console>:12: warning: Unused import
import breeze.linalg.{DenseMatrix, DenseVector}
^
<console>:12: warning: Unused import
import breeze.linalg.{DenseMatrix, DenseVector}
^
filtered: Vector[com.github.jonnylaw.dlm.PfState] = Vector(PfState(1.0,Vector(DenseVector(-2.3625025642777553), DenseVector(-0.4688224197480089), DenseVector(-1.162930600402134), DenseVector(0.9559243041815069), DenseVector(-0.30087496731122293), DenseVector(-0.37320994915993966), DenseVector(-0.9885270620967307), DenseVector(-1.3840414978356388), DenseVector(0.42062251817069607), DenseVector(-0.08129736644848612), DenseVector(1.1339618573822159), DenseVector(1.5243039004645105), DenseVector(0.20724178059092366), DenseVector(0.42514477645067666), DenseVector(0.8783650838210421), DenseVector(0.6805113607226727), DenseVector(-0.5598110723253853), DenseVector(-0.27956417907092973), DenseVector(-0.09300491121358057), DenseVector(-1.9330734556477431), DenseVector(1.484717472568114), DenseVec...