particleFilter: The particle filtering algorithm

In SMC: Sequential Monte Carlo (SMC) Algorithm

Description

Function for doing particle filtering given the state equation (via generateNextStreamFunc), and the observation equation density (via logObsDensFunc).

See the sections Details, Required Functions and Optional Functions for explanation on the arguments and the return values of the arguments that are themselves functions.

Usage

particleFilter(nStreams,                           
               nPeriods,                           
               dimPerPeriod,                       
               generateNextStreamsFunc,             
               logObsDensFunc,                     
               resampCriterionFunc = NULL,   
               resampFunc          = NULL,   
               summaryFunc         = NULL,   
               nMHSteps            = 0,      
               MHUpdateFunc        = NULL,   
               nStreamsPreResamp   = NULL,   
               returnStreams       = FALSE,  
               returnLogWeights    = FALSE,  
               verboseLevel        = 0,      
               ...)                                

Arguments

nStreams	integer > 0.
nPeriods	integer > 0.
dimPerPeriod	integer > 0.
generateNextStreamsFunc	function of six arguments (currentPeriod, lag1Streams, lag1LogWeights, streamIndices, startingStreams, ...).
logObsDensFunc	function of three arguments (currentPeriod, currentStreams, ...).
resampCriterionFunc	function of four arguments (currentPeriod, currentStreams, currentLogWeights, ...).
resampFunc	function of four arguments (currentPeriod, currentStreams, currentLogWeights, ...).
summaryFunc	function of four arguments (currentPeriod, currentStreams, currentLogWeights, ...).
nMHSteps	integer >= 0.
MHUpdateFunc	function of six arguments (currentPeriod, nMHSteps, currentStreams, lag1Streams, lag1LogWeights, ...).
nStreamsPreResamp	integer > 0.
returnStreams	logical.
returnLogWeights	logical.
verboseLevel	integer, a value >= 2 produces a lot of output.
...	optional arguments to be passed to generateNextStreamsFunc, logObsDensFunc, resampCriterionFunc, resampFunc, summaryFunc and MHUpdateFunc.

Details

We introduce the following terms, which will be used in the sections Required Function and Optional Function below:

stream

the state vector also called the particle, the hidden state or the latent variable. Below we will use the terms stream and state vector interchangeably.

dimPerPeriod

the dimension of the space, the state vectors live in.

Value

This function returns a list with the following components:

draws	a list with the following components: summary, propUniqueStreamIds, streams, logWeights, acceptanceRates. See the section Note for more details.
nStreams	the nStreams argument.
nPeriods	the nPeriods argument.
dimPerPeriod	the dimPerPeriod argument.
nStreamsPreResamp	the nStreamsPreResamp argument.
nMHSteps	the nMHSteps argument.
filterType	type of the filter: “particleFilter”.
time	the time taken by the run.

Optional function: generateNextStreamsFunc

Arguments:

The following argument(s) require some explanation:

lag1Streams

a matrix of dimension nStreams x dimPerPeriod of streams for currentPeriod - 1.

lag1LogWeights

a vector of length nStreams of log weights corresponding to the streams in the argument matrix lag1Streams.

streamIndices

a vector of length >= nStreams which are to be updated from currentPeriod - 1 to currentPeriod.

startingStreams

a matrix of dimension nStreams x dimPerPeriod to be used for currentPeriod = 1. If this is NULL, then the function should provide a way to generate streams for currentPeriod = 1.

Return value:

a matrix of dimension nStreamIndices x dimPerPeriod. The rows of this matrix contain the state vectors for period currentPeriod given the state vectors to be found in the streamIndices rows of the argument lag1Streams matrix. Here nStreamIndices is the length of the argument streamIndices.

Note:

this function should distinguish the cases currentPeriod == 1 and currentPeriod > 1 inside of it.

Optional function: logObsDensFunc

Arguments:

The following argument(s) require some explanation:

currentStreams

a matrix with dimPerPeriod columns, the rows containing the streams for currentPeriod.

Return value:

a vector of length nCurrentStreams, where nCurrentStreams refers to the number of rows of the currentStreams matrix argument. This vector contains the observation equation density values for currentPeriod in the log scale, evaluated at the rows of currentStreams.

Note:

nCurrentStreams might be >= nStreams.

Optional function: resampCriterionFunc

Arguments:

The following argument(s) require some explanation:

currentStreams

a matrix with dimPerPeriod columns, the rows containing the updated streams for currentPeriod.

currentLogWeights

a vector of log weights corresponding to the streams in the argument matrix currentStreams.

Return value:

TRUE or FALSE reflecting the decision of the resampling scheme implemented by this function.

Note:

The following points are in order:

–

resampling schemes manily depend on currentLogWeights, the other two arguments might come in handy for implementing period or stream specific resampling schemes.

