A Base R, simple implementation of the No-Underrun Sampler. This package aims to mostly directly implement the algorithm as described by th
first commit
+190
+1
.gitignore
+1
.gitignore
···
1
+
.Rproj.user
+10
DESCRIPTION
+10
DESCRIPTION
···
1
+
Package: R.NURS
2
+
Title: What the Package Does (One Line, Title Case)
3
+
Version: 0.0.0.9000
4
+
Authors@R:
5
+
person("First", "Last", , "first.last@example.com", role = c("aut", "cre"))
6
+
Description: What the package does (one paragraph).
7
+
License: MIT + file LICENSE
8
+
Encoding: UTF-8
9
+
Roxygen: list(markdown = TRUE)
10
+
RoxygenNote: 7.3.2
+21
LICENSE.md
+21
LICENSE.md
···
1
+
# MIT License
2
+
3
+
Copyright (c) 2025 R.NURS authors
4
+
5
+
Permission is hereby granted, free of charge, to any person obtaining a copy
6
+
of this software and associated documentation files (the "Software"), to deal
7
+
in the Software without restriction, including without limitation the rights
8
+
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
9
+
copies of the Software, and to permit persons to whom the Software is
10
+
furnished to do so, subject to the following conditions:
11
+
12
+
The above copyright notice and this permission notice shall be included in all
13
+
copies or substantial portions of the Software.
14
+
15
+
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
16
+
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
17
+
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
18
+
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
19
+
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
20
+
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
21
+
SOFTWARE.
+17
R.NURS.Rproj
+17
R.NURS.Rproj
···
1
+
Version: 1.0
2
+
3
+
RestoreWorkspace: No
4
+
SaveWorkspace: No
5
+
AlwaysSaveHistory: Default
6
+
7
+
EnableCodeIndexing: Yes
8
+
Encoding: UTF-8
9
+
10
+
AutoAppendNewline: Yes
11
+
StripTrailingWhitespace: Yes
12
+
LineEndingConversion: Posix
13
+
14
+
BuildType: Package
15
+
PackageUseDevtools: Yes
16
+
PackageInstallArgs: --no-multiarch --with-keep.source
17
+
PackageRoxygenize: rd,collate,namespace
+92
R/NURS.R
+92
R/NURS.R
···
1
+
log_sum_exp <- function(log_vals) {
2
+
if (requireNamespace("matrixStats", quietly = TRUE)) {
3
+
matrixStats::logSumExp(log_vals)
4
+
} else {
5
+
m <- max(log_vals)
6
+
m + log(sum(exp(log_vals - m)))
7
+
}
8
+
}
9
+
10
+
# No underrun condition
11
+
NURS_stop <- function(log_vals, log_eps_h) {
12
+
max(log_vals[1], log_vals[length(log_vals)]) <=
13
+
log_eps_h + log_sum_exp(log_vals)
14
+
}
15
+
16
+
# Recursive sub stopping
17
+
NURS_sub_stop <- function(log_vals, log_eps_h) {
18
+
n <- length(log_vals)
19
+
if (n < 2) return(FALSE)
20
+
if (NURS_stop(log_vals, log_eps_h)) return(TRUE)
21
+
left_indices <- 1:(n / 2)
22
+
NURS_sub_stop(log_vals[left_indices], log_eps_h) ||
23
+
NURS_sub_stop(log_vals[-left_indices], log_eps_h)
24
+
}
25
+
26
+
NURS_step <- function(logpdf, theta, epsilon, h, M) {
27
+
d <- length(theta)
28
+
# hit
29
+
z <- rnorm(d)
30
+
rho <- z / sqrt(sum(z^2))
31
+
32
+
# run
33
+
s <- runif(1, -h / 2, h / 2)
34
+
log_cur <- logpdf(theta)
35
+
log_shifted <- logpdf(theta + s * rho)
36
+
theta0 <- if (
37
+
(log_shifted - log_cur >= 0) || (runif(1) < exp(log_shifted - log_cur))
38
+
)
39
+
theta + s * rho else theta
40
+
41
+
# helpers
42
+
log_eps_h <- log(epsilon) + log(h)
43
+
orbit_points <- list(theta0)
44
+
log_vals <- logpdf(theta0)
45
+
46
+
# bookkeeping to get ends
47
+
left <- right <- 1
48
+
49
+
# Orbit selection procedure
50
+
B <- sample(c(FALSE, TRUE), M, replace = TRUE)
51
+
for (k in seq_len(M)) {
52
+
B_k <- B[k]
53
+
n_ext <- 2^(k - 1)
54
+
orbit_ext <- lapply(
55
+
seq_len(n_ext),
56
+
\(i)
57
+
(if (B_k) orbit_points[[right]] else orbit_points[[left]]) +
58
+
(2 * B_k - 1) * i * h * rho
59
+
)
60
+
log_ext <- sapply(orbit_ext, logpdf)
61
+
62
+
if (NURS_sub_stop(log_ext, log_eps_h)) break
63
+
64
+
if (B_k) {
65
+
orbit_points <- c(orbit_points, orbit_ext)
66
+
log_vals <- c(log_vals, log_ext)
67
+
} else {
68
+
orbit_points <- c(orbit_ext, orbit_points)
69
+
log_vals <- c(log_ext, log_vals)
70
+
left <- 1
71
+
}
72
+
right <- right + n_ext
73
+
74
+
if (NURS_stop(log_vals, log_eps_h)) break
75
+
}
76
+
77
+
orbit_points[[sample(
78
+
length(log_vals),
79
+
1,
80
+
prob = exp(log_vals - log_sum_exp(log_vals))
81
+
)]]
82
+
}
83
+
84
+
NURS <- function(logpdf, theta_init, n, epsilon, h, M) {
85
+
d <- length(theta_init)
86
+
draws <- matrix(NA_real_, n, d)
87
+
draws[1, ] <- theta_init
88
+
for (i in 2:n) {
89
+
draws[i, ] <- NURS_step(logpdf, draws[i - 1, ], epsilon, h, M)
90
+
}
91
+
draws
92
+
}
+42
README.md
+42
README.md
···
1
+
2
+
# R.NURS
3
+
4
+
Base R implementation of the No-Underrun Sampler.
5
+
6
+
## Installation
7
+
8
+
You can install the development version of R.NURS like so:
9
+
10
+
``` r
11
+
# install.packages("pak")
12
+
pak::pak("VisruthSK/R.NURS")
13
+
```
14
+
15
+
## Example
16
+
17
+
This is a basic example which shows you how to solve a common problem:
18
+
19
+
``` r
20
+
library(R.NURS)
21
+
library(ggplot2)
22
+
23
+
set.seed(0)
24
+
logpdf_funnel <- function(theta) {
25
+
y <- theta[1]
26
+
dnorm(y, 0, 3, log = TRUE) + sum(dnorm(theta[-1], 0, exp(y / 2), log = TRUE))
27
+
}
28
+
samples <- NURS(
29
+
logpdf_funnel,
30
+
theta_init = rep(0, 15),
31
+
n = 5000,
32
+
epsilon = 0.1,
33
+
h = 0.5,
34
+
M = 5
35
+
)
36
+
37
+
data.frame(y = samples[, 1], x1 = samples[, 2]) |>
38
+
ggplot(aes(x = y, y = x1)) +
39
+
geom_point(alpha = 0.3) +
40
+
theme_minimal()
41
+
```
42
+