IMIS algorithm draws samples from the posterior distribution of multivariate variable. The user has to define the following R functions in advance: prior(x) calculates prior density of x, likelihood(x) calculates the likelihood of x, and sample.prior(n) draws n samples from the prior distribution. For multivariate x, the prior or the likelihood of a vector should be a scalar, and the prior or the likelihood of a matrix should be a vector.
Usage
IMIS(B, B.re, number_k, D)
Arguments
B
The incremental sample size at each iteration of IMIS.
B.re
The desired posterior sample size at the resample stage.
number_k
The maximum number of iterations in IMIS.
D
The number of optimizers which could be 0.
Value
resample
The posterior resamples
stat
Diagnostic statistics at each IMIS iteration: Marginal likelihood (1st col), Expected number of unique points among the posterior resamples (2nd col), Maximum importance weight (3rd col), Effective sample size (4th col), Entropy of importance weights relative to the uniform distribution (5th col), Rescaled variance of importance weights (6th col).
center
The centers of Gaussian components
Author(s)
Adrian Raftery and Le Bao <lebao@uw.edu>
References
Raftery A.E. and Bao L. (2010) "Estimating and projecting trends in HIV/AIDS generalized epidemics using incremental mixture importance sampling." Biometrics.
Examples
## Example for multivariate case
likelihood <- function(theta) dmvnorm(theta, c(1,1), matrix(c(1,0.6,0.6,1),2,2))
prior <- function(theta) dmvnorm(theta, c(0,0), diag(3,2))
sample.prior <- function(n) rmvnorm(n, c(0,0), diag(3,2))
result = IMIS(500, 3000, 100, 10)
x1 = x2 = seq(-2, 4, by=0.1)
z = matrix(NA,length(x1),length(x2))
for (i in 1:length(x1))
for (j in 1:length(x2))
z[i,j] = likelihood(c(x1[i],x2[j])) * prior(c(x1[i],x2[j]))
contour(x1, x2, z, drawlabels=FALSE, pty="s")
points(result$resample[,1], result$resample[,2], cex=0.1)
## Example for univariate case
likelihood <- function(theta) exp(-1*sin(3*theta)*sin(theta^2) - 0.1*theta^2)
prior <- function(theta) dnorm(theta, 0, 5)
sample.prior <- function(n) rnorm(n, 0, 5)
result = IMIS(500, 3000, 100, 10)
plot(density(result$resample, adjust=0.3), xlim=c(-5,5), main = "wild function")
x = seq(-5, 5, 0.001)
lines(prior(x)*likelihood(x)~x, xlim=c(-5,5), col="red")
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(IMIS)
Loading required package: mvtnorm
> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/IMIS/IMIS.Rd_%03d_medium.png", width=480, height=480)
> ### Name: IMIS
> ### Title: Incremental Mixture Importance Sampling
> ### Aliases: IMIS
>
> ### ** Examples
>
> ## Example for multivariate case
> likelihood <- function(theta) dmvnorm(theta, c(1,1), matrix(c(1,0.6,0.6,1),2,2))
> prior <- function(theta) dmvnorm(theta, c(0,0), diag(3,2))
> sample.prior <- function(n) rmvnorm(n, c(0,0), diag(3,2))
> result = IMIS(500, 3000, 100, 10)
[1] "5000 likelihoods are evaluated in 0 minutes"
[1] "Stage MargLike UniquePoint MaxWeight ESS"
[1] 1.000 -3.426 1446.263 0.001 1509.554
[1] "maximum posterior= -4.77 , likelihood= -1.69 , prior= -3.08 , time used= 0 minutes, convergence= 0"
[1] "maximum posterior= -4.77 , likelihood= -1.69 , prior= -3.08 , time used= 0 minutes, convergence= 0"
[1] "maximum posterior= -4.77 , likelihood= -1.69 , prior= -3.08 , time used= 0 minutes, convergence= 0"
[1] "maximum posterior= -4.77 , likelihood= -1.69 , prior= -3.08 , time used= 0 minutes, convergence= 0"
[1] "maximum posterior= -4.77 , likelihood= -1.69 , prior= -3.08 , time used= 0 minutes, convergence= 0"
[1] "maximum posterior= -4.77 , likelihood= -1.69 , prior= -3.08 , time used= 0 minutes, convergence= 0"
[1] "maximum posterior= -4.77 , likelihood= -1.69 , prior= -3.08 , time used= 0 minutes, convergence= 0"
[1] "maximum posterior= -4.77 , likelihood= -1.69 , prior= -3.08 , time used= 0 minutes, convergence= 0"
[1] "maximum posterior= -4.77 , likelihood= -1.69 , prior= -3.08 , time used= 0 minutes, convergence= 0"
[1] "maximum posterior= -4.77 , likelihood= -1.69 , prior= -3.08 , time used= 0 minutes, convergence= 0"
[1] 2.000 -3.430 2448.009 0.000 7058.853
> x1 = x2 = seq(-2, 4, by=0.1)
> z = matrix(NA,length(x1),length(x2))
> for (i in 1:length(x1))
+ for (j in 1:length(x2))
+ z[i,j] = likelihood(c(x1[i],x2[j])) * prior(c(x1[i],x2[j]))
> contour(x1, x2, z, drawlabels=FALSE, pty="s")
> points(result$resample[,1], result$resample[,2], cex=0.1)
>
> ## Example for univariate case
> likelihood <- function(theta) exp(-1*sin(3*theta)*sin(theta^2) - 0.1*theta^2)
> prior <- function(theta) dnorm(theta, 0, 5)
> sample.prior <- function(n) rnorm(n, 0, 5)
> result = IMIS(500, 3000, 100, 10)
[1] "5000 likelihoods are evaluated in 0 minutes"
[1] "Stage MargLike UniquePoint MaxWeight ESS"
[1] 1.000 -0.840 1783.181 0.001 2382.389
[1] "maximum posterior= -1.96 , likelihood= 0.61 , prior= -2.57 , time used= 0 minutes, convergence= 0"
[1] "maximum posterior= -2.18 , likelihood= 0.36 , prior= -2.54 , time used= 0 minutes, convergence= 0"
[1] "maximum posterior= -2.55 , likelihood= 0.13 , prior= -2.68 , time used= 0 minutes, convergence= 0"
[1] "maximum posterior= -2.18 , likelihood= 0.36 , prior= -2.54 , time used= 0 minutes, convergence= 0"
[1] "maximum posterior= -2.68 , likelihood= -0.04 , prior= -2.63 , time used= 0 minutes, convergence= 0"
[1] "maximum posterior= -3.24 , likelihood= -0.43 , prior= -2.81 , time used= 0 minutes, convergence= 0"
[1] "maximum posterior= -2.53 , likelihood= 0 , prior= -2.53 , time used= 0 minutes, convergence= 0"
[1] "maximum posterior= -5.48 , likelihood= -2.3 , prior= -3.19 , time used= 0 minutes, convergence= 0"
[1] "maximum posterior= -5.47 , likelihood= -2.35 , prior= -3.13 , time used= 0 minutes, convergence= 0"
[1] "maximum posterior= -10.1 , likelihood= -6.23 , prior= -3.87 , time used= 0 minutes, convergence= 0"
[1] 2.000 -0.823 2338.445 0.000 5577.408
> plot(density(result$resample, adjust=0.3), xlim=c(-5,5), main = "wild function")
> x = seq(-5, 5, 0.001)
> lines(prior(x)*likelihood(x)~x, xlim=c(-5,5), col="red")
>
>
>
>
>
> dev.off()
null device
1
>
providing 
.
