-
Notifications
You must be signed in to change notification settings - Fork 5
/
efficient.closures.R
88 lines (74 loc) · 2.2 KB
/
efficient.closures.R
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
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
# 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
}
}
.make.quadratic.penalty <- function(priorParameters) {
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
}
}