–

if nStreamsPreResamp > nStreams, then this function should always return TRUE.

Optional function: resampFunc

Arguments:

see the sub-section Arguments: for section Optional function: resampCriterionFunc.

Return value:

a named list with the following components:

currentStreams

a matrix of dimension nStreams x dimPerPeriod. The rows of this matrix contain the streams for period currentPeriod + 1 that were resampled from those of the argument currentStreams matrix, which may contain >= nStreams rows.

currentLogWeights

The log weights vector of length nStreams, associated with the streams that were resampled in the returned currentStreams matrix. Note, after the resampling step, usually all the log weights are set to 0.

Note:

the components of the list returned by this function and the arguments to this function have two common names, namely, currentStreams and currentLogWeights. These entities have different meanings, as explained above. For example, the argument matrix currentStreams could possibly have >= nStreams rows, whereas the returned currentStreams has exactly nStreams number of (resampled) streams in its rows.

Optional function: summaryFunc

Arguments:

The following argument(s) require some explanation:

currentStreams

a matrix of dimension nStreams x dimPerPeriod of streams for currentPeriod.

currentLogWeights

a vector of log weights corresponding to the streams in the argument matrix currentStreams.

Return value:

a vector of length of dimSummPerPeriod of summaries for currentPeriod given the currentStreams and the currentLogWeights.

Optional function: MHUpdateFunc

Arguments:

The following argument(s) require some explanation:

nMHSteps

the number of Metropolis Hastings (MH) steps (iterations) to be performed.

currentStreams

a matrix of dimension nStreams x dimPerPeriod of streams for currentPeriod.

lag1Streams

a matrix of dimension nStreams x dimPerPeriod of streams for currentPeriod - 1.

lag1LogWeights

a vector of length nStreams of log weights corresponding to the streams in the argument matrix lag1Streams.

Return value:

a named list with the following components:

currentStreams

a matrix of dimension nStreams x dimPerPeriod. The rows of this matrix contain the streams for period currentPeriod that are (possibly) MH-updated versions of the rows of the argument currentStreams matrix.

acceptanceRates

a vector of length nStreams, representing the acceptance rates of the nMHSteps-many MH steps for each of the streams in the rows of the argument currentStreams matrix.

Note:

a positive value of nMHSteps performs as many MH steps on the rows of the argument currentStreams matrix. This is done to reduce the possible degeneracy after the resampling.

Warning

Using very small values (<= 1e3) for nStreams might not give reliable results.

Note

The effect of leaving the default value NULL for some of the arguments above are as follows:

resampCriterionFunc

the builtin resampling criterion, namely, resample when square of the coefficient of variation of the weights >= 1, is used.

resampFunc

the builtin resampling function, which resamples streams with probability proportional to their weights, is used.

summaryFunc

the builtin summary function, which returns the weighted average of each of the dimPerPeriod dimensions, is used.

MHUpdateFunc

the builtin Metropolis Hastings updating function, which generates proposals for currentPeriod streams using those of currentPeriod - 1, is used.

nStreamsPreResamp

it is set to nStreams.

Also, the following point is worth noting:

resampCriterionFunc, resampFunc, summaryFunc and MHUpdateFunc

are only necessary when user wants to try out new resampling schemes, enhanced summary generation procedures or more efficient MH updating rules, as part of their research. The default builtins take care of the typical problems.

This function returns a list with component called draw. The detailed description of this component, as promised in section Value, is as follows. It is a list itself with the following components:

summary

a matrix of dimension nPeriods x dimSummPerPeriod.

propUniqueStreamIds

a vector of length nPeriods. The values are either proportions of unique stream ids accpeted (at each period) if resampling was done or NA.

streams

an array of dimension nStreams x dimPerPeriod x nPeriods. This is returned if returnStreams = TRUE.

logWeights

a matrix of dimension nStreams x nPeriods. This is returned if returnLogWeights = TRUE.

acceptanceRates

a matrix of dimension nStreams x nPeriods. This is returned if nMHSteps > 0.

Author(s)

Gopi Goswami [e-mail]

References

Jun S. Liu (2001). Monte Carlo strategies for scientific computing. Springer. Page 66.

See Also

auxiliaryParticleFilter

Examples

MSObj  <- MarkovSwitchingFuncGenerator(-13579)
smcObj <-
    with(MSObj,
     {
         particleFilter(nStreams                = 5000,
                        nPeriods                = nrow(yy),
                        dimPerPeriod            = ncol(yy),
                        generateNextStreamsFunc = generateNextStreamsFunc,
                        logObsDensFunc          = logObsDensFunc,
                        returnStreams           = TRUE,
                        returnLogWeights        = TRUE,
                        verboseLevel            = 1)
     })
