Learn R Programming

optismixture (version 0.1)

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

Description

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

Usage

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.p
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 f
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,
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
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
cv
TRUE indicates optimizing alpha and beta at the same time, and estimate with the formula $$\hat{\mu}_{\alpha_{**},\beta} = \frac{1}{n_2}\sum\limits_{i = 1}^{n_2}\frac{f(x_{i})p(x_{i}) - {\beta}^{\mathsf{T}}\bigl(\boldsymbol{q}(x_{i}) - p
opt.info
logical value indicating whether to save the returned value of the optimization procedure. See penoptpersp and penoptpersp.alpha.only for
opt.param
a list of paramters for the damped Newton method with backtracking line search [object Object],[object Object],[object Object] Only need to supply part of the list to change the default value.

Value

  • a list of [object Object],[object Object],[object Object],[object Object],[object Object],[object Object],[object Object]

Details

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

References

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

Examples

Run this code
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)[1], each = length(j.vec)),
                    variance = rep(j.vec, 2))
J <- dim(j.mat)[1] + 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")

Run the code above in your browser using DataLab