# batch.estimation: Two stage estimation, a pilot estimate of mixing alpha and a... In optismixture: Optimal Mixture Weights in Multiple Importance Sampling

## Description

Two stage estimation, a pilot estimate of mixing alpha and a following importance sampling, with or without control variates

## Usage

 1 2 3 4 5 batch.estimation(seed, batch.size, mixture.param, eps = rep(0.1/mixture.param$J, mixture.param$J), fname = "f", rpname = "rp", rqname = "rq", dpname = "dp", dqname = "dq", cv = TRUE, opt.info = FALSE, opt.param = list(reltol = 0.001, relerr = 0.001, rho0 = 1, maxin = 20, maxout = 30)) 

## Arguments

 seed seed for sampling batch.size length two vector of batch sizes mixture.param mixture.param = list(p, J, ...), where p is the dimension of the sample, and J is the number of mixture components, including the defensive one. mixture.param should be compatible with user defined functions f(n, j, mixture.param), rp(n, mixture.param), rq(n, j, mixture.param), dp(xmat, mixture.param), dq(xmat, j, mixture.param) eps the lower bound for optimizing alpha. if eps is of length 1, it is expanded to rep(eps, J), default to be rep(0.1/J, J) fname name of user defined function fname(xmat, j, mixture.param). xmat is an n \times p matrix of n samples with p dimensions. fname returns a vector of function values for each row in xmat. fname is defined for j = 1, \cdots, J. j = 1, \cdots, J - 1 corresponds to different proposal mixture components, and j = J corresponds to the defensive mixture component. rpname name of user definded function rpname(n, mixture.param). It generates n random samples from target distribution pname. Parameters can be specified in mixture.param. rpname returns an n \times p matrix. rqname name of user defined function rqname(n, j, mixture.param). It generate n random samples from the j^{th} mixture component of proposal mixture distribution. rqname returns an n \times p matrix. rqname is defined for j = 1, \cdots, J - 1. dpname name of user defined function dpname(xmat, mixture.param). It returns the density of xmat from the target distribution pname as a vector of length nrow(xmat). Note that only the ratio between dpname and dqname matters. So dpname and dqname can return values of C \timesdpname and C \timesdqname respectively. dqname name of user defined function dqname(xmat, j, mixture.param). It returns the density of xmat from the proposal distribution q_j as a vector of length nrow(xmat). dqname is defined for j = 1, \cdots, J - 1. Note that only the ratio between dpname and dqname matters. So dpname and dqname can return values of C \timesdpname and C \timesdqname respectively. cv TRUE indicates optimizing alpha and beta at the same time, and estimate with the formula \hat{μ}_{α_{**},β} = \frac{1}{n_2}∑\limits_{i = 1}^{n_2}\frac{f(x_{i})p(x_{i}) - {β}^{\mathsf{T}}\bigl(\boldsymbol{q}(x_{i}) - p(x_{i})\boldsymbol{1}\bigr)}{q_{α_{**}}(x_{i})}. FALSE indicates optimizing alpha only, and estimate with the formula \hat{μ}_{α_{*}}= \frac{1}{n_2}∑\limits_{i = 1}^{n_2}\frac{f(x_{i})p(x_{i})}{q_{α_{*}}(x_{i})}. opt.info logical value indicating whether to save the returned value of the optimization procedure. See penoptpersp and penoptpersp.alpha.only for the returned value. opt.param a list of paramters for the damped Newton method with backtracking line search reltolrelative tolerence in dampednewton step, default to be 10^{-2} relerrNewton step stop when within (1+relerr) of minimum variance, default to be 10^{-3} rho0initial value for rho, default to be 1 Only need to supply part of the list to change the default value.

## Details

Estimate E_p f with a two step importance sample procedure. See He & Owen(2014) for details.

## Value

a list of

mu.hat

the estimate for mu

sd.hat

estimated sd of mu.hat

alpha.opt

the estimated optimal alpha

beta.opt

if cv = TRUE, the estimated optimal beta

opt.info

if opt.info = TRUE, also return the list(x=x, y=y, z=z, alpha=alpha,beta=beta,rho=rho,f=f

rhopen=rhopen, outer=log(rho0/rho,mu), relerr = relerr, alphasum = sum(alpha))

from the optimization after batch 1.

## References

Boyd, S. and Vandeberghe, L. (2004). Convex optimization. Cambridge University Press, Cambridge

## Examples

  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 library(optismixture) seed <- 1 p <- 5 rho <- 1/2 gamma <- 2.4 sigma.dvar <- function(rho, p){ sigma <- matrix(0, p, p) for(i in 1:(p-1)){ for(j in (i+1):p){ sigma[i,j] <- rho^(abs(i-j)) } } sigma <- sigma + t(sigma) diag(sigma) <- 1 return(sigma) } sigma <- sigma.dvar(rho, p) batch.size <- c(10^4, 1002) j.vec <- 2^-(seq(1,5,1)) eps.status <- 1 eps.safe <- 0.1 ## initialization and construct derived parameters mu <- rep(0, p) x0 <- matrix(1, 1, p) x0.mat <- rbind(rep(1,p), rep(-1, p)) j.mat <- data.frame(centerid = rep(1:dim(x0.mat), each = length(j.vec)), variance = rep(j.vec, 2)) J <- dim(j.mat) + 1 eps <- rep(0.1/J, J) mixture.param <- list(x0 = x0, x0.mat = x0.mat, p = p, sigma = sigma, gamma = gamma, j.mat = j.mat, J = J) f <- function(x, j, mixture.param){ f1 <- function(x, mixture.param){ x0 <- mixture.param$x0 gamma <- mixture.param$gamma return(sum((x - x0)^2)^(-gamma/2)) } return(apply(x, 1, f1, mixture.param)) } dq <- function(x, j, mixture.param){ centerid <- mixture.param$j.mat[j, 1] j.param <- mixture.param$j.mat[j, 2] return(mvtnorm::dmvnorm(x, mixture.param$x0.mat[centerid,], j.param*diag(mixture.param$p))) } dp <- function(x, mixture.param){ return(mvtnorm::dmvnorm(x, rep(0, mixture.param$p), mixture.param$sigma)) } rq <- function(n, j, mixture.param){ centerid <- mixture.param$j.mat[j, 1] j.param <- mixture.param$j.mat[j,2] return(mvtnorm::rmvnorm(n, mixture.param$x0.mat[centerid, ], j.param*diag(mixture.param$p))) } rp <- function(n, mixture.param){ mu <- rep(0, mixture.param$p) sigma <- mixture.param$sigma return(mvtnorm::rmvnorm(n, mu, sigma)) } a <- batch.estimation(seed, batch.size, mixture.param, eps, cv = FALSE, fname = "f", rpname = "rp", rqname = "rq", dpname = "dp", dqname = "dq") 

optismixture documentation built on May 1, 2019, 10:13 p.m.