stors_package.R
In stors: Step Optimised Rejection Sampling

## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----setup, echo = TRUE-------------------------------------------------------

library(stors)

# Optimize a proposal for the Normal distribution with 4091 steps
proposal <- srnorm_optimize(steps = 4091)

# Visualize the optimized proposal
plot(proposal)

## ----echo = TRUE--------------------------------------------------------------

# Generating 10 samples from the standard Normal distribution
standard_normal_samples <- srnorm(10, mean = 2, sd = 1)
print(standard_normal_samples)


## ----echo = TRUE--------------------------------------------------------------

# Optimize a custom proposal for N(3, 2)
proposal <- srnorm_optimize(mean = 3, sd = 2, steps = 2000)

# Generate a sample of size 10^3 from N(3, 2)
sample <- srnorm_custom(10^3)

# Visualize the generated sample
hist(sample, main = "Histogram of Samples from N(3, 2)", xlab = "Value", col = "skyblue", border = "white")



## ----echo = TRUE--------------------------------------------------------------

# Define the truncation bounds
lower_bound <- -1
upper_bound <- 1

# Optimize the scalable proposal with truncation bounds
proposal <- srnorm_optimize(xl = lower_bound, xr = upper_bound, steps = 4091)

# Generate samples from the truncated standard Normal distribution
sample <- srnorm(10^3)
hist(sample, main = "Truncated Standard Normal (-1, 1)", xlab = "Value", col = "lightblue", border = "white")

# Generate samples from the truncated Normal distribution N(2, 1)
sample <- srnorm(10^3, mean = 2)
hist(sample, main = "Truncated Normal N(2, 1) (-1, 1)", xlab = "Value", col = "lightgreen", border = "white")



## ----echo = TRUE--------------------------------------------------------------

# Define the truncation bounds
lower_bound <- 2
upper_bound <- 6

# Optimize a custom proposal with truncation bounds and mean = 4
proposal <- srnorm_optimize(mean = 4, xl = lower_bound, xr = upper_bound, steps = 4091)

# Generate samples from the truncated Normal distribution
sample <- srnorm_custom(10^3)
hist(sample, main = "Truncated Custom Normal (2, 6)", xlab = "Value", col = "lightblue", border = "white")




## ----echo = TRUE--------------------------------------------------------------

# Optimize a proposal for the Laplace distribution
proposal <- srlaplace_optimize(steps = 4091)

# View the details of the optimized proposal
print(proposal)

# Visualize the proposal
plot(proposal)


## ----echo = TRUE--------------------------------------------------------------

# Optimize a proposal for the Normal distribution
proposal <- srlaplace_optimize(steps = 1020)

# Visualize the entire proposal
plot(proposal, main = "Full Proposal Visualization")

# Focus on a specific region of the proposal (e.g., x-axis range -2 to 2)
plot(proposal, x_min = -1, x_max = 1, main = "Zoomed-In Proposal (-1, 1)")


## ----echo = TRUE--------------------------------------------------------------
# Customizing the proposal for a specific sample size
custom_proposal <- srnorm_optimize(target_sample_size = 10000)

print(custom_proposal)

plot(custom_proposal)

# Customizing the proposal with a specific pre-acceptance threshold
custom_proposal_high_theta <- srnorm_optimize(theta = 0.9)

print(custom_proposal_high_theta)

plot(custom_proposal_high_theta)

# Customizing the proposal within a specific range
custom_proposal_range <- srnorm_optimize(proposal_range = c(-1, 1))

print(custom_proposal_range)

plot(custom_proposal_range)


## ----echo = TRUE--------------------------------------------------------------
# Customizing the proposal with a specific number of steps
custom_proposal_steps <- srnorm_optimize(steps = 50)

plot(custom_proposal_steps)

custom_proposal_steps <- srnorm_optimize(steps = 500)

plot(custom_proposal_steps)


## ----echo = FALSE-------------------------------------------------------------

# Delete all non-symmetric proposals for the Normal distribution
delete_built_in_proposal(sampling_function = "srnorm", proposal_type = "scaled")
delete_built_in_proposal(sampling_function = "srnorm", proposal_type = "custom")

# Create a symmetric proposal for the standard Normal distribution
proposal <- srnorm_optimize(symmetric = TRUE, steps = 2040)

# Visualize the symmetric proposal
plot(proposal)


## ----echo = FALSE-------------------------------------------------------------

# Delete symmetric proposals for the Normal distribution
delete_built_in_proposal(sampling_function = "srnorm", proposal_type = "scaled")

