Description Usage Arguments Details Value References Examples
Two stage estimation, a pilot estimate of mixing alpha and a following importance sampling, with or without control variates
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))
|
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 |
opt.param |
a list of paramters for the damped Newton method with backtracking line search
Only need to supply part of the list to change the default value. |
Estimate E_p f with a two step importance sample procedure. See He & Owen(2014) for details.
a list of
the estimate for mu
estimated sd of mu.hat
the estimated optimal alpha
if cv = TRUE, the estimated optimal beta
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.
Boyd, S. and Vandeberghe, L. (2004). Convex optimization. Cambridge University Press, Cambridge
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)[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")
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.