A Base R, simple implementation of the No-Underrun Sampler. This package aims to mostly directly implement the algorithm as described by th
Implemented memory trick
+38
-25
+27
-20
R/NURS.R
+27
-20
R/NURS.R
···
39
39
#' @param epsilon density threshold
40
40
#' @param h lattice size
41
41
#' @param M maximum number of doublings
42
-
#' @returns One NURS step from theta.
42
+
#' @returns next draw
43
43
NURS_step <- function(logpdf, theta, epsilon, h, M) {
44
44
d <- length(theta)
45
45
# hit
···
57
57
theta + s * rho else theta
58
58
59
59
log_eps_h <- log(epsilon) + log(h)
60
-
orbit_points <- list(theta0)
61
-
log_vals <- logpdf(theta0)
60
+
logW <- logpdf(theta0)
61
+
theta_tilde <- theta0
62
62
63
63
# bookkeeping to get orbit ends
64
64
left <- right <- 1
65
+
left_point <- right_point <- theta0
66
+
left_val <- right_val <- logW
65
67
66
68
# Orbit selection procedure
67
69
B <- sample(c(FALSE, TRUE), M, replace = TRUE)
···
70
72
n_ext <- 2^(k - 1)
71
73
orbit_ext <- lapply(
72
74
seq_len(n_ext),
73
-
\(i)
74
-
(if (B_k) orbit_points[[right]] else orbit_points[[left]]) +
75
-
(2 * B_k - 1) * i * h * rho
75
+
\(i) (if (B_k) right_point else left_point) + i * (2 * B_k - 1) * h * rho
76
76
)
77
77
log_ext <- sapply(orbit_ext, logpdf)
78
78
79
+
# Recursive sub-stopping criterion
79
80
if (NURS_sub_stop(log_ext, log_eps_h)) break
80
81
82
+
for (i in seq_len(n_ext)) {
83
+
theta_i <- orbit_ext[[i]]
84
+
log_i <- log_ext[i]
85
+
logW_new <- log_sum_exp(c(logW, log_i))
86
+
if (runif(1) < exp(log_i - logW_new)) theta_tilde <- theta_i
87
+
logW <- logW_new
88
+
}
89
+
81
90
if (B_k) {
82
-
orbit_points <- c(orbit_points, orbit_ext)
83
-
log_vals <- c(log_vals, log_ext)
91
+
right_point <- orbit_ext[[n_ext]]
92
+
right_val <- log_ext[n_ext]
84
93
} else {
85
-
orbit_points <- c(orbit_ext, orbit_points)
86
-
log_vals <- c(log_ext, log_vals)
87
-
left <- 1
94
+
left_point <- orbit_ext[[n_ext]]
95
+
left_val <- log_ext[n_ext]
88
96
}
89
-
right <- right + n_ext
90
97
91
-
if (NURS_stop(log_vals, log_eps_h)) break
98
+
# No-underrun stopping criterion
99
+
if (max(left_val, right_val) <= log_eps_h + logW) break
92
100
}
93
101
94
-
orbit_points[[sample(
95
-
length(log_vals),
96
-
1,
97
-
prob = exp(log_vals - log_sum_exp(log_vals))
98
-
)]]
102
+
theta_tilde
99
103
}
100
104
101
105
#' NURS draws
···
103
107
#' @param logpdf log (non-normalized) target density
104
108
#' @param theta_init initial state
105
109
#' @param n number of draws
106
-
#' @param epsilon density threshold
107
-
#' @param h lattice size
110
+
#' @param epsilon non-negative density threshold
111
+
#' @param h positive lattice size
108
112
#' @param M maximum number of doublings
113
+
#' @returns a sequence of draws
114
+
#'
109
115
#' @export
110
116
#' @source <https://arxiv.org/abs/2501.18548v2>
111
117
NURS <- function(logpdf, theta_init, n, epsilon, h, M) {
118
+
stopifnot(epsilon >= 0, h > 0)
112
119
d <- length(theta_init)
113
120
draws <- matrix(NA, n, d)
114
121
draws[1, ] <- theta_init
+2
-1
README.Rmd
+2
-1
README.Rmd
···
18
18
<!-- badges: start -->
19
19
<!-- badges: end -->
20
20
21
-
A Base R, simple implementation of the No-Underrun Sampler. This implementation aims to mostly directly implement the algorithm as described by the paper, with at most small changes for code aesthetics and performance.
21
+
A Base R, simple implementation of the No-Underrun Sampler. This implementation aims to mostly directly implement the algorithm as described by the paper, with at most small changes for code aesthetics and performance. This version uses the memory saving technique described in section 2.2.
22
22
23
23
## Installation
24
24
···
40
40
y <- theta[1]
41
41
dnorm(y, 0, 3, log = TRUE) + sum(dnorm(theta[-1], 0, exp(y / 2), log = TRUE))
42
42
}
43
+
43
44
samples <- NURS(
44
45
logpdf_funnel,
45
46
theta_init = rep(0, 15),
+3
-1
README.md
+3
-1
README.md
···
10
10
A Base R, simple implementation of the No-Underrun Sampler. This
11
11
implementation aims to mostly directly implement the algorithm as
12
12
described by the paper, with at most small changes for code aesthetics
13
-
and performance.
13
+
and performance. This version uses the memory saving technique described
14
+
in section 2.2.
14
15
15
16
## Installation
16
17
···
33
34
y <- theta[1]
34
35
dnorm(y, 0, 3, log = TRUE) + sum(dnorm(theta[-1], 0, exp(y / 2), log = TRUE))
35
36
}
37
+
36
38
samples <- NURS(
37
39
logpdf_funnel,
38
40
theta_init = rep(0, 15),
+5
-2
man/NURS.Rd
+5
-2
man/NURS.Rd
···
16
16
17
17
\item{n}{number of draws}
18
18
19
-
\item{epsilon}{density threshold}
19
+
\item{epsilon}{non-negative density threshold}
20
20
21
-
\item{h}{lattice size}
21
+
\item{h}{positive lattice size}
22
22
23
23
\item{M}{maximum number of doublings}
24
+
}
25
+
\value{
26
+
a sequence of draws
24
27
}
25
28
\description{
26
29
NURS draws
+1
-1
man/NURS_step.Rd
+1
-1
man/NURS_step.Rd
man/figures/README-funnel-1.png
man/figures/README-funnel-1.png
This is a binary file and will not be displayed.