# Optimize the proposal for a truncated standard Normal distribution
proposal <- srnorm_optimize(mean = 0, sd = 1, xl = -1, xr = 1, steps = 4020)

# Visualize the truncated proposal
plot(proposal, main = "Truncated Standard Normal Proposal (-1, 1)", xlab = "Value", col = "lightblue")

# Generate 10^4 samples from the truncated distribution
sample <- srnorm_custom(10^4)

# Visualize the generated samples
hist(sample, main = "Histogram of Samples from Truncated Standard Normal", xlab = "Value", col = "lightgreen", border = "white")


## ----echo = FALSE-------------------------------------------------------------
name <- names(stors:::built_in_proposals)
formatted_names <- sprintf("{.fun %s_custom}", name)
cli::cli_ul(formatted_names)

## ----echo = TRUE--------------------------------------------------------------

modes_bi <- c(0.00134865, 3.99865)

f_bi <- function(x) {
  0.5 * (sqrt(2 * pi))^(-1) * exp(-(x^2) / 2) + 0.5 * (sqrt(2 * pi))^(-1) * exp(-((x - 4)^2) / 2)
}

h_bi <- function(x) log(f_bi(x))


h_prime_bi <- function(x) {
  (-(exp(-1 / 2 * (-4 + x)^2) * 0.5 * (-4 + x)) / sqrt(2 * pi) - (exp(-x^2 / 2) * 0.5 * x) / sqrt(2 *  pi)) / ((exp(-x^2 / 2) * 0.5) / sqrt(2 * pi) + (exp(-1 / 2 * (-4 + x)^2) * 0.5) / sqrt(2 * pi))
}

bi_proposal <- build_proposal(lower = -Inf, upper = Inf, modes = modes_bi, f = f_bi, h = h_bi, h_prime = h_prime_bi)

plot(bi_proposal)


## ----echo = TRUE--------------------------------------------------------------

# Printing the properties of the custom proposal
print(bi_proposal)

# Visualizing the custom proposal
plot(bi_proposal)


## ----echo = TRUE--------------------------------------------------------------

# Creating a sampler function using the custom proposal
bi_sampler <- build_sampler(bi_proposal)

plot(bi_proposal)

# Generating samples from the target distribution
sample <- bi_sampler(10^3)

hist(sample, main = "Sample Generated From Multi-modal Distribution ", xlab = "Value", col = "lightblue", border = "white")



## ----echo = TRUE--------------------------------------------------------------

# Define modes for the multi-modal distribution
modes_bi <- c(0.00134865, 3.99865)

# Build a truncated proposal for the multi-modal distribution
bi_trunc_proposal <- build_proposal(
  lower = -1,
  upper = 6,
  modes = modes_bi,
  f = f_bi,
  h = h_bi,
  h_prime = h_prime_bi,
  steps = 2040
)

# Visualize the truncated proposal
plot(bi_trunc_proposal, main = "Truncated Multi-Modal Proposal (-1, 6)", xlab = "Value", col = "skyblue")

# Create a sampler for the truncated distribution
bi_sampler_trunc_sampler <- build_sampler(bi_trunc_proposal)

# Generate 10^3 samples from the truncated distribution
sample <- bi_sampler_trunc_sampler(10^3)

# Generate 10^3 samples from the truncated distribution
sample <- bi_sampler_trunc_sampler(10^3)

# Visualize the generated samples
hist(
  sample,
  main = "Histogram of Samples from Truncated Multi-Modal Distribution",
  xlab = "Value",
  col = "lightgreen",
  border = "white"
)




## ----echo = TRUE--------------------------------------------------------------

save_proposal(bi_proposal, "bimodal_proposal")


## ----include = FALSE----------------------------------------------------------
for (i in 1:3) {
  save_proposal(bi_proposal, paste0("bimodal_proposal_", i))
}


## ----echo = TRUE--------------------------------------------------------------

bi_proposal <- load_proposal("bimodal_proposal")


## ----echo = TRUE--------------------------------------------------------------

print_proposals()


## ----echo = TRUE--------------------------------------------------------------

delete_proposal("bimodal_proposal")

print_proposals()


## ----include = FALSE----------------------------------------------------------

srnorm_optimize() # optimize the scaled non-symmetric proposal, since we have deleted it earlier

for (i in 1:3) {
  delete_proposal(paste0("bimodal_proposal_", i))
}

for (name in names(stors:::built_in_proposals)) {
  fun_name <- paste0(name, "_optimize")
  do.call(fun_name, list(steps = 4091))
}