#' @rdname dgmcm.loglik
#' @param res The resolution at which the inversion of \code{qgmm.marginal} is
#' done. Default is 1000.
#' @param spread The number of marginal standard deviations from the marginal
#' means the \code{pgmm.marginal} is to be evaluated on.
#' @param rule The extrapolation rule used in \code{\link{approxfun}}.
qgmm.marginal <- function (u, theta, res = 1000, spread = 5, rule = 2) {
d <- theta$d
m <- theta$m
n.samples <- round(res * theta$pie)
n.samples[n.samples == 0] <- 2
# Create grid of evalutation
s <- NULL
for (i in 1:d) {
for (j in 1:m) {
m.ij <- theta$mu[[j]][i]
sd.ij <- sqrt(theta$sigma[[j]][i,i])
s <- c(s, seq(m.ij-spread*sd.ij, m.ij+spread*sd.ij, l = n.samples[j]))
}
}
dim(s) <- c(sum(n.samples), d)
# Evaluate on cdf on the grid
eval <-
pgmm_marginal(z = s, mus = theta$mu, sigmas = theta$sigma, pie = theta$pie)
# Invert function
z.out <- NULL
for (j in 1:d) {
z.out <- c(z.out, approxfun(eval[, j], s[, j], rule = rule)(u[, j]))
}
z.out.is.na <- is.na(z.out)
if (any(z.out.is.na)) {
z.out[z.out.is.na & u >= 1 - 1e-5] <- Inf
z.out[z.out.is.na & u <= 0 + 1e-5] <- -Inf
}
dim(z.out) <- c(nrow(u),d)
return(z.out)
}
# qgmm.marginal2 <- function (u, theta, res = 1000, spread = 5)
# {
#
# d <- theta$d
# m <- theta$m
# n.samples <- round(res * theta$pie)
# n.samples[n.samples == 0] <- 2
#
# # Create grid of evalutation
# s <- NULL
# for (i in 1:d) {
# for (j in 1:m) {
# m.ij <- theta$mu[[j]][i]
# sd.ij <- sqrt(theta$sigma[[j]][i,i])
# s <- c(s, seq(m.ij-spread*sd.ij, m.ij+spread*sd.ij, l = n.samples[j]))
# }
# }
# dim(s) <- c(sum(n.samples), d)
#
# # Evaluate on cdf on the grid
# eval <- pgmm.marginal(z = s, theta = theta)
#
# # Invert function
# z.out <- NULL
# for (j in 1:d) {
# z.out <- c(z.out, approxfun(eval[, j], s[, j], rule = 2)(u[, j]))
# }
# if (any(is.na(z.out))) {
# z.out[is.na(z.out) & u >= 1] <- Inf
# z.out[is.na(z.out) & u <= 0] <- -Inf
# }
# dim(z.out) <- c(nrow(u),d)
# return(z.out)
# }
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.