IMIS: Incremental Mixture Importance Sampling

In IMIS: Increamental Mixture Importance Sampling

Description

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 <[e-mail]>

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")

Example output

Loading required package: mvtnorm
[1] "5000 likelihoods are evaluated in 0 minutes"
[1] "Stage   MargLike   UniquePoint   MaxWeight   ESS"
[1]    1.000   -3.435 1458.273    0.001 1532.030
[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.426 2448.024    0.000 7060.457
[1] "5000 likelihoods are evaluated in 0 minutes"
[1] "Stage   MargLike   UniquePoint   MaxWeight   ESS"
[1]    1.000   -0.816 1804.709    0.001 2454.909
[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= -2.18 , likelihood= 0.36 , prior= -2.54 , 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= -5.34 , likelihood= -2.29 , prior= -3.05 , 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= -8.88 , likelihood= -5.13 , prior= -3.75 , time used= 0 minutes, convergence= 0"
[1]    2.000   -0.812 2322.303    0.000 5339.194

IMIS documentation built on May 2, 2019, 4:05 p.m.