Nothing
nlm.control <- function(d, typsize = rep(1, d), fscale = 1, print.level = 0,
ndigit = 12, gradtol = 1e-6, stepmax = max(1000 * sqrt(sum((d/typsize)^2)), 1000),
steptol = 1e-6, iterlim = 200, check.analyticals = TRUE){
list(typsize = typsize, fscale = fscale, print.level = print.level,
ndigit = ndigit, gradtol = gradtol, stepmax = stepmax,
steptol = steptol, iterlim = iterlim,
check.analyticals = check.analyticals)
}
nlopt <- function(par, fn, ..., method = c("nlminb", "nlm", "Nelder-Mead", "BFGS"),
control = list(), hesscontrol = list(), lower = -Inf, upper = Inf){
## optimizer
method <- match.arg(method)
if(method == "nlm"){
if (!is.null(control)) {
if (!is.list(control)) {
stop("when specified, 'nlmcontrol' must be a list")
}
pardim <- length(par)
nlmctrl <- do.call("nlm.control", c(d = pardim, control))
}
nlmcall <- list(f = fn, p = par, ..., hessian = TRUE)
nlmcall <- c(nlmcall, nlmctrl)
optres <- do.call(nlm, nlmcall)
par <- optres$estimate
hess <- optres$hessian
minimum <- optres$minimum
if(optres$code < 3){
conv <- 0
} else {
conv <- 1
}
} else {
if(method == "nlminb"){
optres <- nlminb(start = par, objective = fn, ..., control=control,
lower=lower, upper=upper)
par <- optres$par
hess <- optimHess(par, fn, gr=NULL, ..., control = hesscontrol)
hess <- 0.5 * (hess + t(hess))
minimum <- optres$objective
} else {
optres <- optim(par, fn, ..., method = method,
control = control, hessian = TRUE)
par <- optres$par
hess <- optres$hessian
minimum <- optres$value
}
conv <- optres$convergence
}
list(par=par, hessian=hess, minimum=minimum, conv = conv)
}
iterLap.control <- function(pardim, gridSize = ceiling(50*pardim^1.25),
delta = 0.01, maxDim = 20, info = 0, eps = 0.005){
## gridSize - number of grid points for new components
## delta - stopping criterion based on maximum error on the grid
## maxDim - Maximum dimension for iterative process
## info - Plot information during function execution (possible values 0,1,2)
## eps - percentage for stopping criterion based on normalization constant
list(gridSize=gridSize, delta=delta, maxDim=maxDim, eps=eps, info=info)
}
GRApprox <- function(post, start, grad,
method = c("nlminb", "nlm", "Nelder-Mead", "BFGS"),
control = list(), ...){
## post - (non-normalized) log-posterior density
## start - vector of starting values if dimension=1
## otherwise matrix of starting values with
## the starting values in the rows
## grad - gradient (currently not implemented)
## control - control parameters for nlopt
##
## ... - additional arguments for post function
method <- match.arg(method)
if(!is.matrix(start)){
start <- as.matrix(start, ncol=1)
}
nStart <- nrow(start)
d <- ncol(start)
post2 <- function(pars, ...){
-post(pars, ...)
}
obj <- list()
for(j in 1:nStart){
nlcall <- list(par=start[j,], fn = post2, ..., method = method)
nlcall <- c(nlcall, control)
obj[[j]] <- do.call(nlopt, nlcall)
}
means <- matrix(0, nrow = nStart, ncol = d)
sigmainv <- sigmas <- eigenHess <- list()
weights <- dets <- numeric(nStart)
z <- 1
for(j in 1:nStart){
eg <- eigen(obj[[j]]$hessian, symmetric = TRUE)
detSig <- 1/prod(eg$values)
if(all(eg$values > 0)){ # only include maxima of posterior
m <- obj[[j]]$par
sw <- sweep(means, 2, m) # check whether found optimum coincides with already found optimum
mabs <- max(sum(abs(m)), 1)
an <- apply(sw, 1, function(x) sum(abs(x)))
ind1 <- j == 1
ind2 <- all(an[1:z] > 0.01*mabs)
if( ind1 | (ind2 & ind2) ){
means[z,] <- m
sigmas[[z]] <- solve(obj[[j]]$hessian)
sigmainv[[z]] <- obj[[j]]$hessian
eigenHess[[z]] <- eg
dets[z] <- detSig
weights[z] <- detSig^0.5*exp(-obj[[j]]$minimum)
z <- z + 1
}
}
}
means <- means[1:(z-1), , drop=FALSE]
weights <- weights[1:(z-1)]/sum(weights[1:(z-1)])
out <- list(weights=weights, means=means,
sigmas=sigmas, eigenHess=eigenHess,
dets = dets[1:(z-1)], sigmainv = sigmainv)
class(out) <- "mixDist"
out
}
print.mixDist <- function(x, digits = 4, ...){
dims <- dim(x$sigmas[[1]])[1]
ncomp <- length(x$weights)
cat("Mixture Approximation\n\n")
cat("Dimensions:", dims, "\n")
cat("Number of Components:", ncomp, "\n")
cat("Means:\n")
mns <- x$means
nams <- paste("Component ", 1:ncomp, ":", sep="")
dimnames(mns) <- list(nams, 1:dims)
print(round(mns, digits))
}
summary.mixDist <- function(object, digits = 4, ...){
class(object) <- "summary.mixDist"
object
}
print.summary.mixDist <- function(x, digits = 4, ...){
dims <- dim(x$sigmas[[1]])[1]
ncomp <- length(x$weights)
cat("Mixture Approximation\n\n")
cat("Number of Components:", ncomp, "\n")
for(m in 1:ncomp){
cat("------------------------------\n")
cat("Component", m, "\n")
cat("Weight:", signif(x$weights[m], digits), "\n")
cat("Means and Standard Variations:\n")
mns <- x$means[m,]
sdev <- sqrt(diag(x$sigmas[[m]]))
pmat <- rbind(means=mns, sdevs = sdev)
dimnames(pmat) <- list(c("means", "sdev"), 1:dims)
print(round(pmat, digits))
if(dims > 1){
cat("Correlation Matrix:\n")
corMat <- cov2cor(x$sigmas[[m]])
dimnames(corMat) <- list(1:dims, 1:dims)
corMat <- round(corMat, digits)
cf <- format(corMat, digits = digits)
cf[row(cf) > col(cf)] <- ""
print(cf, quote = FALSE)
}
cat("\n")
}
}
getGrid <- function (n, mean, eigenSigma, hess = TRUE,
method = c("Pseudo-Rand", "Quasi-Rand"),
df = Inf, ctrl = list(seed=4711, scrambling=1)){
if (length(mean) != nrow(eigenSigma$vectors)) {
stop("mean and sigma have non-conforming size")
}
if(hess){
ord <- order(eigenSigma$values)
eg <- eigenSigma
eigenSigma$values <- 1/eg$values[ord]
eigenSigma$vectors <- eg$vectors[,ord]
}
dim <- length(mean)
method <- match.arg(method)
if(method == "Pseudo-Rand"){
if(df == Inf){
vals <- rnorm(n * dim)
} else {
vals <- rt(n * dim, df = df)
}
Mat <- matrix(vals, nrow = n)
}
if(method == "Quasi-Rand"){
if(method == "Quasi-Rand")
Mat <- sobol(n, dim = dim,
scrambling = ctrl$scrambling,
seed = ctrl$seed)
if(df == Inf){
Mat <- qnorm(Mat)
} else {
Mat <- qt(Mat, df = df)
}
}
ev <- eigenSigma
retval <- ev$vectors %*% diag(sqrt(ev$values), length(ev$values)) %*% t(ev$vectors)
retval <- Mat %*% retval
retval <- sweep(retval, 2, mean, "+")
retval
}
dmvtAdapted <- function (x, obj, ind = NULL, mat = FALSE, df = 4){
if (is.vector(x)) {
x <- matrix(x, nrow = length(x))
}
dm <- length(obj$sigmas)
if(is.null(ind)){
ind <- 1:dm
} else {
if(!mat)
stop("Cannot calculate density.")
}
out <- matrix(nrow = nrow(x), ncol = length(ind))
z <- 1
for(i in ind){
sigmaInv <- obj$sigmainv[[i]]
detSigma <- obj$dets[i]
mean <- obj$means[i,]
m <- NCOL(sigmaInv)
centr <- sweep(x, 2, mean)
distval <- rowSums((centr %*% sigmaInv) * centr)
logdet <- log(detSigma)
logretval <- lgamma((m + df)/2) - (lgamma(df/2) + 0.5 * (logdet +
m * logb(pi * df))) - 0.5 * (df + m) * logb(1 + distval/df)
out[,z] <- exp(logretval)
z <- z + 1
}
if(mat){
return(out)
} else {
out%*%obj$weights
}
}
dmvtAdaptedC <- function(x, obj, df){
mns <- as.double(t(obj$means))
sigmainv <- as.double(unlist(obj$sigmainv))
dets <- as.double(obj$dets)
weights <- as.double(obj$weights)
ncomp <- as.integer(length(weights))
dm <- as.integer(length(mns)/ncomp)
df <- as.integer(df)
out <- .C("calcmixtdens", as.double(x), dm, ncomp, mns,
sigmainv, dets, weights, df, double(dm), double(1))[[10]]
out
}
dmvnormAdapted <- function (x, obj, ind = NULL, mat = FALSE){
if (is.vector(x)) {
x <- matrix(x, ncol = length(x))
}
dm <- length(obj$sigmas)
if(is.null(ind)){
ind <- 1:dm
} else {
if(!mat)
stop("Cannot calculate density.")
}
out <- matrix(nrow = nrow(x), ncol = length(ind))
z <- 1
for(i in ind){
sigmaInv <- obj$sigmainv[[i]]
detSigma <- obj$dets[i]
mean <- obj$means[i,]
m <- NCOL(sigmaInv)
centr <- sweep(x, 2, mean)
distval <- rowSums((centr %*% sigmaInv) * centr)
logdet <- log(detSigma)
logretval <- -m/2*log(2*pi)-0.5*logdet-0.5*distval
out[,z] <- exp(logretval)
z <- z + 1
}
if(mat){
return(out)
} else {
out%*%obj$weights
}
}
dmvnormAdaptedC <- function(x, obj){
mns <- as.double(t(obj$means))
sigmainv <- as.double(unlist(obj$sigmainv))
dets <- as.double(obj$dets)
weights <- as.double(obj$weights)
ncomp <- as.integer(length(weights))
dm <- as.integer(length(mns)/ncomp)
out <- .C("calcmixmvdens", as.double(x), dm, ncomp, mns,
sigmainv, dets, weights, double(dm), double(1))[[9]]
out
}
LSfitgeq0 <- function(y, X){
D <- crossprod(X)
d <- crossprod(y, X)
Amat <- diag(dim(X)[2])
bvec <- numeric(dim(X)[2])
out <- try(solve.QP(D, d, Amat, bvec)$solution, silent = TRUE)
if(!inherits(out, "try-error")){
return(out)
} else { # if solve.QP does not work use heuristic approach
cf <- qr.coef(qr(X), y)
ind <- cf < 0
if(any(ind)){
cf <- qr.coef(qr(X[,!ind]), y)
out <- numeric(ncol(X))
out[!ind] <- cf
out[out < 0] <- 0
}
return(out)
}
}
samplingMixDist <- function(nSim, obj, df = 4){
## nSim - number of simulations
## obj - mixDist object
if(!inherits(obj, "mixDist"))
stop("obj needs to be of class mixDist")
weights <- obj$weights
means <- obj$means
egHess <- obj$eigenHess
k <- length(weights)
kmat <- rmultinom(1, nSim, weights)
samples <- NULL
for(i in 1:k){
if(kmat[i,1] > 0){
mat0 <- rmvtAdapt(kmat[i,1], means[i,], egHess[[i]], df = df)
samples <- rbind(samples, mat0)
}
}
ind <- sample(1:nSim)
samples[ind,]
}
IS <- function(obj, nSim, df = 4, post, vectorized = FALSE,
cores = 1, ...){
if(!inherits(obj, "mixDist"))
stop("obj needs to be of class mixDist")
samp <- samplingMixDist(nSim, obj, df = df)
if(vectorized){
vals1 <- post(samp,...)
} else {
sampList <- split(samp, 1:nSim)
if(cores == 1){
vals1 <- lapply(sampList, function(x) post(x, ...))
} else {
vals1 <- mclapply(sampList, function(x) post(x, ...), mc.cores = cores)
}
vals1 <- unlist(vals1)
}
scal <- max(vals1)
vals1 <- exp(vals1 - scal)
if(df == Inf){
vals2 <- dmvnormAdapted(samp, obj)
} else {
vals2 <- dmvtAdapted(samp, obj, df = df)
}
w <- vals1/vals2
const <- sum(w)
w <- w/const
out <- list(samp = samp, w = w, normconst = const/nSim*exp(scal),
ESS = 1/sum(w^2))
class(out) <- "IS"
out
}
IMH <- function(obj, nSim, df = 4, post, vectorized = FALSE, cores = 1, ...){
if(!inherits(obj, "mixDist"))
stop("obj needs to be of class mixDist")
samp <- samplingMixDist(nSim, obj, df = df)
if(is.vector(samp)){
samp <- matrix(samp, nrow = length(samp))
}
if(vectorized){
vals1 <- post(samp,...)
} else {
sampList <- split(samp, 1:nSim)#as.list(data.frame(t(samp)))
if(cores == 1){
vals1 <- lapply(sampList, function(x) post(x, ...))
} else {
#vals1 <- mclapply(sampList, function(x) post(x, ...), mc.cores = cores)
stop("currently only cores = 1 possible")
}
vals1 <- unlist(vals1)
}
scal <- max(vals1)
vals1 <- exp(vals1 - scal)
vals2 <- dmvtAdapted(samp, obj, df = df)
w <- vals1/vals2
sampOut <- sampOut <- matrix(0, nrow = nSim, ncol = ncol(samp))
wcur <- w[1]
sampCur <- sampOut[1,] <- samp[1,]
u <- runif(nSim)
z <- 1
for(i in 2:nSim){
if(u[i] < w[i]/wcur){
sampOut[i,] <- sampCur <- samp[i,]
wcur <- w[i]
z <- z + 1
} else {
sampOut[i,] <- sampCur
}
}
const <- sum(w)
w <- w/const
out <- list(samp = sampOut, w = w, normconst = const/nSim*exp(scal),
accept = z/nSim)
class(out) <- "IMH"
out
}
rmvtAdapt <- function(n, mean = NULL, eigenSigma = NULL, df){
## random numbers from multivariate t-distribution
dm <- length(eigenSigma$values)
grd <- getGrid(n, mean = rep(0, dm), eigenSigma, hess = TRUE,
method = "Pseudo-Rand")
if(df != Inf){
sqchsq <- sqrt(rchisq(n, df)/df)
} else {
sqchsq <- rep(1,n)
}
grd <- grd/sqchsq
mnmat <- matrix(mean, byrow = T, nrow = n, ncol = dm)
grd + mnmat
}
resample <- function(n, obj){
if(!inherits(obj, "IS"))
stop("obj needs to be of class IS")
samples <- obj$samp
wgts <- obj$w
## performs residual resampling
if(is.vector(samples)){
samples <- matrix(samples, nrow = length(samples))
}
nS <- nrow(samples)
reps <- floor(wgts*n)
if(!all(reps==0)){
inds <- rep(1:nS, reps)
probAdjusted <- (wgts*n-reps)/(n-sum(reps))
inds2 <- sample(1:nS, n-sum(reps), prob=probAdjusted, replace=TRUE)
inds <- c(inds, inds2)
} else {
inds <- sample(1:nS, n, prob=wgts, replace=TRUE)
}
samples[inds,]
}
iterLap <- function(post, ..., GRobj = NULL,
vectorized = FALSE, startVals = NULL,
method = c("nlminb", "nlm", "Nelder-Mead", "BFGS"),
control = NULL, nlcontrol = list()){
## GRobj - object of class mixDist
## post - log-posterior
## control - control parameters for iterated Laplace approximation
## nlcontrol - control parameters for nlm
method <- match.arg(method)
## extract dimension and convert startVals to matrix
if(!is.null(startVals)){
if(!is.matrix(startVals)){
startVals <- as.matrix(startVals, ncol=1)
}
pardim <- ncol(startVals)
}
## calculcate (or extract) GR object
if(is.null(GRobj)){
if(is.null(startVals)){
stop("need to specify starting matrix if GRobj is missing")
}
GRobj <- GRApprox(post, startVals, method = method, control=nlcontrol, ...)
}
obj <- GRobj # obj contains current approximation
pardim <- ncol(obj$means)
## check control list
if (!is.null(control)) {
if(!is.list(control)) {
stop("when specified, 'control' must be a list")
}
ctrl <- do.call("iterLap.control", c(pardim=pardim,control))
} else {
ctrl <- iterLap.control(pardim = pardim)
}
propnorm <- numeric(ctrl$maxDim)
## calculate starting grid (from GR approximation)
dm <- length(obj$weights)
grid <- NULL
for(i in 1:dm){
grd <- getGrid(ctrl$gridSize, obj$means[i,], obj$eigenHess[[i]], hess = TRUE,
method = "Quasi-Rand", df=Inf)
grid <- rbind(grid, grd)
}
## evaluate kernel and mixture of t-distributions at starting grid
if(vectorized){
vals <- post(grid, ...)
} else {
vals <- apply(grid, 1, function(x) post(x, ...))
}
lgMax <- max(vals)
X <- dmvnormAdapted(grid, obj, mat = TRUE)
z <- dm
## determine weights
obj$weights <- LSfitgeq0(exp(vals-lgMax), X)
propnorm[z] <- sum(obj$weights)*exp(lgMax)
## iterate number of components
while(z < ctrl$maxDim){
if(ctrl$info > 0){
cat("Add Component: ", z+1, "\n")
}
pred <- X%*%obj$weights # calculate predictions
if(checkStop(vals, lgMax, pred, ctrl, propnorm, z)){ # check stop criterion
break
} else { # calculate starting values for optimizer of residual
startOpt <- getStartVal(vals, lgMax, pred, grid, obj$means[z,], z)
}
## run optimizer
optres <- opt(obj, z, startOpt, method, ctrl, nlcontrol, post, lgMax, ...)
if(!optres$add){ # optimizer output does not contain local optimum
if(ctrl$info > 0){
cat("Cannot improve approximation.\n")
}
break
} else {
## add component to approximation
obj <- optres$obj
z <- z + 1
## evaluate new component on current grid (ind = z)
Xgrd <- dmvnormAdapted(grid, obj, ind = z, mat = TRUE)
X <- cbind(X, Xgrd)
## calculate grid of new component
grdNew <- getGrid(ctrl$gridSize, obj$means[z,], obj$eigenHess[[z]],
hess = TRUE, method = "Quasi-Rand", df = Inf)
## evaluate all components on new grid
Xnew2 <- dmvnormAdapted(grdNew, obj, mat = TRUE)
X <- rbind(X, Xnew2)
grid <- rbind(grid, grdNew)
## evaluate target function at new component
if(vectorized){
valsnew <- post(grdNew, ...)
} else {
valsnew <- apply(grdNew, 1, function(x) post(x, ...))
}
if(any(is.na(valsnew)) | any(abs(valsnew)==Inf)){
warning("Inf/NA produced.")
}
vals <- c(vals, valsnew)
lgMax <- max(vals)
obj$weights <- LSfitgeq0(exp(vals-lgMax), X)
propnorm[z] <- sum(obj$weights)*exp(lgMax)
if(ctrl$info > 0){
cat("Norm. Constant of Approx.: ", signif(propnorm[z], 5), "\n\n")
}
}
}
obj$weights <- abs(obj$weights/sum(abs(obj$weights)))
class(obj) <- "mixDist"
obj
}
opt <- function(obj, z, start, method, control, nlcontrol, post, lgMax, ...){
## optimizer
post2 <- function(pars, obj, eps = 0.0001, ...){
val <- exp(post(pars, ...)-lgMax)-dmvnormAdaptedC(pars, obj)
if(val < eps | is.na(val)){
val <- exp(val-eps)*eps
}
return(-log(val))
}
j <- 0
while(j < nrow(start)){
j <- j + 1
nlcall <- list(par=start[j,], fn = post2, ...,
method = method, obj = obj, control = nlcontrol)
optobj <- do.call(nlopt, nlcall)
eg <- eigen(optobj$hessian, symmetric = TRUE)
detSig <- 1/prod(eg$values)
if(all(eg$values > 0)){ # only include maxima of posterior
m <- optobj$par
sw <- sweep(obj$means, 2, m) # check whether found optimum coincides with already found optimum
mabs <- sum(abs(m))
an <- apply(sw, 1, function(x) sum(abs(x)))
ind2 <- all(an[1:z] > 0.005*mabs)
if(ind2){
if(control$info > 1){
cat("Added Component:", signif(m,5), "\n")
}
obj$means <- rbind(obj$means, m)
obj$sigmas[[z+1]] <- solve(optobj$hessian)
obj$sigmainv[[z+1]] <- optobj$hessian
obj$eigenHess[[z+1]] <- eg
obj$dets <- c(obj$dets, detSig)
add <- TRUE
break
} else {
if(control$info > 1){
cat("Too close to already used components\n")
}
add <- FALSE
}
} else {
if(control$info > 1){
cat("Negative Eigenvalues\n")
}
add <- FALSE
}
}
list(obj = obj, add = add)
}
calcDist <- function(start, oldVal){
dfs <- sweep(start, 2, oldVal)
rev(order(rowSums(dfs^2)))
}
extractMoments <- function(obj){
mn <- obj$weights%*%obj$means
sigm <- mn2 <- matrix(0, nrow = ncol(obj$means), ncol = ncol(obj$means))
for(j in 1:length(obj$weights)){
sigm <- sigm + obj$weights[j]*obj$sigmas[[j]]
mn2 <- mn2 + obj$weights[j]*obj$means[j,]%*%t(obj$means[j,])
}
sigm <- sigm + mn2 - t(mn)%*%mn
list(mean = as.numeric(mn), sigma = sigm)
}
checkStop <- function(vals, scal, pred, control, propnorm, z){
out <- FALSE
if(max(exp(vals-scal) - pred) < control$delta){
out <- TRUE
if(control$info > 0){
cat("Maximum Error on grid smaller than delta\n")
}
}
if(z >= 3){
dif <- abs(propnorm[z]-mean(propnorm[(z-1):(z-2)]))
std <- propnorm[z]
if(dif/std < control$eps){
if(control$info > 0){
cat("Approximation of normalization constant does not improve\n")
}
out <- TRUE
}
}
out
}
getStartVal <- function(vals, scal, pred, grid, oldVal, init){
## determine 10 grid point with largest importance weight
ind <- rev(order(exp(vals-scal) / pred))[1:10]
## select starting values for k-means
vals <- c(1,9,3,7,2,8,4,6,5,6,1,5,8,4,7,3,9,2,4,6,5,2,8,3,7,1,9)
val <- vals[(init-1)%%27+1] # emulate some randomness for ind0
ind0 <- c(val, val+3, val+6)%%10+1 # this increases diversity
startkm <- grid[ind[ind0],]
## reduce this to 3 possibly distinct values
start <- kmeans(grid[ind,], startkm)$centers
## sort those 3 values according to their distance to last point
ind2 <- calcDist(start, oldVal)
start[ind2, , drop=FALSE]
}
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.