Nothing
### dnvmix() ###################################################################
##' @title Merge two matrices that are sorted in one of their columns
##' @param A (n, k) matrix, column 'sort.by' sorted increasingly
##' @param B (m, k) matrix, column 'sort.by' sorted increasingly
##' @param sort.by column which is sorted and should be sorted by
##' @return (n+m, k) matrix with rows from A and B so that
##' column 'sort.by' is sorted increasingly
##' @author Erik Hintz
##' @note k = 1 is allowed
merge_m <- function(A, B, sort.by = 1){
## Checks
if(!is.matrix(A)) A <- as.matrix(A)
if(!is.matrix(B)) B <- as.matrix(B)
dimA <- dim(A)
dimB <- dim(B)
stopifnot(dimA[2] == dimB[2]) # A and B have the same number of columns
n <- dimA[1] # number of rows of A
m <- dimB[1] # number of rows of B
res <- matrix(NA, ncol = dimA[2], nrow = n + m) # result matrix
p.A <- 1 # pointer to current element in 'A'
p.B <- 1 # pointer to current element in 'B'
for(i in 1:(n+m)){
res[i, ] <- if(A[p.A, sort.by] < B[p.B, sort.by]){
p.A <- p.A + 1
A[p.A-1, , drop = FALSE]
} else {
p.B <- p.B + 1
B[p.B-1, , drop = FALSE]
}
## Arrived at the end of 'A'
if(p.A == n + 1){
res[(i+1):(n+m), ] <- B[p.B:m, , drop = FALSE]
break
}
## Arrived at the end of 'B'
if(p.B == m + 1){
res[(i+1):(n+m), ] <- A[p.A:n, , drop = FALSE]
break
}
}
res
}
### dnvmix() ###################################################################
##' @title Merge two matrices that are sorted in one of their columns
##' @param A (n, k) matrix, column 'sort.by' sorted increasingly
##' @param B (m, k) matrix, nothing sorted
##' @param sort.by column which is sorted and should be sorted by
##' @return (n+m, k) matrix with rows from A and B so that
##' column 'sort.by' is sorted increasingly
##' @author Erik Hintz
##' @note k = 1 is allowed
merge_m2 <- function(A, B, sort.by = 1){
## Checks
if(!is.matrix(A)) A <- as.matrix(A)
if(!is.matrix(B)) B <- as.matrix(B)
dimA <- dim(A)
dimB <- dim(B)
stopifnot(dimA[2] == dimB[2]) # A and B have the same number of columns
n <- dimA[1] # number of rows of A
m <- dimB[1] # number of rows of B
res <- matrix(NA, ncol = dimA[2], nrow = n + m) # result matrix
p.A <- 1 # pointer to current element in 'A'
p.B <- 1 # pointer to current element in 'B'
for(i in 1:(n+m)){
res[i, ] <- if(A[p.A, sort.by] < B[p.B, sort.by]){
p.A <- p.A + 1
A[p.A-1, , drop = FALSE]
} else {
p.B <- p.B + 1
B[p.B-1, , drop = FALSE]
}
## Arrived at the end of 'A'
if(p.A == n + 1){
res[(i+1):(n+m), ] <- B[p.B:m, , drop = FALSE]
break
}
## Arrived at the end of 'B'
if(p.B == m + 1){
res[(i+1):(n+m), ] <- A[p.A:n, , drop = FALSE]
break
}
}
res
}
##' @title Exp - log trick for log(sum_i exp(a_i))
##' @param M (n1, n2) matrix
##' @return n2-vector log(colSums(exp(M)))
##' @author Erik Hintz
##' @note NO checking is done for efficiency reasons
logsumexp <- function(M)
{
n.rows <- nrow(M)
n.cols <- ncol(M)
if(!is.matrix(M)) M <- rbind(M)
cmax <- apply(M, 2, max)
cmax + log(.colSums(exp( M - rep(cmax, each = n.rows)), m = n.rows, n = n.cols))
}
##' @title Adaptive RQMC Method to estimate the log-density of an NVMIX
##' distribution
##' @param qW function of one variable specifying the quantile function of W.
##' @param maha2.2 squared maha-distances divided by 2. Assumed to be *sorted*
##' @param lconst see ?densmix_
##' @param d dimension of the underlying NVM_d distribution
##' @param k k = d for estimating log-density; k = d+2 for log-weights for fitnvmix
##' @param control see ?dnvmix
##' @param UsWs matrix; each row is of the form (u, qW(u)) where u in (0,1)
##' @return list of three ($ldensities, $error, $numiter); each vector of
##' length(maha2.2)
##' @author Erik Hintz and Marius Hofert
densmix_adaptrqmc <- function(qW, maha2.2, lconst, d, k = d, control, UsWs)
{
## 1 Set-up #################################################################
numiter <- 0 # counter for the number of iterations
ONE <- 1-.Machine$double.neg.eps
ZERO <- .Machine$double.neg.eps
tol.bisec <- control$dnvmix.tol.bisec # vector of length 3
## [1]/[2]/[3] => tolerance on 'u' / 'W' / 'log-integrand' for bisections
## Initialize output object
n <- length(maha2.2)
ldensities <- rep(NA, n)
error <- rep(NA, n)
numiters <- rep(NA, n)
## Grab 'UsWs'; store them in vectors and sort
stopifnot(is.matrix(UsWs), dim(UsWs)[2] == 2) # sanity check
ordering.U <- order(UsWs[, 1, drop = FALSE])
numObs <- dim(UsWs)[1] + 2 # will add 'ZERO' and 'ONE'
## Set up matrix of the form U | qW(U) | l.integrand(qW((U))
U.W.lint <- matrix(NA, ncol = 3, nrow = numObs)
U.W.lint[, 1] <- c(ZERO, UsWs[ordering.U, 1], ONE)
U.W.lint[, 2] <- c(max(qW(ZERO), ZERO), UsWs[ordering.U, 2], qW(ONE))
## Check if W is bounded
W.max <- qW(1) # <= inf
W.min <- qW(0) # >= 0
W.min_gen <- U.W.lint[1, 2] # qW(ZERO) >= 0; smallest 'W' we can generate
W.max_gen <- U.W.lint[numObs, 2] # qW(ONE) <= Inf; largest 'W' we can generate
isWfin <- is.finite(W.max)
## 2 Main loop over Mahalanobis distances ##################################
## For *each* maha2.2 value separately, find optimal integration region
## around the 'peek', integrate there via RQMC and outside via crude
## trapezoidal rules
for(ind in 1:n) {
curr.maha2.2 <- max(maha2.2[ind], ZERO) # avoid maha2.2 = 0
curr.lconst <- lconst[ind]
## Initialize various quantities
curr.error <- NA
ldens.right <- NA # log-density in (u_r, 1)
ldens.left <- NA # log-density in (0, u_l)
ldens.stratum <- NA # log-density in (u_l, u_r)
uLuR <- rep(NA, 2) # c('u.left', 'u.right') for later
## Determine location of maximum
int.argmax.u <- NA # u* such that qmix(u*) = W* where integrand(u*) is max
int.argmax.w <- 2*curr.maha2.2/k # value of W = qmix(u*) at max assuming W is *unbounded*
## Check special cases
peekRight <- FALSE # logical if peek beyond ONE (can't sample there)
peekLeft <- FALSE # logical if peek below ZERO (can't sample there)
if(int.argmax.w <= W.min_gen) {
## Case 1: 'peek' too close to 0
int.argmax.w <- if(int.argmax.w <= W.min){
## Can only happen when W.min > 0, ie when support(W) != (0, Inf)
W.min
} else {
peekLeft <- TRUE
W.min_gen
}
uLuR[1] <- 0 # u.left = 0
int.argmax.u <- 0
ldens.left <- -Inf # (0, u_l) = (0, 0) => empty region here
} else if(int.argmax.w >= W.max_gen){
## Case 2: 'peek' too close to 1
int.argmax.w <- if(int.argmax.w >= W.max){
W.max # W has bounded support!
} else {
peekRight <- TRUE
W.max_gen
}
uLuR[2] <- 1 # u.right = 1
int.argmax.u <- 1
ldens.right <- -Inf # (u_r, 1) = (1, 1) => empty region here
}
outofreach.w <- peekLeft || peekRight
## Realizations of log-integrand
U.W.lint[,3] <- curr.lconst - log(U.W.lint[,2])*k/2 -
curr.maha2.2/U.W.lint[,2] # sorted according to ordering(U)!
## Value of the theoretical maximum of the log-integrand
l.int.theo.max <- curr.lconst - log(int.argmax.w)*k/2 -
curr.maha2.2/int.argmax.w
## Threshold: Will only use RQMC where l.integrand > l.tol.int.lower
l.tol.int.lower <-
min(max(log(control$dnvmix.tol.int.lower), l.int.theo.max -
control$dnvmix.order.lower*log(10)), 0)
## 2.1 Find u* = argmax_u{integrand(u)} ################################
## If !NA => Already set above
if(is.na(int.argmax.u)) {
## Find u* such that g(u*) ~= g_max where g is the integrand
## <=> find u* such that qW(U*) = int.argmax.w via binary search
## => find starting values (u1, u2): qW(u1) < int.argmax.w & qW(u2) >= int.argmax.w
ind.first.bigger <- which(U.W.lint[, 2] >= int.argmax.w)[1]
u1u2 <- c(U.W.lint[ind.first.bigger - 1, 1],
U.W.lint[ind.first.bigger, 1])
## Set up bisection
convd <- FALSE
numiter <- 0
## Matrix to store additional u's and W's generated:
addvalues <- matrix(NA, ncol = 2, nrow = control$dnvmix.max.iter.bisec)
index.first.u <- which(U.W.lint[,1] >= u1u2[1])[1] # will insert after this index
while(!convd && numiter < control$dnvmix.max.iter.bisec) {
numiter <- numiter + 1
## Next point to check (midpoint of u1u2):
u.next <- mean(u1u2)
diff <- ((w.next <- max(ZERO, qW(u.next))) - int.argmax.w)
## Store u.next and quanitle:
addvalues[numiter, ] <- c(u.next, w.next)
## Update u1u2 depending on sign of 'diff' and check convergence:
if(diff > 0) u1u2[2] <- u.next else u1u2[1] <- u.next
convd <- (abs(diff) < tol.bisec[2]) || (diff(u1u2) < tol.bisec[1])
}
int.argmax.u <- u.next
## Add additional values of 'U', 'W' and 'l.integrand' in 'U.W.lint'
## while *preserving the order*.
## Order of first column = order of all columns
order.add.us <- order(addvalues[1:numiter, 1, drop = FALSE])
WsToAdd <- addvalues[order.add.us, 2, drop = FALSE] # needed twice
U.W.lint <-
rbind(U.W.lint[1:index.first.u,],
cbind(addvalues[order.add.us, 1, drop = FALSE], WsToAdd,
curr.lconst - log(WsToAdd)*k/2 - curr.maha2.2/WsToAdd),
U.W.lint[(index.first.u+1):numObs,])
## Update length (=nrow) of 'U.W.lint'
numObs <- numObs + numiter
}
## 2.2 Find stratum (u.left, u.right) over which RQMC is applied #######
## 2.2.1 Find "candidates" for bisection ###############################
## Need to distinguish multiple special cases:
greater.thshold <- (U.W.lint[, 3] > l.tol.int.lower)
if(any(greater.thshold)) {
## In this case there is at least one obs exceeding the threshold
ind.greater <- which(greater.thshold)
candid.left <- c(max(0, U.W.lint[ind.greater[1] - 1, 1]), # < thshold
min(int.argmax.u, U.W.lint[ind.greater[1], 1]) ) # > thshold
## Note: If ind.greater[1] = 1 => U[ind.greater[1] - 1] = numeric(0)
## => max(0, U[ind.greater[1] - 1]) = 0
## Need several cases for 'candid.right'.
last.ind.greater <- ind.greater[length(ind.greater)] # last index: g > thshold
candid.right <- if(last.ind.greater < numObs) {
## Let u^ be largest u >= u* such that g(u^) > threshold (exists)
## Case 1:
## There is u' with u'>=u^ and g(u') < threshold
## => 'u.right' in (u^, u')
c(max(U.W.lint[last.ind.greater, 1], int.argmax.u),
U.W.lint[last.ind.greater + 1, 1] )
} else {
## Case 2: the last point is such that g(u) > threshold => g(ONE) > threshold
## => u.right = 1 no matter where the peek is, but the handling of the region
## (ONE, 1) *does* depend on where the peek is.
ldens.right <- if(peekRight) {
## Area with peek right of ONE => can't simulate there
## => stratify from u.left to u.right = 1 (=ONE)
## => use *crude* approximation for region (ONE, 1) where max is
log(ZERO) + l.int.theo.max - log(2)
} else {
## Area with the peek is left of ONE but theoretical 'u.right' is
## out of reach due to machine precision.
## => stratify from u.left to u.right = 1 (=ONE)
## => do nothing in (ONE, 1)
-Inf
}
## See above: u.right = 1
c(1, 1)
}
} else if(peekRight) {
## In this case, none of the observations is > threshold and peek close to 1:
## Will use *all* obs to estimate in (0, ONE) and crude approximation in (ONE, 1)
## where max is. NO stratification!
candid.left <- c(1, 1)
candid.right <- c(1, 1)
ldens.right <- log(.Machine$double.neg.eps) + l.int.theo.max - log(2)
ldens.stratum <- -Inf
} else if(peekLeft) {
candid.left <- c(0, 0)
candid.right <- c(0, 0)
ldens.left <- log(.Machine$double.neg.eps) + l.int.theo.max - log(2)
ldens.stratum <- -Inf
} else {
candid.left <- c(0, int.argmax.u)
candid.right <- c(int.argmax.u, 1)
}
## 2.2.2 Bisection to find 'u.left' and 'u.right' ######################
candids <- rbind(candid.left, candid.right)
for(i in 1:2) {
curr.candid <- candids[i, ]
if(!is.na(uLuR[i])) next # we already set u.left or u.right
uLuR[i] <- if(diff(curr.candid) <= tol.bisec[1]) {
curr.candid[2]
} else {
## Use bisection similar to the one used to find u*
convd <- FALSE
numiter <- 0
## Matrix to store additional u's, W's, l.integrand's generated:
addvalues <- matrix(NA, ncol = 3,
nrow = control$dnvmix.max.iter.bisec)
while(!convd && numiter < control$dnvmix.max.iter.bisec) {
numiter <- numiter + 1
## Next point to check
u.next <- mean(curr.candid)
w.next <- max(qW(u.next), ZERO)
diff <- (l.int.next <- (curr.lconst -log(w.next)*k/2 -
curr.maha2.2/w.next)) - l.tol.int.lower
## Store values generated
addvalues[numiter, ] <- c(u.next, w.next, l.int.next)
## Update 'curr.candid' depending on sign of 'diff' and check convergence:
if(diff > 0) curr.candid[2] <- u.next else curr.candid[1] <- u.next
convd <-
(abs(diff) < tol.bisec[3]) || (diff(curr.candid) < tol.bisec[1])
}
## Update U.W.lint[]
U.W.lint <- rbind(U.W.lint, addvalues[1:numiter,, drop = FALSE])
## Sort again
U.W.lint <- U.W.lint[order(U.W.lint[, 1]), , drop = FALSE]
## Update length (=nrow) of 'U.W.lint'
numObs <- numObs + numiter
u.next
}
}
## For readability
u.left <- if(uLuR[1] <= ZERO) 0 else uLuR[1]
u.right <- if(uLuR[2] >= ONE) 1 else uLuR[2]
## 2.3 Estimate integral in the three regions ##########################
## 2.3.1 Estimate log of integral over (0, u.left)
ldens.left <- if(is.na(ldens.left)) {
## 'ldens.left' not set yet => peek is *not* in (0, u.left) => Riemann sums
if(u.left == 0) {
-Inf
} else if(u.left == 1) {
## => Special case where no obs > threshold
## => Use all obs and Riemann for this region
weights <- abs(c(U.W.lint[1, 1], U.W.lint[2:numObs, 1] -
U.W.lint[1:(numObs-1), 1]))
logsumexp(as.matrix(log(weights) +
(c(-Inf, U.W.lint[1:(numObs-1), 3]) +
U.W.lint[1:numObs, 3])/2, ncol = 1))
} else {
## 0 < u.left < 1 => Find obs in (0, u.left)
u_sml <- (U.W.lint[, 1] <= u.left)
sum_u_sml <- sum(u_sml)
if(sum_u_sml > 1) {
## Case 1: We have >1 observations in (0, u.left)
last_sml <- which(u_sml)[sum_u_sml]
weights <- abs(c(U.W.lint[1, 1],
U.W.lint[2:last_sml, 1] -
U.W.lint[1:(last_sml-1), 1]))
res <- logsumexp(as.matrix(log(weights) +
(c(-Inf, U.W.lint[1:(last_sml-1), 3]) +
U.W.lint[1:last_sml, 3])/2, ncol = 1))
if(is.nan(res)) -Inf else res
} else {
## Case 2: No observations in (0, u.left)
log(u.left) + l.tol.int.lower - log(2)
}
}
} else {
ldens.left
}
## 2.3.2 Estimate log of integral over (u.right, 1)
ldens.right <- if(is.na(ldens.right)) {
## 'ldens.right' not set yet => peek is *not* in (u.right, 1) => Riemann sums
if(u.right == 1) {
-Inf
} else if(u.right == 0) {
## => Special case where no obs > threshold
## => Use all obs and Riemann for this region
weights <-
abs(c(U.W.lint[1, 1], U.W.lint[2:numObs, 1]-U.W.lint[1:(numObs-1), 1]))
logsumexp(as.matrix(log(weights) +
(c(-Inf, U.W.lint[1:(numObs-1), 3]) +
U.W.lint[1:numObs, 3])/2, ncol = 1))
} else {
## 0 < u.right < 1 => Find obs in (u.right, 1)
u_gtr <- (U.W.lint[,1] >= u.right)
sum_u_gtr <- sum(u_gtr)
if(sum_u_gtr > 1) {
## Case 1: We have >1 observations in (u.right, 1)
first_gtr <- which(u_gtr)[1]
weights <- abs(c(U.W.lint[ (first_gtr+1):numObs, 1] -
U.W.lint[first_gtr:(numObs-1), 1],
.Machine$double.neg.eps))
res <- logsumexp(as.matrix(log(weights) +
(c(U.W.lint[(first_gtr+1):numObs, 3], -Inf) +
U.W.lint[first_gtr:numObs, 3])/2, ncol = 1))
if(is.nan(res)) -Inf else res
} else {
## Case 2: No or only one observations in (u.right, 1)
log1p(-u.right) + l.tol.int.lower - log(2)
}
}
} else {
ldens.right
}
## 2.3.3 Estimate log of integral over (u.left, u.right)
## We use RQMC in this region, unless 'ldens.strat' was already set.
stratlength <- u.right - u.left
rqmc.numiter <- 0
ldens.stratum <- if(is.na(ldens.stratum)) {
if(stratlength > control$dnvmix.tol.stratlength) {
ldens.obj <-
densmix_rqmc(qW, maha2.2 = curr.maha2.2, lconst = curr.lconst,
d = d, k = k, control = control, u.left = u.left,
u.right = u.right, return.all = FALSE,
max.iter.rqmc = control$max.iter.rqmc -
control$dnvmix.max.iter.rqmc.pilot)
curr.error <- ldens.obj$error
rqmc.numiter <- ldens.obj$numiter
ldens.obj$ldensities + log(stratlength)
} else {
log(max(stratlength, ZERO)) + l.int.theo.max - log(2)
}
} else {
ldens.stratum
}
## 2.4 Combine ###########################################################
ldensities[ind] <-
logsumexp(rbind(ldens.left, ldens.right, ldens.stratum, deparse.level = 0))
error[ind] <- curr.error
numiters[ind] <- rqmc.numiter
}
## 3. Return ################################################################
list(ldensities = ldensities, error = error, numiter = numiters)
}
##' @title Density of a Multivariate Normal Variance Mixture for restricted W
##' @param qW function of one variable specifying the quantile function of W.
##' @param maha2.2 squared maha-distances divided by 2. Assumed to be *sorted*
##' @param lrdet log(sqrt(det(scale))) where 'scale' is the scale matrix of
##' the normal variance mixture distribution.
##' @param d dimension of the Normal Variance Mixture
##' @param k power of qW^{-1} in the root (k=d => density, k=d+2 =>unnormalized
##' conditional expectation of 1/W)
##' @param control see ?dnvmix
##' @param u.left numeric in (0,1)
##' @param u.right numeric in (0,1), > u.left. Density will be estimated
##' conditional on W being between its 'u.left' and 'u.right' quantile.
##' @param max.iter.rqmc maximum number of iterations
##' @param return.all logical; if true, matrix (U, qW(U)) also returned.
##' @return List of three:
##' $ldensities n-vector with computed log-density values
##' $numiter numeric, number of iterations needed
##' $error n-vector of error estimates for log-densities; either
##' relative error or absolte error depending on is.na(control$dnvmix.reltol)
##' $UsWs (B, n) matrix (U, qW(U)) where U are uniforms
##' (only if return.all = TRUE)
##' @author Erik Hintz and Marius Hofert
densmix_rqmc <- function(qW, maha2.2, lconst, d, k = d, control, u.left = 0,
u.right = 1, max.iter.rqmc, return.all)
{
## 1. Setup #################################################################
dblng <- (control$increment == "doubling")
B <- control$B # number of randomizations
n <- length(maha2.2) # sample size
current.n <- control$fun.eval[1] # initial sample size
numiter <- 0 # counter for the number of iterations
total.fun.evals <- 0
ZERO <- .Machine$double.neg.eps
## Absolte/relative precision?
if(is.na(control$dnvmix.reltol)) {
## Use absolute error
tol <- control$dnvmix.abstol
do.reltol <- FALSE
} else {
## Use relative error
tol <- control$dnvmix.reltol
do.reltol <- TRUE
}
## Store seed if 'sobol' is used to get the same shifts later:
if(control$method == "sobol") {
seeds_ <- sample(1:(1e5*B), B) # B seeds for 'sobol()'
}
## Additional variables needed if the increment chosen is "dblng"
if(dblng) {
if(control$method == "sobol") useskip <- 0
denom <- 1
}
## Matrix to store RQMC estimates
rqmc.estimates <- matrix(-Inf, ncol = n, nrow = B)
## Will be needed a lot:
CI.factor.sqrt.B <- control$CI.factor / sqrt(B)
## Define trafo-function that maps u to (q,1) or (1,q) depending on 'up'
trafo <- function(u) u.left + (u.right - u.left)*u
## Initialize 'max.error' to > tol so that we can enter the while loop
max.error <- tol + 42
## Matrix to store U, W values => nrows = maximal number of funevals
if(return.all) {
max.nrow <- if(dblng) current.n*B*2^(max.iter.rqmc-1) else
max.iter.rqmc*B*current.n
UsWs <- matrix(NA, ncol = 2, nrow = max.nrow)
curr.lastrow <- 0 # will count row-index additional points are being inserted after
}
## 2. Main loop #############################################################
## while() runs until precision abstol is reached or the number of function
## evaluations exceed fun.eval[2]. In each iteration, B RQMC estimates of
## the desired log-densities are calculated.
while(max.error > tol & numiter < max.iter.rqmc &
total.fun.evals < control$fun.eval[2])
{
## In each randomization ...
for(b in 1:B) {
## Get the point set
U <- sort(switch(
control$method,
"sobol" = {
if(dblng) {
qrng::sobol(n = current.n, d = 1, randomize = "digital.shift",
seed = seeds_[b], skip = (useskip * current.n))
} else {
qrng::sobol(n = current.n, d = 1, randomize = "digital.shift",
seed = seeds_[b], skip = (numiter * current.n))
}
},
"ghalton" = {
qrng::ghalton(n = current.n, d = 1, method = "generalized")
},
"PRNG" = {
runif(current.n)
})) # sorted for later!
## Evaluate the integrand at the (next) point set
W <- qW(U <- trafo(U)) # realizations of the mixing variable; sorted!
## Need to replace values < ZERO by ZERO. W is *sorted*, so check using
## loop instd of 'pmax' (more efficient)
for(ind in 1:current.n) if(W[ind] < ZERO) W[ind] <- ZERO else break
## Update 'UsWs'
if(return.all) {
UsWs[(curr.lastrow + 1) : (curr.lastrow + current.n), ] <- cbind(U, W)
curr.lastrow <- curr.lastrow + current.n
}
next.estimate <- .Call("eval_densmix_integrand",
W = as.double(W),
maha2_2 = as.double(maha2.2),
current_n = as.integer(current.n),
n = as.integer(n),
d = as.integer(d),
k = as.integer(k),
lconst = as.double(lconst))
## Update RQMC estimates
rqmc.estimates[b,] <-
if(dblng) {
## In this case both, rqmc.estimates[b,] and
## next.estimate depend on n.current points
.Call("logsumexp2",
a = as.double(rqmc.estimates[b, ]),
b = as.double(next.estimate),
n = as.integer(n)) - log(denom)
} else {
## In this case, rqmc.estimates[b,] depends on
## numiter * current.n points whereas next.estimate
## depends on current.n points
.Call("logsumexp2",
a = as.double(rqmc.estimates[b,] + log(numiter)),
b = as.double(next.estimate),
n = as.integer(n)) - log(numiter + 1)
}
} # end for(b in 1:B)
## Update of various variables
## Double sample size and adjust denominator in averaging as well as useskip
if(dblng) {
## Change denom and useksip (exactly once, in the first iteration)
if(numiter == 0) {
denom <- 2
useskip <- 1
} else {
## Increase sample size n. This is done in all iterations
## except for the first two
current.n <- 2 * current.n
}
}
## Total number of function evaluations
total.fun.evals <- total.fun.evals + B * current.n
numiter <- numiter + 1
## Update error. The following is slightly faster than 'apply(..., 2, var)'
ldensities <- logsumexp(rqmc.estimates) - log(B) # performs better than .colMeans
vars <- .colMeans((rqmc.estimates - rep(ldensities, each = B))^2, B, n, 0)
error <- if(!do.reltol) { # absolute error
sqrt(vars)*CI.factor.sqrt.B
} else { # relative error
sqrt(vars)/abs(ldensities)*CI.factor.sqrt.B
}
max.error <- max(error)
} # while()
## 3. Return ################################################################
if(return.all) {
list(ldensities = ldensities, numiter = numiter, error = error,
UsWs = UsWs[1:curr.lastrow,])
} else {
list(ldensities = ldensities, numiter = numiter, error = error)
}
}
##' @title Integration Routine of log h(u) on the unit interval
##' @param qW function of one variable specifying the quantile function of W.
##' @param maha2.2 squared maha-distances divided by 2. Assumed to be *sorted*
##' @param lconst vector of same length as 'maha2.2', see below under 'Note'
##' @param d dimension of the Normal Variance Mixture
##' @param control see ?get_set_param()
##' @param verbose see ?dnvmix()
##' @return List of three:
##' $ldensities n-vector with computed log-density values
##' $error n-vector of *absolute* error estimates for log-densities
##' $numiter n-vector of number of iterations needed
##' @author Erik Hintz and Marius Hofert
##' @note Estimates int_0^1 h(u) du where log h(u) = lconst - d/2 log qW(u) - maha2.2/qW(u).
##' For density of 'nvmix': lconst = rep(-lrdet -d/2log(2pi), length(maha2.2));
##' For density of 'gammamix': lconst = -lgamma(d/2) - d/2 log(2) + (d/2-1)log(2*maha2.2)
densmix_ <- function(qW, maha2.2, lconst, d, control, verbose)
{
## Absolte/relative precision?
tol <- if(is.na(control$dnvmix.reltol)) { # if 'reltol = NA' use absolute precision
do.reltol <- FALSE
control$dnvmix.abstol
} else { # otherwise use relative precision (default)
do.reltol <- TRUE
control$dnvmix.reltol
}
## Call RQMC procedure without any stratification
rqmc.obj <- densmix_rqmc(qW, maha2.2 = maha2.2, lconst = lconst, d = d,
control = control, u.left = 0, u.right = 1,
max.iter.rqmc = control$dnvmix.max.iter.rqmc.pilot,
return.all = TRUE)
## Extract results
ldens <- rqmc.obj$ldensities
numiter <- rep(rqmc.obj$numiter, length(maha2.2))
error <- rqmc.obj$error
if(any(error > tol)) {
## Accuracy not reached for at least one 'maha2.2' value
## => Use adaptive approach
if(control$dnvmix.doAdapt) {
notRchd <- which(error > tol)
rqmc.obj <- densmix_adaptrqmc(
qW, maha2.2 = maha2.2[notRchd], lconst = lconst[notRchd], d = d,
UsWs = rqmc.obj$UsWs,control = control)
ldens[notRchd] <- rqmc.obj$ldensities
numiter[notRchd] <- numiter[notRchd] + rqmc.obj$numiter
error[notRchd] <- rqmc.obj$error
}
## Handle warnings:
if(as.logical(verbose)) { # verbose can be passed by 'fitnvmix' => between 0 and 3
if(any(is.na(error))) {
## At least one error is NA
warning("Estimation unreliable, corresponding error estimate NA")
}
whichNA <- which(is.na(error))
if(any(error[setdiff(1:length(error), whichNA)] > tol)) # 'setdiff' needed if 'whichNA' is empty
warning("Tolerance not reached for all inputs; consider increasing 'max.iter.rqmc' in the 'control' argument.")
}
}
## Compute relative/absolute error for return
if(do.reltol){
relerror <- error
abserror <- error * abs(ldens)
} else {
abserror <- error
relerror <- error/abs(ldens)
}
## Return
list(ldensities = ldens, numiter = numiter, abserror = abserror,
relerror = relerror)
}
##' @title Density of a Multivariate Normal Variance Mixture
##' @param x (n, d)-matrix of evaluation points
##' @param qmix specification of the (mixture) distribution of W. This can be:
##' 1) a character string specifying a supported distribution (additional
##' arguments of this distribution are passed via '...').
##' 2) a list of length at least one; the first argument specifies
##' the base name of an existing distribution which can be sampled
##' with prefix "q", the other elements denote additional parameters
##' passed to this "rmix" random number generator.
##' 3) a function being interpreted as the quantile function F_W^-.
##' @param loc d-vector (location vector)
##' @param scale (d, d)-covariance matrix (scale matrix)
##' @param factor Cholesky factor (lower triangular matrix) of 'scale';
##' important here so that det(scale) is computed correctly!
##' @param control list; see ?get_set_param()
##' @param log logical indicating whether the logarithmic density is to be computed
##' @param verbose logical indicating whether warnings shall be thrown.
##' @param ... additional arguments passed to the underlying mixing distribution
##' @return n-vector with computed density values and attributes 'error'
##' (error estimate) and 'numiter' (number of while-loop iterations)
##' @author Erik Hintz and Marius Hofert
dnvmix <- function(x, qmix, loc = rep(0, d), scale = diag(d),
factor = NULL, # needs to be lower triangular!
control = list(), log = FALSE, verbose = TRUE, ...)
{
## 1 Setup ##################################################################
if(!is.matrix(x)) x <- rbind(x)
d <- ncol(x) # dimension
if(!is.matrix(scale)) scale <- as.matrix(scale)
stopifnot(length(loc) == d, dim(scale) == c(d, d))
verbose <- as.logical(verbose)
## Deal with algorithm parameters
control <- get_set_param(control)
## If 'factor' is not provided, determine it here as a *lower* triangular matrix
if(is.null(factor)) factor <- t(chol(scale)) # lower triangular
## Deal with 'qmix'
mix_list <- get_mix_(qmix = qmix, callingfun = "dnvmix", ... )
qW <- mix_list[[1]] # function(u)
special.mix <- mix_list[[2]] # string or NA
## Build result object (log-density)
lres <- rep(-Inf, (n <- nrow(x))) # n-vector of results
abserror <- rep(NA, n)
relerror <- rep(NA, n)
notNA <- rowSums(is.na(x)) == 0
lres[!notNA] <- NA
x <- x[notNA,, drop = FALSE] # non-missing data (rows)
## 2 Actual computation ####################################################
## Recall that 'factor' is *lower triangular*. For short, let 'factor' = L
## Solve L * z = x_i - mu for z, so z = L^{-1} * (x_i - mu) (d vector)
## => z^2 (=> componentwise) = z^T z = (x_i - mu)^T * (L^{-1})^T L^{-1} (x_i - mu)
## = z^T z = (x_i - mu)^T * (L L^T )^{-1} (x_i - mu)
## = (x_i - mu)^T * scale^{-1} * (x_i - mu)
## = quadratic form
## Now do this for *all* x_i simultaneously using that L is lower triangular:
## Forwardsolve: "right-hand sides" of equation must be in the *columns*, thus t(x)
z <- forwardsolve(factor, t(x) - loc, transpose = FALSE)
maha2 <- colSums(z^2) # = sum(z^T z); n-vector of squared Mahalanobis distances from x to mu w.r.t. scale
## Note: could probably be done with mahalanobis() but unclear how we would
## get det(scale) then.
## log(sqrt(det(scale))) = log(det(scale))/2 = log(det(R^T R))/2 = log(det(R)^2)/2
## = log(prod(diag(R))) = sum(log(diag(R)))
lrdet <- sum(log(diag(factor)))
if(!is.finite(lrdet)) stop(paste("Density not defined for singular 'scale' "))
## Counter
numiter <- 0 # initialize counter
## Deal with the different distributions
if(!is.na(special.mix)) {
lres[notNA] <- switch(
special.mix,
"inverse.gamma" = {
df <- mix_list$param
lgamma((df + d) / 2) - lgamma(df/2) - (d/2) * log(df*pi) - lrdet -
(df+d)/2 * log1p(maha2/df)
},
"constant" = {
-(d/2) * log(2 * pi) - lrdet - maha2/2
},
"pareto" = {
alpha <- mix_list$param
log(alpha) - d/2*log(2*pi) - lrdet - (alpha+d/2)*log(maha2/2) +
pgamma(maha2/2, scale = 1, shape = alpha+d/2, log.p = TRUE) +
lgamma(alpha+d/2)
})
if(!log) lres <- exp(lres) # already exponentiate
abserror <- rep(0, length(maha2))
relerror <- rep(0, length(maha2)) # no error here
} else {
## General case of a multivariate normal variance mixture (=> RQMC)
## Prepare inputs for densmix_()
## Sort maha-distance and divide by 2; store ordering to recover original
## ordering later
ordering.maha <- order(maha2)
maha2.2 <- maha2[ordering.maha]/2
## Define log-constant for the integration
lconst <- rep(-lrdet - d/2*log(2*pi) , length(maha2.2))
## Call internal densmix_ (which itself calls C-Code and handles warnings)
ests <- densmix_(qW, maha2.2 = maha2.2, lconst = lconst,
d = d, control = control, verbose = verbose)
## Grab results
lres[notNA] <- ests$ldensities[order(ordering.maha)]
relerror[notNA] <- ests$relerror[order(ordering.maha)]
abserror[notNA] <- ests$abserror[order(ordering.maha)]
## Correct results and error if 'log = FALSE'
if(!log){
lres[notNA] <- exp(lres[notNA])
## CI for mu: exp(logmu_hat +/- abserr(logmu_hat))) = (lower, upper)
## => compute max. error on mu_hat as max( (upper - mu), (mu - lower) )
relerror[notNA] <- max( (exp(abserror[notNA]) - 1), (1 - exp(-abserror[notNA])) )
abserror[notNA] <- lres[notNA] * relerror[notNA]
}
numiter <- ests$numiter
}
## 3. Return ################################################################
## Note that 'lres' was exponentiated already if necessary.
attr(lres, "abs. error") <- abserror
attr(lres, "rel. error") <- relerror
attr(lres, "numiter") <- numiter
lres
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.