<div>
<p>This example shows how to construct and conduct inference on a state space model using particle filtering algorithms. <code>nimble</code> currently has versions of the bootstrap filter, the auxiliary particle filter, the ensemble Kalman filter, and the Liu and West filter implemented. Additionally, particle MCMC samplers are available and can be specified for both univariate and multivariate parameters.</p>
<div>
<h1>Model Creation</h1>
<p>
Assume <img title="x_{i}" alt="x_{i}" class="latex" src="http://s0.wp.com/latex.php?latex=x_%7Bi%7D&bg=ffffff&%23038;fg=000000&%23038;s=0"> is the latent state and <img title="y_{i}" alt="y_{i}" class="latex" src="http://s0.wp.com/latex.php?latex=y_%7Bi%7D&bg=ffffff&%23038;fg=000000&%23038;s=0"> is the observation at time <img title="i" alt="i" class="latex" src="http://s0.wp.com/latex.php?latex=i&bg=ffffff&%23038;fg=000000&%23038;s=0"> for <img title="i=1,ldots,t" alt="i=1,ldots,t" class="latex" src="http://s0.wp.com/latex.php?latex=i%3D1%2C%5Cldots%2Ct&bg=ffffff&%23038;fg=000000&%23038;s=0">. We define our state space model as </p>
<p><img title="x_{i} sim N(a cdot x_{i-1} + b, sigma_{PN})" alt="x_{i} sim N(a cdot x_{i-1} + b, sigma_{PN})" class="latex" src="http://s0.wp.com/latex.php?latex=x_%7Bi%7D+%5Csim+N%28a+%5Ccdot+x_%7Bi-1%7D+%2B+b%2C+%5Csigma_%7BPN%7D%29&bg=ffffff&%23038;fg=000000&%23038;s=0"><br>
<img title=" y_{i} sim t_{5}(x_{i}, sigma_{OE})" alt=" y_{i} sim t_{5}(x_{i}, sigma_{OE})" class="latex" src="http://s0.wp.com/latex.php?latex=+y_%7Bi%7D+%5Csim+t_%7B5%7D%28x_%7Bi%7D%2C+%5Csigma_%7BOE%7D%29&bg=ffffff&%23038;fg=000000&%23038;s=0"></p>
<p> with initial states </p>
<p><img title="  x_{1} sim N(frac{b}{1-a}, frac{sigma_{PN}}{sqrt{1-a^2}})" alt="  x_{1} sim N(frac{b}{1-a}, frac{sigma_{PN}}{sqrt{1-a^2}})" class="latex" src="http://s0.wp.com/latex.php?latex=++x_%7B1%7D+%5Csim+N%28%5Cfrac%7Bb%7D%7B1-a%7D%2C+%5Cfrac%7B%5Csigma_%7BPN%7D%7D%7B%5Csqrt%7B1-a%5E2%7D%7D%29&bg=ffffff&%23038;fg=000000&%23038;s=0"><br>
<img title=" y_{1} sim t_{5}(x_{1}, sigma_{OE})" alt=" y_{1} sim t_{5}(x_{1}, sigma_{OE})" class="latex" src="http://s0.wp.com/latex.php?latex=+y_%7B1%7D+%5Csim+t_%7B5%7D%28x_%7B1%7D%2C+%5Csigma_%7BOE%7D%29&bg=ffffff&%23038;fg=000000&%23038;s=0"></p>
<p> and prior distributions </p>
<p> <img title=" a sim Unif(-0.999, 0.999) " alt=" a sim Unif(-0.999, 0.999) " class="latex" src="http://s0.wp.com/latex.php?latex=+a+%5Csim+Unif%28-0.999%2C+0.999%29+&bg=ffffff&%23038;fg=000000&%23038;s=0"><br>
<img title=" b sim N(0, 1000) " alt=" b sim N(0, 1000) " class="latex" src="http://s0.wp.com/latex.php?latex=+b+%5Csim+N%280%2C+1000%29+&bg=ffffff&%23038;fg=000000&%23038;s=0"></p>
<p> where <img title="N(mu, sigma)" alt="N(mu, sigma)" class="latex" src="http://s0.wp.com/latex.php?latex=N%28%5Cmu%2C+%5Csigma%29&bg=ffffff&%23038;fg=000000&%23038;s=0"> denotes a normal distribution with mean <img title="mu" alt="mu" class="latex" src="http://s0.wp.com/latex.php?latex=%5Cmu&bg=ffffff&%23038;fg=000000&%23038;s=0"> and standard deviation <img title="sigma" alt="sigma" class="latex" src="http://s0.wp.com/latex.php?latex=%5Csigma&bg=ffffff&%23038;fg=000000&%23038;s=0">, and <img title="t_{nu}(mu, sigma)" alt="t_{nu}(mu, sigma)" class="latex" src="http://s0.wp.com/latex.php?latex=t_%7B%5Cnu%7D%28%5Cmu%2C+%5Csigma%29&bg=ffffff&%23038;fg=000000&%23038;s=0"> is a shifted, scaled <img title="t" alt="t" class="latex" src="http://s0.wp.com/latex.php?latex=t&bg=ffffff&%23038;fg=000000&%23038;s=0">-distribution with center parameter <img title="mu" alt="mu" class="latex" src="http://s0.wp.com/latex.php?latex=%5Cmu&bg=ffffff&%23038;fg=000000&%23038;s=0">, scale parameter <img title="sigma" alt="sigma" class="latex" src="http://s0.wp.com/latex.php?latex=%5Csigma&bg=ffffff&%23038;fg=000000&%23038;s=0">, and <img title="nu" alt="nu" class="latex" src="http://s0.wp.com/latex.php?latex=%5Cnu&bg=ffffff&%23038;fg=000000&%23038;s=0"> degrees of freedom. </p>
<p> We specify and build our state space model below, using <img title="t=10" alt="t=10" class="latex" src="http://s0.wp.com/latex.php?latex=t%3D10&bg=ffffff&%23038;fg=000000&%23038;s=0"> time points:</p>
<div>
<div>
<div>
<pre><span>## load the nimble library and set seed</span>
<span>library</span><span>(</span><span>'nimble'</span><span>)</span>
<span>set.seed</span><span>(</span><span>1</span><span>)</span>

