|
| 1 | +############################################################################## |
| 2 | +# |
| 3 | +# Example of fully-adapted particle filtering |
| 4 | +# in a linear Gaussian state space model |
| 5 | +# |
| 6 | +# This program is free software; you can redistribute it and/or modify |
| 7 | +# it under the terms of the GNU General Public License as published by |
| 8 | +# the Free Software Foundation; either version 2 of the License, or |
| 9 | +# (at your option) any later version. |
| 10 | +# |
| 11 | +# This program is distributed in the hope that it will be useful, |
| 12 | +# but WITHOUT ANY WARRANTY; without even the implied warranty of |
| 13 | +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
| 14 | +# GNU General Public License for more details. |
| 15 | +# |
| 16 | +# You should have received a copy of the GNU General Public License along |
| 17 | +# with this program; if not, write to the Free Software Foundation, Inc., |
| 18 | +# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. |
| 19 | +# |
| 20 | +############################################################################## |
| 21 | + |
| 22 | +#' State estimation in a linear Gaussian state space model |
| 23 | +#' |
| 24 | +#' @description |
| 25 | +#' Minimal working example of state estimation in a linear Gaussian state |
| 26 | +#' space model using Kalman filtering and a fully-adapted particle filter. |
| 27 | +#' The code estimates the bias and mean squared error (compared with the |
| 28 | +#' Kalman estimate) while varying the number of particles in the particle |
| 29 | +#' filter. |
| 30 | +#' @details |
| 31 | +#' The Kalman filter is a standard implementation without an input. The |
| 32 | +#' particle filter is fully adapted (i.e. takes the current observation into |
| 33 | +#' account when proposing new particles and computing the weights). |
| 34 | +#' @return |
| 35 | +#' Returns a plot with the generated observations y and the difference in the |
| 36 | +#' state estimates obtained by the Kalman filter (the optimal solution) and |
| 37 | +#' the particle filter (with 20 particles). Furthermore, the function returns |
| 38 | +#' plots of the estimated bias and mean squared error of the state estimate |
| 39 | +#' obtained using the particle filter (while varying the number of particles) |
| 40 | +#' and the Kalman estimates. |
| 41 | +#' |
| 42 | +#' The function returns a list with the elements: |
| 43 | +#' \itemize{ |
| 44 | +#' \item{y: The observations generated from the model.} |
| 45 | +#' \item{x: The states generated from the model.} |
| 46 | +#' \item{xhatKF: The estimate of the state from the Kalman filter.} |
| 47 | +#' \item{xhatPF: The estimate of the state from the particle filter with |
| 48 | +#' 20 particles.} |
| 49 | +#' } |
| 50 | +#' @references |
| 51 | +#' Dahlin, J. & Schoen, T. B. "Getting started with particle |
| 52 | +#' Metropolis-Hastings for inference in nonlinear dynamical models." |
| 53 | +#' pre-print, arXiv:1511.01707, 2015. |
| 54 | +#' @author |
| 55 | +#' Johan Dahlin <johan.dahlin@liu.se> |
| 56 | +#' @note |
| 57 | +#' See Section 4.2 in the reference for more details. |
| 58 | +#' @keywords |
| 59 | +#' misc |
| 60 | +#' @export |
| 61 | +#' @importFrom grDevices col2rgb |
| 62 | +#' @importFrom grDevices rgb |
| 63 | +#' @importFrom graphics abline |
| 64 | +#' @importFrom graphics hist |
| 65 | +#' @importFrom graphics layout |
| 66 | +#' @importFrom graphics lines |
| 67 | +#' @importFrom graphics par |
| 68 | +#' @importFrom graphics plot |
| 69 | +#' @importFrom graphics points |
| 70 | +#' @importFrom stats acf |
| 71 | +#' @importFrom stats density |
| 72 | +#' @importFrom stats sd |
| 73 | +#' @importFrom stats var |
| 74 | + |
| 75 | +example1_lgss <- function() { |
| 76 | + |
| 77 | + # Set the random seed to replicate results in tutorial |
| 78 | + set.seed( 10 ) |
| 79 | + |
| 80 | + ############################################################################## |
| 81 | + # Define the model |
| 82 | + ############################################################################## |
| 83 | + |
| 84 | + # Here, we use the following model |
| 85 | + # |
| 86 | + # x[tt+1] = phi * x[tt] + sigmav * v[tt] |
| 87 | + # y[tt] = x[tt] + sigmae * e[tt] |
| 88 | + # |
| 89 | + # where v[tt] ~ N(0,1) and e[tt] ~ N(0,1) |
| 90 | + |
| 91 | + # Set the parameters of the model |
| 92 | + phi = 0.75 |
| 93 | + sigmav = 1.00 |
| 94 | + sigmae = 0.10 |
| 95 | + |
| 96 | + # Set the number of time steps to simulate |
| 97 | + T = 250 |
| 98 | + |
| 99 | + # Set the initial state |
| 100 | + x0 = 0 |
| 101 | + |
| 102 | + ############################################################################## |
| 103 | + # Generate data |
| 104 | + ############################################################################## |
| 105 | + |
| 106 | + data <- generateData( phi, sigmav, sigmae, T, x0 ) |
| 107 | + x <- data$x |
| 108 | + y <- data$y |
| 109 | + |
| 110 | + # Plot the latent state and observations |
| 111 | + layout( matrix( c(1,1,2,2,3,4), 3, 2, byrow = TRUE ) ) |
| 112 | + par ( mar = c(4,5,0,0) ) |
| 113 | + |
| 114 | + grid <- seq( 0, T ) |
| 115 | + |
| 116 | + plot( grid, y, col="#1B9E77", type="l", xlab="time", ylab="observation", ylim=c(-6,6), bty="n") |
| 117 | + |
| 118 | + ################################################################################### |
| 119 | + # State estimation using the particle filter and Kalman filter |
| 120 | + ################################################################################### |
| 121 | + |
| 122 | + # Using N = 20 particles and plot the estimate of the latent state |
| 123 | + N <- 20 |
| 124 | + outPF <- sm( y, phi, sigmav, sigmae, N, T, x0 ) |
| 125 | + outKF <- kf( y, phi, sigmav, sigmae, x0, 0.01 ) |
| 126 | + diff <- outPF$xh - outKF$xh[-(T+1)] |
| 127 | + |
| 128 | + grid <- seq( 0, T-1 ) |
| 129 | + plot( grid, diff, col="#7570B3", type="l", xlab="time", ylab="difference in state estimate", ylim=c(-0.1,0.1), bty="n") |
| 130 | + |
| 131 | + # Compute bias and MSE |
| 132 | + logBiasMSE = matrix( 0, nrow = 7, ncol = 2 ) |
| 133 | + gridN = c( 10, 20, 50, 100, 200, 500, 1000) |
| 134 | + |
| 135 | + for ( ii in 1:length(gridN) ) { |
| 136 | + smEstimate <- sm(y,phi,sigmav,sigmae,gridN[ii],T,x0)$xh |
| 137 | + kmEstimate <- outKF$xh[-(T+1)] |
| 138 | + |
| 139 | + logBiasMSE[ ii, 1 ] <- log( mean( abs( smEstimate - kmEstimate ) ) ) |
| 140 | + logBiasMSE[ ii, 2 ] <- log( mean( ( smEstimate - kmEstimate )^2 ) ) |
| 141 | + } |
| 142 | + |
| 143 | + # Plot the bias and MSE for comparison |
| 144 | + plot( gridN, logBiasMSE[ , 1 ], col="#E7298A", type="l", xlab="no. particles (N)", ylab="log-bias", ylim=c(-7,-3), bty="n") |
| 145 | + points( gridN, logBiasMSE[ , 1 ], col="#E7298A", pch=19 ) |
| 146 | + |
| 147 | + plot( gridN, logBiasMSE[ , 2 ], col="#66A61E", lwd=1.5, type="l", xlab="no. particles (N)", ylab="log-MSE", ylim=c(-12,-6), bty="n") |
| 148 | + points( gridN, logBiasMSE[ , 2 ], col="#66A61E", pch=19 ) |
| 149 | + |
| 150 | + list(y=y, x=x, xhatKF=outKF$xh[-(T+1)], xhatPF=outPF$xh) |
| 151 | +} |
| 152 | + |
| 153 | +################################################################################### |
| 154 | +# End of file |
| 155 | +################################################################################### |
0 commit comments