-
Notifications
You must be signed in to change notification settings - Fork 0
/
compute_posterior_probs.Rmd
197 lines (176 loc) · 7.6 KB
/
compute_posterior_probs.Rmd
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
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
---
title: "Computation of posterior probability errors for different approximations of the Matérn covariance function"
author: "David Bolin, Vaibhav Mehandiratta, and Alexandre B. Simas"
date: "Created: 2024-05-20. Last modified: `r Sys.Date()`."
output: rmarkdown::html_vignette
vignette: >
%\VignetteIndexEntry{Computation of posterior probability errors for different approximations of the Matérn covariance function}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteEncoding{UTF-8}
references:
- id: BMS2024
title: "Linear cost and exponentially convergent approximation of Gaussian Matérn stochastic processes"
author:
- family: Bolin
given: David
- family: Mehandiratta
given: Vaibhav
- family: Simas
given: Alexandre B.
container-title: Preprint
- id: DBFG
title: "Hierarchical nearest-neighbor Gaussian process models for large geostatistical datasets"
author:
- family: Datta
given: A.
- family: Benerjee
given: S.
- family: Finley
given: A. O.
- family: Gelfand
given: A. E.
container-title: J. Amer. Statist. Assoc.
type: article
issue: 111
pages: 800-812
issued:
year: 2016
- id: KS
title: "Approximate state-space Gaussian processes via spectral transformation"
author:
- family: Karvonen
given: T.
- family: Sarkkä
given: S.
container-title: 2016 IEEE 26th International Workshop on Machine Learning for Signal Processing (MLSP)
type: article
pages: 1-6
issued:
year: 2016
- id: RR
title: "Random features for large-scale kernel machines"
author:
- family: Rahimi
given: A.
- family: Recht
given: B.
container-title: Advances in neural information processing systems
type: article
issue: 20
issued:
year: 2007
---
```{r setup, include = FALSE}
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
set.seed(1)
```
## Introduction
This vignette contains posterior probability errors.
Our goal is to generate the data as
$$Y_i = u(t_i) + \varepsilon_i,\quad i=1,\ldots,N,$$
where $u(\cdot)$ is a stationary Gaussian process with a Mat\'ern covariance function, $t_1,\ldots,t_N$ are locations on a bounded interval on the real line, and $\varepsilon_1,\ldots,\varepsilon_N$ are measurement errors, such that $\varepsilon_1,\ldots,\varepsilon_N$ are i.i.d., independent of the field $u(\cdot)$ and follow standard normal distribution. Then, we will compute the posterior mean using the exact Matérn covariance function and we will compare to the posterior means obtained by using the different approximations of the covariance function.
We will then compute approximations of the probability
$$P(u(t_i) \geq y_i: i\in I | Y),$$
where $t_i$ are some fixed locations and $I$ is an index set.
## Computing the errors
Let us load all the auxiliary functions:
```{r, message=FALSE}
source("aux_functions/aux_functions_cov.R")
source("aux_functions/probability_computations.R")
```
We will now consider the basic setup:
```{r}
N <- 100
N_prob <- 30
obs_ind <- sort(sample(1:N, N-N_prob))
prob_ind <- which(!(1:N %in% obs_ind))
L <- 10
range <- 1
sigma <- 1
nu_vec=c(0.6,0.8,1,1.2)#seq(0.6, 1.2, by = 0.1)
idx <- (nu_vec + 0.5)%%1 > 1e-10
nu_vec <- nu_vec[idx]
samples = 10
m <- 1:6
# Std. dev. of the measurement error
sigma_e <- 0.01
loc <- seq(0,L, length.out = N)
m_nngp_fun <- function(m, alpha){
if(alpha<1) {
mn <- m
} else if (alpha < 2) {
mn <- round(m*(ceil(alpha)+1)^1.5)
} else {
mn <- round(m*(ceil(alpha)+1)^1.5)
}
return(mn)
}
m_pca_fun <- function(m, alpha){
return(pmax(9*round((m+1)*sqrt(m) * ceil(alpha)^(3/2)),2))
}
m_statespace_fun <- function(m, alpha){
return(m+floor(alpha) + 1)
}
m_fourier_fun <- function(m, alpha){
return(m+floor(alpha) + 1)
}
```
We will now compute the probability estimates based on the above calibration, which is taken with respect to prediction timings.
It is important to observe that we are computing the posterior probabilities for nearest neighbor Gaussian models on locations that are not supporting points. That is, we assume that the location in which we want to do prediction is a new location to the model.
Note that it is possible to have unobserved locations that are supporting points for the nearest neighbor Gaussian processes, however in this case, a practitioner would need to know he/she would be interested in studying that location at the time they were building the model, which is an unrealistic assumption.
```{r}
prob_ests <- compute_prob_ests(type = "unobserved", loc_full = loc, obs_ind = obs_ind, loc_prob_idx = prob_ind, L_statespace = L, range = range, sigma = sigma, nu_vec = nu_vec, sigma_e = sigma_e, m_rat = m, m_nngp_fun = m_nngp_fun,
m_pca_fun = m_pca_fun, m_fourier_fun = m_fourier_fun, m_statespace_fun = m_statespace_fun, samples = 10)
```
```{r}
err.rat <- unlist(prob_ests[[1]]$rat) - prob_ests[[1]]$true
err.nn <- unlist(prob_ests[[1]]$nngp) - prob_ests[[1]]$true
plot(m, err.rat, ylim = c(min(c(min(err.rat), min(err.nn))),
max(c(max(err.rat), max(err.nn)))), main = nu_vec[1])
lines(m, err.nn, col = 2)
err.rat <- unlist(prob_ests[[2]]$rat) - prob_ests[[2]]$true
err.nn <- unlist(prob_ests[[2]]$nngp) - prob_ests[[2]]$true
plot(m, err.rat, ylim = c(min(c(min(err.rat), min(err.nn))),
max(c(max(err.rat), max(err.nn)))), main = nu_vec[2])
lines(m, err.nn, col = 2)
err.rat <- unlist(prob_ests[[3]]$rat) - prob_ests[[3]]$true
err.nn <- unlist(prob_ests[[3]]$nngp) - prob_ests[[3]]$true
plot(m, err.rat, ylim = c(min(c(min(err.rat), min(err.nn))),
max(c(max(err.rat), max(err.nn)))), main = nu_vec[3])
lines(m, err.nn, col = 2)
err.rat <- unlist(prob_ests[[4]]$rat) - prob_ests[[4]]$true
err.nn <- unlist(prob_ests[[4]]$nngp) - prob_ests[[4]]$true
plot(m, err.rat, ylim = c(min(c(min(err.rat), min(err.nn))),
max(c(max(err.rat), max(err.nn)))), main = nu_vec[4])
lines(m, err.nn, col = 2)
```
Let us now compute the posterior probabilities with respect to all locations, that is, both observed and unobserved locations:
```{r}
prob_ests_all <- compute_prob_ests(type = "all", loc_full = loc, obs_ind = obs_ind, loc_prob_idx = prob_ind, L_statespace = L, range = range, sigma = sigma, nu_vec = nu_vec, sigma_e = sigma_e, m_rat = m, m_nngp_fun = m_nngp_fun,
m_pca_fun = m_pca_fun, m_fourier_fun = m_fourier_fun, m_statespace_fun = m_statespace_fun, samples = 10)
```
```{r}
err.rat <- unlist(prob_ests_all[[1]]$rat) - prob_ests_all[[1]]$true
err.nn <- unlist(prob_ests_all[[1]]$nngp) - prob_ests_all[[1]]$true
plot(m, err.rat, ylim = c(min(c(min(err.rat), min(err.nn))),
max(c(max(err.rat), max(err.nn)))), main = nu_vec[1])
lines(m, err.nn, col = 2)
err.rat <- unlist(prob_ests_all[[2]]$rat) - prob_ests_all[[2]]$true
err.nn <- unlist(prob_ests_all[[2]]$nngp) - prob_ests_all[[2]]$true
plot(m, err.rat, ylim = c(min(c(min(err.rat), min(err.nn))),
max(c(max(err.rat), max(err.nn)))), main = nu_vec[2])
lines(m, err.nn, col = 2)
err.rat <- unlist(prob_ests_all[[3]]$rat) - prob_ests_all[[3]]$true
err.nn <- unlist(prob_ests_all[[3]]$nngp) - prob_ests_all[[3]]$true
plot(m, err.rat, ylim = c(min(c(min(err.rat), min(err.nn))),
max(c(max(err.rat), max(err.nn)))), main = nu_vec[3])
lines(m, err.nn, col = 2)
err.rat <- unlist(prob_ests_all[[4]]$rat) - prob_ests_all[[4]]$true
err.nn <- unlist(prob_ests_all[[4]]$nngp) - prob_ests_all[[4]]$true
plot(m, err.rat, ylim = c(min(c(min(err.rat), min(err.nn))),
max(c(max(err.rat), max(err.nn)))), main = nu_vec[4])
lines(m, err.nn, col = 2)
```