<span>## define the model</span>
<span>stateSpaceCode</span> <span>&lt;-</span> <span>nimbleCode</span><span>({</span>
    <span>a</span> <span>~</span> <span>dunif</span><span>(</span><span>-</span><span>0.9999</span><span>,</span> <span>0.9999</span><span>)</span>
    <span>b</span> <span>~</span> <span>dnorm</span><span>(</span><span>0</span><span>,</span> <span>sd</span> <span>=</span> <span>1000</span><span>)</span>
    <span>sigPN</span> <span>~</span> <span>dunif</span><span>(</span><span>1e-04</span><span>,</span> <span>1</span><span>)</span>
    <span>sigOE</span> <span>~</span> <span>dunif</span><span>(</span><span>1e-04</span><span>,</span> <span>1</span><span>)</span>
    <span>x[</span><span>1</span><span>]</span> <span>~</span> <span>dnorm</span><span>(b</span><span>/</span><span>(</span><span>1</span> <span>-</span> <span>a),</span> <span>sd</span> <span>= sigPN</span><span>/</span><span>sqrt</span><span>((</span><span>1</span><span>-</span><span>a</span><span>*</span><span>a)))</span>
    <span>y[</span><span>1</span><span>]</span> <span>~</span> <span>dt</span><span>(</span><span>mu</span> <span>= x[</span><span>1</span><span>],</span> <span>sigma</span> <span>= sigOE,</span> <span>df</span> <span>=</span> <span>5</span><span>)</span>
    <span>for</span> <span>(i</span> <span>in</span> <span>2</span><span>:</span><span>t) {</span>
        <span>x[i]</span> <span>~</span> <span>dnorm</span><span>(a</span> <span>*</span> <span>x[i</span> <span>-</span> <span>1</span><span>]</span> <span>+</span> <span>b,</span> <span>sd</span> <span>= sigPN)</span>
        <span>y[i]</span> <span>~</span> <span>dt</span><span>(</span><span>mu</span> <span>= x[i],</span> <span>sigma</span> <span>= sigOE,</span> <span>df</span> <span>=</span> <span>5</span><span>)</span>
    <span>}</span>
<span>})</span>

