R Graphical Manual

Browse All

Last data update: 2014.03.03

R: Importance Sampling of Exponential Family State Space Model

importanceSSM

R Documentation

Importance Sampling of Exponential Family State Space Model

Description

Function importanceSSM simulates states or signals of the exponential family state space model conditioned with the observations, returning the simulated samples of the states/signals with the corresponding importance weights.

Usage

importanceSSM(model, type = c("states", "signals"), filtered = FALSE,
  nsim = 1000, save.model = FALSE, theta, antithetics = FALSE,
  maxiter = 50)

Arguments

`model`	Exponential family state space model of class `SSModel`.
`type`	What to simulate, `"states"` or `"signals"`. Default is `"states"`
`filtered`	Simulate from p(α_t\|y_{t-1},...,y_1) instead of p(α\|y). Note that for large models this can be very slow. Default is FALSE.
`nsim`	Number of independent samples. Default is 1000.
`save.model`	Return the original model with the samples. Default is FALSE.
`theta`	Initial values for the conditional mode theta.
`antithetics`	Logical. If TRUE, two antithetic variables are used in simulations, one for location and another for scale. Default is FALSE.
`maxiter`	Maximum number of iterations used in linearisation. Default is 50.

Details

Function can use two antithetic variables, one for location and other for scale, so output contains four blocks of simulated values which correlate which each other (ith block correlates negatively with (i+1)th block, and positively with (i+2)th block etc.).

Value

A list containing elements

`samples`	Simulated samples.
`weights`	Importance weights.
`model`	Original model in case of `save.model==TRUE`.

Examples

data("sexratio")
model <- SSModel(Male ~ SSMtrend(1, Q = list(NA)), u = sexratio[,"Total"], data = sexratio,
                distribution = "binomial")
fit <- fitSSM(model, inits = -15, method = "BFGS")
fit$model$Q #1.107652e-06
# Computing confidence intervals for sex ratio
# Uses importance sampling on response scale (1000 samples with antithetics)
set.seed(1)
imp <- importanceSSM(fit$model, nsim = 250, antithetics = TRUE)
sexratio.smooth <- numeric(length(model$y))
sexratio.ci <- matrix(0, length(model$y), 2)
w <- imp$w/sum(imp$w)
for(i in 1:length(model$y)){
  sexr <- exp(imp$sample[i,1,])
  sexratio.smooth[i]<-sum(sexr*w)
  oo <- order(sexr)
  sexratio.ci[i,] <- c(sexr[oo][which.min(abs(cumsum(w[oo]) - 0.05))],
                   sexr[oo][which.min(abs(cumsum(w[oo]) - 0.95))])
}

## Not run: 
# Filtered estimates
impf <- importanceSSM(fit$model, nsim = 250, antithetics = TRUE,filtered=TRUE)
sexratio.filter <- rep(NA,length(model$y))
sexratio.fci <- matrix(NA, length(model$y), 2)
w <- impf$w/rowSums(impf$w)
for(i in 2:length(model$y)){
  sexr <- exp(impf$sample[i,1,])
  sexratio.filter[i] <- sum(sexr*w[i,])
  oo<-order(sexr)
  sexratio.fci[i,] <- c(sexr[oo][which.min(abs(cumsum(w[i,oo]) - 0.05))],
                    sexr[oo][which.min(abs(cumsum(w[i,oo]) - 0.95))])
}

ts.plot(cbind(sexratio.smooth,sexratio.ci,sexratio.filter,sexratio.fci),
        col=c(1,1,1,2,2,2),lty=c(1,2,2,1,2,2))

## End(Not run)

Results


R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(KFAS)
> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/KFAS/importanceSSM.Rd_%03d_medium.png", width=480, height=480)
> ### Name: importanceSSM
> ### Title: Importance Sampling of Exponential Family State Space Model
> ### Aliases: importanceSSM
> 
> ### ** Examples
> 
> data("sexratio")
> model <- SSModel(Male ~ SSMtrend(1, Q = list(NA)), u = sexratio[,"Total"], data = sexratio,
+                 distribution = "binomial")
> fit <- fitSSM(model, inits = -15, method = "BFGS")
> fit$model$Q #1.107652e-06
, , 1

             [,1]
[1,] 1.107654e-06

> # Computing confidence intervals for sex ratio
> # Uses importance sampling on response scale (1000 samples with antithetics)
> set.seed(1)
> imp <- importanceSSM(fit$model, nsim = 250, antithetics = TRUE)
> sexratio.smooth <- numeric(length(model$y))
> sexratio.ci <- matrix(0, length(model$y), 2)
> w <- imp$w/sum(imp$w)
> for(i in 1:length(model$y)){
+   sexr <- exp(imp$sample[i,1,])
+   sexratio.smooth[i]<-sum(sexr*w)
+   oo <- order(sexr)
+   sexratio.ci[i,] <- c(sexr[oo][which.min(abs(cumsum(w[oo]) - 0.05))],
+                    sexr[oo][which.min(abs(cumsum(w[oo]) - 0.95))])
+ }
> 
> ## Not run: 
> ##D # Filtered estimates
> ##D impf <- importanceSSM(fit$model, nsim = 250, antithetics = TRUE,filtered=TRUE)
> ##D sexratio.filter <- rep(NA,length(model$y))
> ##D sexratio.fci <- matrix(NA, length(model$y), 2)
> ##D w <- impf$w/rowSums(impf$w)
> ##D for(i in 2:length(model$y)){
> ##D   sexr <- exp(impf$sample[i,1,])
> ##D   sexratio.filter[i] <- sum(sexr*w[i,])
> ##D   oo<-order(sexr)
> ##D   sexratio.fci[i,] <- c(sexr[oo][which.min(abs(cumsum(w[i,oo]) - 0.05))],
> ##D                     sexr[oo][which.min(abs(cumsum(w[i,oo]) - 0.95))])
> ##D }
> ##D 
> ##D ts.plot(cbind(sexratio.smooth,sexratio.ci,sexratio.filter,sexratio.fci),
> ##D         col=c(1,1,1,2,2,2),lty=c(1,2,2,1,2,2))
> ## End(Not run)
> 
> 
> 
> 
> 
> dev.off()
null device 
          1 
>