print(smcObj)
print(names(smcObj))
with(c(smcObj, MSObj),
 {
     par(mfcol = c(2, 1))
     plot(as.ts(yy),
          main     = expression('The data and the underlying regimes'),
          cex.main = 0.8,
          xlab     = 'period',
          ylab     = 'data and the regime means',
          cex.lab  = 0.8)
     lines(as.ts(mu), col = 2, lty = 2)
     plot(as.ts(draws$summary[1, ]),
          main     = expression('The underlying regimes and their estimates'),
          cex.main = 0.8,
          xlab     = 'period',
          ylab     = 'regime means',
          cex.lab  = 0.8)
     lines(as.ts(mu), col = 2, lty = 2)        
 })

MSObj  <- MarkovSwitchingFuncGenerator(-97531)
smcObj <-
    with(MSObj,
     {
         particleFilter(nStreams                = 5000,
                        nPeriods                = nrow(yy),
                        dimPerPeriod            = ncol(yy),
                        generateNextStreamsFunc = generateNextStreamsFunc,
                        logObsDensFunc          = logObsDensFunc,
                        nMHSteps                = 10,
                        returnStreams           = TRUE,
                        returnLogWeights        = TRUE,
                        verboseLevel            = 1)
     })
print(smcObj)
print(names(smcObj))
with(c(smcObj, MSObj),
 {
     par(mfcol = c(2, 1))
     plot(as.ts(yy),
          main     = expression('The data and the underlying regimes'),
          cex.main = 0.8,
          xlab     = 'period',
          ylab     = 'data and the regime means',
          cex.lab  = 0.8)
     lines(as.ts(mu), col = 2, lty = 2)
     plot(as.ts(draws$summary[1, ]),
          main     = expression('The underlying regimes and their estimates'),
          cex.main = 0.8,
          xlab     = 'period',
          ylab     = 'regime means',
          cex.lab  = 0.8)
     lines(as.ts(mu), col = 2, lty = 2)        
 })

Example output

##
## Sequential Monte Carlo Package (SMC)
##
## Functionality: sequential Monte Carlo (SMC) or sequential importance
## sampling (SIS) or hidden Markov models (HMM), particle filter (PF)
## and auxiliary particle filter (APF)
##
## Use: "help(package = SMC)" at the R prompt for more info
##
## Copyright (C) 2006-2020 Gopi Goswami
##
##    Created by: Gopi Goswami <[e-mail]>
## Maintained by: Gopi Goswami <[e-mail]>
##


BEGIN: SMC
..........
[Time to finish (est):        1 secs, this period:       10]....................
[Time to finish (est):        1 secs, this period:       30]....................
[Time to finish (est):        1 secs, this period:       50]....................
[Time to finish (est):        1 secs, this period:       70]....................
[Time to finish (est):        1 secs, this period:       90]..........
[Total time:          0 secs (usr),     0 secs (sys), this period:        100]
E N D: SMC

The settings for this particleFilter run:
                       
nStreams:          5000
nStreamsPreResamp: 5000
nPeriods:           100
dimPerPeriod:         1

The summary of the resampling proportions:
[Note: resampling was done for 12 periods]

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
 0.0450  0.1195  0.4140  0.3033  0.4244  0.4268      88 

 [1] "draws"                   "nStreams"               
 [3] "nPeriods"                "dimPerPeriod"           
 [5] "nStreamsPreResamp"       "nMHSteps"               
 [7] "summaryFunc"             "time"                   
 [9] "generateNextStreamsFunc" "logObsDensFunc"         
[11] "filterType"             

BEGIN: SMC
..........
[Time to finish (est):        1 secs, this period:       10]....................
[Time to finish (est):        1 secs, this period:       30]....................
[Time to finish (est):        1 secs, this period:       50]....................
[Time to finish (est):        1 secs, this period:       70]....................
[Time to finish (est):        1 secs, this period:       90]..........
[Total time:          1 secs (usr),     0 secs (sys), this period:        100]
E N D: SMC

The settings for this particleFilter run:
                       
nStreams:          5000
nStreamsPreResamp: 5000
nPeriods:           100
dimPerPeriod:         1
nMHSteps:            10

The summary of the resampling proportions:
[Note: resampling was done for 13 periods]

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
 0.0392  0.3220  0.4158  0.3419  0.4226  0.5042      87 

 [1] "draws"                   "nStreams"               
 [3] "nPeriods"                "dimPerPeriod"           
 [5] "nStreamsPreResamp"       "nMHSteps"               
 [7] "summaryFunc"             "time"                   
 [9] "generateNextStreamsFunc" "logObsDensFunc"         
[11] "filterType"             

SMC documentation built on May 1, 2019, 9:12 p.m.