<span>## define data, constants, and initial values </span>
<span>data</span> <span>&lt;-</span> <span>list</span><span>(</span>
    <span>y</span> <span>=</span> <span>c</span><span>(</span><span>0.213</span><span>,</span> <span>1.025</span><span>,</span> <span>0.314</span><span>,</span> <span>0.521</span><span>,</span> <span>0.895</span><span>,</span> <span>1.74</span><span>,</span> <span>0.078</span><span>,</span> <span>0.474</span><span>,</span> <span>0.656</span><span>,</span> <span>0.802</span><span>)</span>
<span>)</span>
<span>constants</span> <span>&lt;-</span> <span>list</span><span>(</span>
    <span>t</span> <span>=</span> <span>10</span>
<span>)</span>
<span>inits</span> <span>&lt;-</span> <span>list</span><span>(</span>
    <span>a</span> <span>=</span> <span>0</span><span>,</span>
    <span>b</span> <span>=</span> <span>.5</span><span>,</span>
    <span>sigPN</span> <span>=</span> <span>.1</span><span>,</span>
    <span>sigOE</span> <span>=</span> <span>.05</span>
<span>)</span>

<span>## build the model</span>
<span>stateSpaceModel</span> <span>&lt;-</span> <span>nimbleModel</span><span>(stateSpaceCode,</span>
                              <span>data</span> <span>= data,</span>
                              <span>constants</span> <span>= constants,</span>
                              <span>inits</span> <span>= inits,</span>
                              <span>check</span> <span>=</span> <span>FALSE</span><span>)</span>
</pre>
</div>
<div>
<pre>## defining model...
</pre>
</div>
<div>
<pre>## building model...
</pre>
</div>
<div>
<pre>## setting data and initial values...
</pre>
</div>
<div>
<pre>## running calculate on model (any error reports that follow may simply
## reflect missing values in model variables) ...
</pre>
</div>
<div>
<pre>## 
</pre>
</div>
<div>
<pre>## checking model sizes and dimensions...
</pre>
</div>
<div>
<pre>## note that missing values (NAs) or non-finite values were found in model
## variables: x, lifted_a_times_x_oBi_minus_1_cB_plus_b. This is not an error,
## but some or all variables may need to be initialized for certain algorithms
## to operate properly.
</pre>
</div>
<div>
<pre>## 
</pre>
</div>
<div>
<pre>## model building finished.
</pre>
</div>
</div>
</div>
</div>
<div>
<h1>Construct and run a bootstrap filter</h1>
<p>We next construct a bootstrap filter to conduct inference on the latent states of our state space model. Note that the bootstrap filter, along with the auxiliary particle filter and the ensemble Kalman filter, treat the top-level parameters <code>a, b, sigPN</code>, and <code>sigOE</code> as fixed. Therefore, the bootstrap filter below will proceed as though <code>a = 0, b = .5, sigPN = .1</code>, and <code>sigOE = .05</code>, which are the initial values that were assigned to the top-level parameters.</p>
<p>The bootstrap filter takes as arguments the name of the model and the name of the latent state variable within the model. The filter can also take a control list that can be used to fine-tune the algorithm’s configuration.</p>
<div>
<div>
<div>
<pre><span>## build bootstrap filter and compile model and filter</span>
<span>bootstrapFilter</span> <span>&lt;-</span> <span>buildBootstrapFilter</span><span>(stateSpaceModel,</span> <span>nodes</span> <span>=</span> <span>'x'</span><span>)</span>
<span>compiledList</span> <span>&lt;-</span> <span>compileNimble</span><span>(stateSpaceModel, bootstrapFilter)</span>
</pre>
</div>
<div>
<pre>## compiling... this may take a minute. Use 'showCompilerOutput = TRUE' to see C++ compiler details.
</pre>
</div>
<div>
<pre>## compilation finished.
</pre>
</div>
<div>
<pre><span>## run compiled filter with 10,000 particles. </span>
<span>## note that the bootstrap filter returns an estimate of the log-likelihood of the model.</span>
<span>compiledList</span><span>$</span><span>bootstrapFilter</span><span>$</span><span>run</span><span>(</span><span>10000</span><span>)</span>
</pre>
</div>
<div>
<pre>## [1] -28.13009
</pre>
</div>
</div>
</div>
<p>Particle filtering algorithms in <code>nimble</code> store weighted samples of the filtering distribution of the latent states in the <code>mvSamples</code> modelValues object. Equally weighted samples are stored in the <code>mvEWSamples</code> object. By default, <code>nimble</code> only stores samples from the final time point.</p>
<div>
<div>
<div>
<pre><span>## extract equally weighted posterior samples of x[10] and create a histogram</span>
<span>posteriorSamples</span> <span>&lt;-</span> <span>as.matrix</span><span>(compiledList</span><span>$</span><span>bootstrapFilter</span><span>$</span><span>mvEWSamples)</span>
<span>hist</span><span>(posteriorSamples)</span>
</pre>
</div>
</div>
<div><img title="plot of chunk plot-bootstrap" alt="plot of chunk plot-bootstrap" src="https://i1.wp.com/r-nimble.org/wp-content/uploads/2017/01/plot-bootstrap-1.png?w=456&ssl=1" data-recalc-dims="1"></div>
</div>
<p>The auxiliary particle filter and ensemble Kalman filter can be constructed and run in the same manner as the bootstrap filter.</p>
</div>
<div>
<h1>Conduct inference on top-level parameters using particle MCMC</h1>
<p>Particle MCMC can be used to conduct inference on the posterior distribution of both the latent states and any top-level parameters of interest in a state space model. The particle marginal Metropolis-Hastings sampler can be specified to jointly sample the <code>a, b, sigPN</code>, and <code>sigOE</code> top level parameters within <code>nimble</code>‘s MCMC framework as follows:</p>
<div>
<div>
<div>
<pre><span>## create MCMC specification for the state space model</span>
<span>stateSpaceMCMCconf</span> <span>&lt;-</span> <span>configureMCMC</span><span>(stateSpaceModel,</span> <span>nodes</span> <span>=</span> <span>NULL</span><span>)</span>

