Description Usage Arguments Details Value Required function: propagateFunc Optional function: resampCriterionFunc Optional function: resampFunc Optional function: summaryFunc Optional function: MHUpdateFunc Warning Note Author(s) References See Also Examples
Function for the doing sequential Monte Carlo algorithm given the
propagation rule over time (via propagateFunc
). This is the
most general interface for implementing a new SMC strategy, by
providing a new propagation rule.
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.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 
nStreams 

nPeriods 

dimPerPeriod 

propagateFunc 

resampCriterionFunc 

resampFunc 

summaryFunc 

nMHSteps 

MHUpdateFunc 

nStreamsPreResamp 

returnStreams 

returnLogWeights 

verboseLevel 

... 
optional arguments to be passed to 
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.
This function returns a list with the following components:
draws 
a list with the following components: 
nStreams 
the 
nPeriods 
the 
dimPerPeriod 
the 
nStreamsPreResamp 
the 
nMHSteps 
the 
filterType 
type of the filter: “sequentialMonteCarlo”. 
time 
the time taken by the run. 
The following argument(s) require some explanation:
nStreamsToGenerate
the number of streams to generate
for propagating from currentPeriod  1
to
currentPeriod
. This function is usally called by setting
nStreamsToGenerate
to nStreamsPreResamp
.
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
.
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
.
a named list with the following components:
currentStreams
a matrix of dimension
nStreamsToGenerate
x dimPerPeriod
. The
rows of this matrix contain the propagated (updated) streams for
period currentPeriod
, given the argument lag1Streams
matrix and the argument lag1LogWeights
vector for
currentPeriod  1
.
currentLogWeights
the propagated (updated) log weights
vector of length nStreamsToGenerate
, associated with the
streams in the rows of the returned currentStreams
matrix.
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
.
TRUE
or FALSE
reflecting the
decision of the resampling scheme implemented by this function.
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
.
see the subsection Arguments: for section Optional function: resampCriterionFunc.
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.
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.
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
.
a vector of length of dimSummPerPeriod
of summaries for currentPeriod
given the
currentStreams
and the currentLogWeights
.
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
.
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) MHupdated 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.
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.
Using very small values (<= 1e3
) for nStreams
might not give reliable results.
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
unlike,
particleFilter
, there is no builtin Metropolis
Hastings updating function, which generates proposals for
currentPeriod
streams using those of currentPeriod 
1
. The user needs to implement this function if
nMHSteps
> 0.
nStreamsPreResamp
it is set to nStreams
.
Also, the following point is worth noting:
resampCriterionFunc
, resampFunc
,
summaryFunc
are only necessary when user wants to try out new resampling schemes or enhanced summary generation procedures, 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
.
Gopi Goswami goswami@stat.harvard.edu
Jun S. Liu (2001). Monte Carlo strategies for scientific computing. Springer. Chapter 3.
Jun S. Liu and Rong Chen (1998). Sequential Monte Carlo methods for dynamical systems. Journal of the American Statistical Association 98(443): 10321044.
particleFilter
,
auxiliaryParticleFilter
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65  MSObj < MarkovSwitchingFuncGenerator(12345)
smcObj <
with(MSObj,
{
sequentialMonteCarlo(nStreams = 5000,
nPeriods = nrow(yy),
dimPerPeriod = ncol(yy),
propagateFunc = propagateFunc,
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(54321)
smcObj <
with(MSObj,
{
sequentialMonteCarlo(nStreams = 5000,
nPeriods = nrow(yy),
dimPerPeriod = ncol(yy),
propagateFunc = propagateFunc,
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)
})

##
## 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) 20062019 Gopi Goswami
##
## Created by: Gopi Goswami <goswami@stat.harvard.edu>
## Maintained by: Gopi Goswami <grgoswami@gmail.com>
##
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 sequentialMonteCarlo run:
nStreams: 5000
nStreamsPreResamp: 5000
nPeriods: 100
dimPerPeriod: 1
The summary of the resampling proportions:
[Note: resampling was done for 11 periods]
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
0.0920 0.2510 0.3152 0.3104 0.4170 0.4556 89
[1] "draws" "nStreams" "nPeriods"
[4] "dimPerPeriod" "nStreamsPreResamp" "nMHSteps"
[7] "summaryFunc" "time" "propagateFunc"
[10] "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: 0 secs (usr), 0 secs (sys), this period: 100]
E N D: SMC
The settings for this sequentialMonteCarlo run:
nStreams: 5000
nStreamsPreResamp: 5000
nPeriods: 100
dimPerPeriod: 1
The summary of the resampling proportions:
[Note: resampling was done for 16 periods]
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
0.0636 0.1475 0.2528 0.2329 0.3118 0.4294 84
[1] "draws" "nStreams" "nPeriods"
[4] "dimPerPeriod" "nStreamsPreResamp" "nMHSteps"
[7] "summaryFunc" "time" "propagateFunc"
[10] "filterType"
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.