# Custom Density Functions

To sample from a custom density function, create said function, then use it as the custom_density parameter with one of the samplers, instead of giving distr_name, distr_params (see Supported Distributions and How to Sample) and weights, if applicable.

In this example we create a 2D uniform distribution, which is not supported by Rcpp and RcppDist at the time of writing.

# Create function
customDensity_r <- function(x){
if (x[1] > 0 && x[1] < 3 && x[2] < 0 && x[2] > -3){
return (1)
} else {
return (0)                    }
}
# Sample
Y <- samplr::sampler_mh(start = c(0,0), custom_density = customDensity_r)
#> Warning in .checkSigmaProp(sigma_prop, start_dim): The variance of the proposal
#> distribution was not given and defaulted to the identity matrix

# Plot results
x <- Y[[1]][,1]
y <- Y[[1]][,2]
plot(x,y, xlim = c(-5,5), ylim = c(-5,5))

#### On Timing

Using a custom function, rather than sampling from a supported distribution, will considerably increase the running time (as the C++ code under the hood has to constantly switch between C++ and R). Creating a density function with Rcpp::cppFunction doesn’t change this much (it’s actually probably worse).

Rcpp::cppFunction("double customDensity_cpp(NumericVector x){
if (x(0) > 0 && x(0) < 3 && x(1) < 0 && x(1) > -3){
return 1;
} else {
return 0;
}
}")
X <- bench::mark(
"CPP Density" = {
samplr::sampler_mh(
start = c(0,0),
custom_density = customDensity_cpp,
sigma_prop = diag(2))
},
"R Density" = {
samplr::sampler_mh(
start = c(0,0),
custom_density = customDensity_r,
sigma_prop = diag(2))
},
"Supported Density" = {
samplr::sampler_mh(
start = c(0,0),
distr_name = "mvnorm",
distr_params = list(c(0,1), diag(2)),
sigma_prop = diag(2))
},
check = FALSE,
)
knitr::kable(as.data.frame(X[,c("expression", "min", "median")]))
expression min median
CPP Density 2.4ms 3.4ms
R Density 1.83ms 1.93ms
Supported Density 935.5µs 1.09ms

Note, however, than this is still much faster than doing everything in R.

bespoke_mh <- function(pdf, start, iterations=1024, sigma_prop){
# Initialize variables ---------------------------------
acceptances <- 0
n_dim <- length(start)
chain <- matrix(0, nrow=iterations, ncol = n_dim)
chain[1,] <- start
probabilities <- c(pdf(start))

# Run the sampler ------------------------------------------------

for (i in 2:iterations){
current_x <- chain[i-1,]
# generate proposal
proposal <- mvtnorm::rmvnorm(1, mean = current_x, sigma = sigma_prop)
# calculate current and proposal probabilities
prob_curr <- probabilities[i-1]
prob_prop <- pdf(proposal)
accept <- FALSE

# proposal is accepted with probability prob_prop / prob_curr
if (prob_curr != 0){
ratio <- prob_prop / prob_curr
if (ratio >= 1){
accept <- TRUE
# The beta parameter (temperature), beta <= 1,
# increases the value of the ratio making hotter
# chains more likely to accept proposals
} else if (stats::runif(1) < ratio){
accept <- TRUE
}
} else {
# in case the current probability is 0
# (in which case the ratio cannot be calculated),
# the step is accepted if the probability of the proposal is not 0
if (prob_prop > 0) {
accept <- TRUE
}
}

if (accept){
chain[i,] <- proposal
acceptances <- acceptances +  1
probabilities[i] <- prob_prop
} else {
chain[i,] <- current_x
probabilities[i] <- prob_curr
}

}

return(list(
chain = chain,
acceptance_ratio = acceptances/(iterations - 1)))
}

X <- bench::mark(
"Bespoke MH" = {
bespoke_mh(
customDensity_r,
c(0,0),
iterations = 1024,
sigma_prop = diag(2)
)
}
)
#> Warning: Some expressions had a GC in every iteration; so filtering is
#> disabled.
knitr::kable(as.data.frame(X[,c("expression", "min", "median")]))
expression min median
Bespoke MH 113ms 115ms