-
Notifications
You must be signed in to change notification settings - Fork 1
/
metropolisGP.R
165 lines (145 loc) · 5.89 KB
/
metropolisGP.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
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
##' Function to perform Metropolis-Hastings for GP hyperparameters with different priors
##'
##' @title Perform metropolis update for GP hyperparameters
##' @param inith initial hyperparamters
##' @param X The data
##' @param tau The indexing parameters
##' @param nk Number of observations
##' @param D number of samples
##' @param niter Number of MH iteractions
##' @param hyperMean A vector indicating the log-normal means. Default is `c(0,0,0)`.
##' @param hyperSd A vector indicating the log-normal standard deviations. Default is `c(1,1,1)`
##' @return Returns new hyperparamters and the acceptance rate
##'
##' @md
##' @rdname bandle-gp
metropolisGP <- function(inith,
X,
tau,
nk,
D,
niter,
hyperMean = c(0,0,0),
hyperSd = c(1,1,1)
){
h <- matrix(0, ncol = niter + 1, nrow = length(inith))
h[, 1] <- inith
ar <- 0
for(i in seq.int(niter)){
xi <- rnorm(length(inith), mean = 0, sd = 0.1) # sample random walk steps
oldHypers <- h[, i]
proposedHypers <- h[, i] + xi #random walk proposal
#compute metropolis ratio, likelihoodGPcpp return negative logliklihood
mhratio <- -likelihoodGPcpp(X, tau, proposedHypers, nk, D) + likelihoodGPcpp(X, tau, oldHypers, nk, D) +
sum(dnorm(proposedHypers, mean = hyperMean, sd = hyperSd, log = TRUE)) - sum(dnorm(oldHypers, mean = hyperMean,
sd = hyperSd, log = TRUE))
if(mhratio > log(runif(1, 0, 1))){
h[, i + 1] <- proposedHypers
ar <- ar + 1
}else{
h[, i + 1] <- oldHypers
}
}
ar <- ar/niter
return(list(h = h, ar = ar))
}
##' @title Perform metropolis update for GP hyperparameters with matern covariance
##' @param inith initial hyperparamters
##' @param X The data
##' @param tau The indexing parameters
##' @param nk Number of observations
##' @param D number of samples
##' @param niter Number of MH iteractions
##' @param nu Smoothness of the matern covariance
##' @param hyppar A vector indicating the penalised complexity prior hyperparameters.
##' Default is `c(1,1,1)`
##' @param propSd The proposal standard deviation. Default is `c(0.3,0.1,0.1)`. Do not
##' change unless you know what you are doing.
##' @md
##' @rdname bandle-gp
metropolisGPmatern <- function(inith,
X,
tau,
nk,
D,
niter,
nu = 2,
hyppar = c(1, 1, 1),
propSd = c(0.3,0.1,0.1)
){
h <- matrix(0, ncol = niter + 1, nrow = length(inith))
h[, 1] <- inith
ar <- 0
propSd <- propSd
for(i in seq.int(niter)){
# repeat proposal so proposals are positive, (truncated sampler)
repeat {
xi <- rnorm(length(inith), mean = 0, sd = propSd) # sample random walk steps
oldHypers <- h[, i]
proposedHypers <- h[, i] + xi #random walk proposal
if ((all((proposedHypers > 0)))){
break
}
}
#compute metropolis ratio, likelihoodGPcpp return negative logliklihood, densities from PC priors + truncated sampler correction
mhratiolike <- -likelihoodGPmatern(X, tau, proposedHypers, nk, D, materncov = TRUE, nu = nu) +
likelihoodGPmatern(X, tau, oldHypers, nk, D, materncov = TRUE, nu = nu)
mhratioprior <- PCrhomvar(rho = proposedHypers[1], a = proposedHypers[2], lambda1 = hyppar[1],
lambda2 = hyppar[2], log = TRUE) +
Gumbel(1/proposedHypers[3], lambda = hyppar[3], log = TRUE) + sum(dnorm(oldHypers, mean = 0,
sd = propSd, log = TRUE)) -
(PCrhomvar(rho = oldHypers[1], a = oldHypers[2], lambda1 = hyppar[1],
lambda2 = hyppar[2], log = TRUE) + Gumbel(1/oldHypers[3], lambda = hyppar[3], log = TRUE) +
sum(dnorm(proposedHypers, mean = 0, sd = propSd, log = TRUE)))
mhratio <- mhratiolike + mhratioprior
if(mhratio > log(runif(1, 0, 1))){
h[, i + 1] <- proposedHypers
ar <- ar + 1
}else{
h[, i + 1] <- oldHypers
}
}
ar <- ar/niter
# returns parameter values on true scale note that variance rather than precision is returned
return(list(h = h, ar = ar))
}
##' @title Type-2 Gumbel distribution
##' @param x observation
##' @param lambda scale parameter of the type-2 Gumbel distribution
##' @param log `logical` indicating whether to return `log`. Default is `TRUE`
##' @return Returns the likelihood of the type-2 GUmbel distribution
##' @md
##'
##' @rdname bandle-gp
Gumbel <- function(x,
lambda,
log = TRUE){
if (log == FALSE) {
gumbel <- (lambda/2) * x^{-3/2} * exp(-lambda * x^{-1/2})
} else {
gumbel <- log(lambda/2) - 3*log(x)/2 - lambda * x^{-1/2}
}
return(gumbel)
}
##' @title Bivariate penalized complexity prior for length-scale and amplitude
##' @param rho length-scale parameter
##' @param a amplitude
##' @param lambda1 first parameter of distribution
##' @param lambda2 second parameter of distribution
##' @param log `logical` indicating whether to return `log`. Default is `TRUE`
##' @return Returns the likelihood of the bivariate penalised complexity prior
##' @md
##'
##' @rdname bandle-gp
PCrhomvar <- function(rho,
a,
lambda1,
lambda2,
log = TRUE){
if (log == FALSE) {
res <- (lambda1 * lambda2/2) * rho^{- 3/2} * exp(-lambda1 * rho^{-1/2} - lambda2 * a)
} else {
res <- log(lambda1 * lambda2/2) - 3 * log(rho)/2 - lambda1 * rho^{-1/2} - lambda2 * a
}
return(res)
}