Online Parameter Inference
Given a DLM, with observations \(Y_{1:N}\) and latent state \(x_{0:N}\) with observation and state equations:
\[\begin{align*} Y_t &= F_t^Tx_t + v_t, & v_t &\sim \mathcal{N}(0, V_t), \\ X_t &= G_tx_{t-1} + w_t, & w_t &\sim \mathcal{N}(0, W_t), \\ x_0 &\sim \mathcal{N}(m_0, C_0) \end{align*}\]We can joint posterior distribution of the parameters and latent-state, \(p(\psi_t, x_t|y_{1:N})\) using on of the following online algorithms. These are useful for streaming data as they donโt require re-running of costly MCMC algorithms as each new data point is observed.
The online filters will be demonstrated for a first order polynomial DLM, with \(V = 2, W = 3\), plotted in the figure below.
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 = Dlm(
| f = (t: Double) => DenseMatrix((1.0)),
| g = (t: Double) => DenseMatrix((1.0))
| )
<console>:12: warning: Unused import
import breeze.linalg.{DenseMatrix, DenseVector}
^
mod: com.github.jonnylaw.dlm.Dlm = Dlm($$Lambda$7176/1387060605@3561f0e4,$$Lambda$7177/873481534@dfde1d2)
scala> val p = DlmParameters(v = DenseMatrix(2.0),
| w = DenseMatrix(3.0),
| DenseVector(0.0),
| DenseMatrix(1.0))
p: com.github.jonnylaw.dlm.DlmParameters = DlmParameters(2.0 ,3.0 ,DenseVector(0.0),1.0 )
scala> val data = Dlm.simulateRegular(mod, p, 1.0).
| steps.
| take(100).
| toVector.
| map(_._1)
<console>:12: warning: Unused import
import breeze.linalg.{DenseMatrix, DenseVector}
^
<console>:12: warning: Unused import
import breeze.linalg.{DenseMatrix, DenseVector}
^
data: scala.collection.immutable.Vector[com.github.jonnylaw.dlm.Data] = Vector(Data(1.0,DenseVector(Some(-2.5293023867229123))), Data(2.0,DenseVector(Some(-1.6779925288016404))), Data(3.0,DenseVector(Some(4.226481325906972))), Data(4.0,DenseVector(Some(4.6912843209107535))), Data(5.0,DenseVector(Some(4.783935988541616))), Data(6.0,DenseVector(Some(5.870163707104027))), Data(7.0,DenseVector(Some(5.86367115767524))), Data(8.0,DenseVector(Some(5.321197241005569))), Data(9.0,DenseVector(Some(5.775779586003766))), Data(10.0,DenseVector(Some(5.914323679662503))), Data(11.0,DenseVector(Some(6.769272425760471))), Data(12.0,DenseVector(Some(8.32019363232675))), Data(13.0,DenseVector(Some(6.430201264460451))), Data(14.0,DenseVector(Some(5.440527789674329))), Data(15.0,DenseVector(Some(8.339597663...
Conjugate Filter
The Conjugate filter is an exact filter, which augments the Kalman Filter to learn the observation variance using a conjugate update of the Gamma distribution.
scala> val prior = InverseGamma(3.0, 4.0)
<console>:12: warning: Unused import
import breeze.linalg.{DenseMatrix, DenseVector}
^
<console>:12: warning: Unused import
import breeze.linalg.{DenseMatrix, DenseVector}
^
prior: com.github.jonnylaw.dlm.InverseGamma = InverseGamma(3.0,4.0)
scala> val filtered = ConjugateFilter(prior, ConjugateFilter.advanceState(p, mod.g)).filter(mod, data, p)
<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.InverseGammaState] = Vector(InverseGammaState(KfState(0.0,DenseVector(0.0),1.0 ,DenseVector(0.0),1.0 ,None,None),Vector(InverseGamma(3.0,4.0))), InverseGammaState(KfState(1.0,DenseVector(-1.6862015911486081),1.3333333333333335 ,DenseVector(0.0),4.0 ,Some(DenseVector(0.0)),Some(6.0 )),Vector(InverseGamma(4.0,6.132456854494073))), InverseGammaState(KfState(2.0,DenseVector(-1.680623749718091),1.3889475828905695 ,DenseVector(-1.6862015911486081),4.333333333333334 ,Some(DenseVector(-1.6862015911486081)),Some(6.377485618164691 )),Vector(InverseGamma(5.0,6.1324784543506246))), InverseGammaState(KfState(3.0,DenseVector(2.697235143703967),1.1362217616070254 ,DenseVector(-1.680623749718091),4.3889475828905695 ,Some(DenseVector(-1.680623749718091)...
Liu and West Filter
The Liu and West filter can be used for general state space models:
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 n0 = 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}
^
n0: Int = 500
scala> // smoothing parameter for the mixture of gaussians, equal to (3 delta - 1) / 2 delta
| val a = (3 * 0.95 - 1) / 2 * 0.95
<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}
^
a: Double = 0.8787499999999998
scala> val prior = for {
| v <- InverseGamma(3.0, 4.0)
| w <- InverseGamma(3.0, 10.0)
| } yield DlmParameters(DenseMatrix(v), DenseMatrix(w), p.m0, p.c0)
<console>:12: warning: Unused import
import breeze.linalg.{DenseMatrix, DenseVector}
^
prior: breeze.stats.distributions.Rand[com.github.jonnylaw.dlm.DlmParameters] = FlatMappedRand(InverseGamma(3.0,4.0),$$Lambda$7192/1145516293@3abaa626)
scala> val filtered = LiuAndWestFilter(n, prior, a, n0).filter(mod, data, p)
<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.PfStateParams] = Vector(PfStateParams(0.0,Vector(DenseVector(1.0577222738067933), DenseVector(0.65907619291104), DenseVector(-0.17823751304632046), DenseVector(0.035659525912833376), DenseVector(1.2117753447354738), DenseVector(-1.5310063344627534), DenseVector(1.3436117931016343), DenseVector(-0.5682146672689058), DenseVector(-0.007146604723922701), DenseVector(1.695918803155607), DenseVector(-0.9668585825692261), DenseVector(0.2717232377169513), DenseVector(-0.21173366560881243), DenseVector(0.15021711747480504), DenseVector(0.4104937700280198), DenseVector(0.43038621577598307), DenseVector(-0.6693304793505566), DenseVector(-0.2848665305125076), DenseVector(-0.5314941308853947), DenseVector(-1.83509945985119), DenseVector(0.5844060246342688), D...
RB Filter
The Rao-Blackwell filter marginalises the Gaussian part of the DLM in order to perform filtering using the Kalman filter conditional on sampled values of the static parameters. Parameter distributions at each time point are approximated by a weighted particle cloud. A Kalman filter is performed for conditional on the value of each particle. Thereby determining the join posterior of the parameters and latent state.
scala> // smoothing parameter for the mixture of gaussians
| val delta = 0.99
<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}
^
delta: Double = 0.99
scala> val a = (3 * delta - 1) / 2 * delta
<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}
^
a: Double = 0.9751499999999999
scala> val prior = for {
| v <- InverseGamma(3.0, 4.0)
| w <- InverseGamma(3.0, 10.0)
| } yield DlmParameters(DenseMatrix(v), DenseMatrix(w), p.m0, p.c0)
<console>:12: warning: Unused import
import breeze.linalg.{DenseMatrix, DenseVector}
^
prior: breeze.stats.distributions.Rand[com.github.jonnylaw.dlm.DlmParameters] = FlatMappedRand(InverseGamma(3.0,4.0),$$Lambda$7245/573816379@30852497)
scala> val filtered = RaoBlackwellFilter(500, prior, a, 250).filter(mod, data, p)
<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.RbState] = Vector(RbState(0.0,Vector(DlmParameters(0.07099092072515194 ,1.7673292858471434 ,DenseVector(-Infinity),0.0 ), DlmParameters(0.45044177682675507 ,1.312731169788882 ,DenseVector(-Infinity),0.0 ), DlmParameters(-0.14838537656764245 ,2.728102242428066 ,DenseVector(-Infinity),0.0 ), DlmParameters(-0.21394137993869727 ,1.2170529583303318 ,DenseVector(-Infinity),0.0 ), DlmParameters(0.8110134577366094 ,1.2526742132506905 ,DenseVector(-Infinity),0.0 ), DlmParameters(-0.9807177792059666 ,1.3161701464877373 ,DenseVector(-Infinity),0.0 ), DlmParameters(0.2426307104889515 ,2.798466333356446 ,DenseVector(-Infinity),0.0 ), DlmParameters(0.19361330542672445 ,1.4046223486343579 ,DenseVector(-Infinity),0.0 ), DlmParameters(0.55...
Storvik Filter
Using conjugate updating for \(V\) and \(W\):
scala> val priorV = InverseGamma(3.0, 4.0)
<console>:12: warning: Unused import
import breeze.linalg.{DenseMatrix, DenseVector}
^
<console>:12: warning: Unused import
import breeze.linalg.{DenseMatrix, DenseVector}
^
priorV: com.github.jonnylaw.dlm.InverseGamma = InverseGamma(3.0,4.0)
scala> val priorW = InverseGamma(3.0, 10.0)
<console>:12: warning: Unused import
import breeze.linalg.{DenseMatrix, DenseVector}
^
<console>:12: warning: Unused import
import breeze.linalg.{DenseMatrix, DenseVector}
^
priorW: com.github.jonnylaw.dlm.InverseGamma = InverseGamma(3.0,10.0)
scala> val filtered = StorvikFilter.filterTs(mod, data, p, priorV, priorW, n, n)
<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.StorvikState] = Vector(StorvikState(0.0,Vector(DenseVector(1.1054674388251258), DenseVector(-1.1288317927649705), DenseVector(0.07231666137262859), DenseVector(-1.8847571229918934), DenseVector(0.42279777154662246), DenseVector(1.199060658542842), DenseVector(0.3376020298001149), DenseVector(0.3129548590826655), DenseVector(-0.14847668805838232), DenseVector(1.7992087122295282), DenseVector(1.206819106098928), DenseVector(0.4021683104747289), DenseVector(0.9713168560299944), DenseVector(0.690403791858766), DenseVector(1.4464506920705553), DenseVector(-0.5865709893639524), DenseVector(0.17738054742852583), DenseVector(1.488091877582659), DenseVector(0.8462074579529251), DenseVector(-0.14495894712479473), DenseVector(-0.051931569469057796), DenseVe...