Introduction of SIR sampling

Sampling Importance Resampling (SIR) is a statistical approach used to generate samples from an unknown or difficult-to-sample probability distribution via importance density sampling, weight adjustment, and resampling. It is a type of Monte Carlo sampling.The algorithm contains two-step process: importance sampling and resampling.

Define the target and proposal distributions

# The target distribution define as a mixture of two normal    distributions 
 target_dist <- function(x) {
  dnorm(x, mean = 15, sd = 1) + dnorm(x, mean = 20, sd = 2)
 }

# The proposal distribution define as a normal distribution
 proposal_dist <- function(x) {
  dnorm(x, mean = 17, sd = 2)
 }

Set the number of particles and iterations

 n_particles <- 1000
 n_iter <- 20

Initialize particles

 particles <- rnorm(n_particles, mean = 17, sd = 2)

Initialize weights

 weights <- rep(1/n_particles, n_particles)

Initialize plot

 plot_densities <- function(iter) {
 plot(target_dist, xlim = c(0, 25), ylim = c(0, 0.5), main =  paste("Iteration", iter), xlab = "x", ylab = "Density",las=1)
 curve(proposal_dist, add = TRUE, col = "red", lty = 2)
 lines(density(particles), col = "blue")
 }

Loop through iterations

 for (i in 1:n_iter) {

 # Generate new particles from proposal distribution
  new_particles <- rnorm(n_particles, mean = 17, sd = 2)

 # Calculate importance weights
  importance_weights <- target_dist(new_particles) / proposal_dist(new_particles)

  # Normalize importance weights
  importance_weights <- importance_weights / sum(importance_weights)

  # Resample particles
  resample_index <- sample.int(n_particles, n_particles, replace = TRUE, prob =  importance_weights)
  particles <- new_particles[resample_index]

  # Update weights
  weights <- rep(1/n_particles, n_particles)

  # Plot results
  plot_densities(i)
  legend("topright",legend=c("Proposal_dist","Target_dis","Estimate_dis"),col=c("red","black","blue"),lty=c(2,1,1))
}

Compute estimate of target distribution

target_estimate <- sum(weights * target_dist(particles))