<span>## add a block pMCMC sampler for a, b, sigPN, and sigOE </span>
<span>stateSpaceMCMCconf</span><span>$</span><span>addSampler</span><span>(</span><span>target</span> <span>=</span> <span>c</span><span>(</span><span>'a'</span><span>,</span> <span>'b'</span><span>,</span> <span>'sigPN'</span><span>,</span> <span>'sigOE'</span><span>),</span>
                              <span>type</span> <span>=</span> <span>'RW_PF_block'</span><span>,</span> <span>control</span> <span>=</span> <span>list</span><span>(</span><span>latents</span> <span>=</span> <span>'x'</span><span>))</span>

<span>## build and compile pMCMC sampler</span>
<span>stateSpaceMCMC</span> <span>&lt;-</span> <span>buildMCMC</span><span>(stateSpaceMCMCconf)</span>
<span>compiledList</span> <span>&lt;-</span> <span>compileNimble</span><span>(stateSpaceModel, stateSpaceMCMC,</span> <span>resetFunctions</span> <span>=</span> <span>TRUE</span><span>)</span>
</pre>
</div>
<div>
<pre>## compiling... this may take a minute. Use 'showCompilerOutput = TRUE' to see C++ compiler details.
</pre>
</div>
<div>
<pre>## compilation finished.
</pre>
</div>
<div>
<pre><span>## run compiled sampler for 5000 iterations</span>
<span>compiledList</span><span>$</span><span>stateSpaceMCMC</span><span>$</span><span>run</span><span>(</span><span>5000</span><span>)</span>
</pre>
</div>
<div>
<pre>## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
</pre>
</div>
<div>
<pre>## NULL
</pre>
</div>
<div>
<pre><span>## create trace plots for each parameter</span>
<span>library</span><span>(</span><span>'coda'</span><span>)</span>
</pre>
</div>
<div>
<pre><span>par</span><span>(</span><span>mfrow</span> <span>=</span> <span>c</span><span>(</span><span>2</span><span>,</span><span>2</span><span>))</span>
<span>posteriorSamps</span> <span>&lt;-</span> <span>as.mcmc</span><span>(</span><span>as.matrix</span><span>(compiledList</span><span>$</span><span>stateSpaceMCMC</span><span>$</span><span>mvSamples))</span>
<span>traceplot</span><span>(posteriorSamps[,</span><span>'a'</span><span>],</span> <span>ylab</span> <span>=</span> <span>'a'</span><span>)</span>
<span>traceplot</span><span>(posteriorSamps[,</span><span>'b'</span><span>],</span> <span>ylab</span> <span>=</span> <span>'b'</span><span>)</span>
<span>traceplot</span><span>(posteriorSamps[,</span><span>'sigPN'</span><span>],</span> <span>ylab</span> <span>=</span> <span>'sigPN'</span><span>)</span>
<span>traceplot</span><span>(posteriorSamps[,</span><span>'sigOE'</span><span>],</span> <span>ylab</span> <span>=</span> <span>'sigOE'</span><span>)</span>
</pre>
</div>
</div>
<div><img title="plot of chunk run-PMCMC" alt="plot of chunk run-PMCMC" src="https://i0.wp.com/r-nimble.org/wp-content/uploads/2017/01/run-PMCMC-1.png?w=456&ssl=1" data-recalc-dims="1"></div>
</div>
<p>The above <code>RW_PF_block</code> sampler uses a multivariate normal proposal distribution to sample vectors of top-level parameters. To sample a scalar top-level parameter, use the <code>RW_PF</code> sampler instead. </p>
</div>
<p></p>
<p> <img height="1" width="1" src="https://r-nimble.org/?feed-stats-post-id=503"></p>
<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="https://r-nimble.org/building-particle-filters-and-particle-mcmc-in-nimble-2"> R – NIMBLE</a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a title="The R Project for Statistical Computing" rel="nofollow" href="https://www.r-project.org/">R</a> news and <a title="R tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/tutorial">tutorials</a> on topics such as: <a title="Data science" rel="nofollow" href="https://www.r-bloggers.com/search/data%20science">Data science</a>, <a title="Big Data" rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data">Big Data, <a title="R jobs" rel="nofollow" href="https://www.r-users.com/">R jobs</a>, visualization (<a title="ggplot and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2">ggplot2</a>, <a title="Boxplots using lattice and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/boxplot">Boxplots</a>, <a title="Maps and gis" rel="nofollow" href="https://www.r-bloggers.com/search/map">maps</a>, <a title="Animation in R" rel="nofollow" href="https://www.r-bloggers.com/search/animation">animation</a>), programming (<a title="RStudio IDE for R" rel="nofollow" href="https://www.r-bloggers.com/search/RStudio">RStudio</a>, <a title="Sweave and literate programming" rel="nofollow" href="https://www.r-bloggers.com/search/sweave">Sweave</a>, <a title="LaTeX in R" rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX">LaTeX</a>, <a title="SQL and databases" rel="nofollow" href="https://www.r-bloggers.com/search/SQL">SQL</a>, <a title="Eclipse IDE for R" rel="nofollow" href="https://www.r-bloggers.com/search/eclipse">Eclipse</a>, <a title="git and github, Version Control System" rel="nofollow" href="https://www.r-bloggers.com/search/git">git</a>, <a title="Large data in R using Hadoop" rel="nofollow" href="https://www.r-bloggers.com/search/hadoop">hadoop</a>, <a title="Web Scraping of google, facebook, yahoo, twitter and more using R" rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping">Web Scraping</a>) statistics (<a title="Regressions and ANOVA analysis tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/regression">regression</a>, <a title="principal component analysis tutorial" rel="nofollow" href="https://www.r-bloggers.com/search/PCA">PCA</a>, <a title="Time series" rel="nofollow" href="https://www.r-bloggers.com/search/time+series">time series</a>, <a title="finance trading" rel="nofollow" href="https://www.r-bloggers.com/search/trading">trading</a>) and more...
<img height="1" alt="" width="1" src="http://feeds.feedburner.com/~r/RBloggers/~4/yfGNWqsH9HQ"></a></div></div>