# R/efficient.closures.R In texmex: Statistical Modelling of Extreme Values

```# Herewith a bunch of factory functions that are used for making priors

.make.mvn.prior <- function(priorParameters) {
mean <- priorParameters[[1]]
cov <- priorParameters[[2]]
if (!isSymmetric(cov, tol=sqrt(.Machine\$double.eps),
check.attributes = FALSE)) {
stop("Covariance must be a symmetric matrix")
}
# we should also check for a positive definite matrix
# but hey
factors <- svd(as.matrix(cov))
if (min(factors\$d) <= 0) {
stop("Covariance must be nonsingular.")
}
num.vars <- length(mean)
if (ncol(cov) != num.vars) {
stop("Covariance and mean must be the same size.")
}

logdet <- sum(log(factors\$d))
base.ret <- -(num.vars * log(2 * pi) + logdet) / 2
stacked <- rbind(t(factors\$u), t(factors\$v))

function(param) {
delta <- param - mean
prod <- as.numeric(stacked %*% delta)
first <- prod[1:num.vars]
last <- prod[(1 + num.vars):(2*num.vars)]
mahab <- as.numeric(first %*% (last / factors\$d))
base.ret - 0.5 * mahab
}
}

centre <- priorParameters[[1]]
cov <- priorParameters[[2]]

factors <- svd(as.matrix(cov))
if (min(factors\$d) == 0) {
stop("Singular covariance matrix: quadratic penalty impossible")
}

stacked <- rbind(t(factors\$u), t(factors\$v))
n.prior <- length(centre)
function(param) {
delta <- param - centre
prod <- as.numeric(stacked %*% delta)
first <- prod[1:n.prior]
last <- prod[(1 + n.prior):(2*n.prior)]
as.numeric(first %*% (last / factors\$d))
}
}

.make.lasso.penalty <- function(priorParameters) {
centre <- priorParameters[[1]]
cov <- diag(priorParameters[[2]])
function(param) {
sum(abs(param - centre) * cov)
}
}

.make.dummy.penalty <- function(priorParameters) {
function(param) {
0
}
}

.random.spd.matrix <- function(dimn) {
# generate a random symmetric positive definite matrix
rank <- 2 * dimn
while (TRUE) {
mat <- matrix(rexp(dimn * rank), nrow=rank)
cov <- t(mat) %*% mat
# cov is almost surely (technical sense) SPD
# but let's check
res <- eigen(cov, TRUE, only.values=TRUE)
if (min(res\$values) > 0) {
return(cov)
}
# otherwise go round again
}
}
```

## Try the texmex package in your browser

Any scripts or data that you put into this service are public.

texmex documentation built on May 2, 2019, 5:41 a.m.