如何在R语言中从自定义分布执行蒙特卡洛模拟抽样

Got it, let's break this down step by step. First, let's clarify what your custom distribution actually is: you've defined an unnormalized probability density function (PDF) with this line:

points <- Rmax * dexGAUS(x, mu = mu, sigma = sigma, nu = tau) * pgamma(x, shape = kappa, rate = rate)

This is a scaled product of the exponentially modified Gaussian (exGAUS) density and the Gamma cumulative distribution function (CDF), restricted to x ∈ [1, 20]. The key issue here is that points isn't a valid PDF yet—its integral over [1,20] doesn't equal 1. That's almost certainly why your first sampling attempt gave a histogram that didn't match expectations.

Below are two robust Monte Carlo methods to sample from this distribution, with code examples and explanations.

Method 1: Accept-Reject Sampling (Intuitive & Efficient for Simple Supports) #

Accept-Reject is a great starting point for distributions with a finite support (like your [1,20] range). Here's how to implement it:

Step 1: Formalize the Unnormalized PDF #

First, wrap your distribution in a reusable function:

library(gamlss)

# Your original parameters
mu <- 1
sigma <- 2
tau <- 3
kappa <- 3
rate <- 1
Rmax <- 20

# Define the unnormalized PDF
f <- function(x) {
  Rmax * dexGAUS(x, mu = mu, sigma = sigma, nu = tau) * pgamma(x, shape = kappa, rate = rate)
}

Step 2: Set Up the Proposal Distribution #

We'll start with a simple uniform proposal over [1,20] (easy to implement). Later, we'll show a more efficient proposal using the exGAUS distribution (since your target has an exGAUS component).

First, calculate the maximum value of f(x) (needed for the acceptance criterion):

x_grid <- seq(1, 20, 0.01)
f_max <- max(f(x_grid))

Step 3: Implement the Accept-Reject Sampler #

accept_reject_uniform <- function(n_samples) {
  samples <- numeric(n_samples)
  accepted <- 0
  
  while (accepted < n_samples) {
    # Draw a candidate from the uniform proposal
    x_candidate <- runif(1, min = 1, max = 20)
    # Calculate acceptance probability (uniform PDF = 1/(20-1))
    accept_prob <- f(x_candidate) / (f_max * (1/19))
    # Draw a uniform random number to decide acceptance
    u <- runif(1)
    
    if (u <= accept_prob) {
      accepted <- accepted + 1
      samples[accepted] <- x_candidate
    }
  }
  return(samples)
}

# Generate 1000 samples (set seed for reproducibility)
set.seed(123)
samples_uniform <- accept_reject_uniform(1000)

Step 4: Validate the Samples #

Compare the histogram to the normalized target PDF (divide f(x) by its integral over [1,20]):

# Calculate normalization constant
norm_const <- integrate(f, lower = 1, upper = 20)$value

# Plot
hist(samples_uniform, freq = FALSE, main = "Samples vs Target PDF", 
     xlab = "x", col = "lightgray")
lines(x_grid, f(x_grid)/norm_const, col = "red", lwd = 2)

Boost Efficiency with a Better Proposal #

The uniform proposal works but can have high rejection rates. Instead, use a truncated exGAUS distribution (matches the core of your target distribution):

# Truncated exGAUS proposal PDF
g <- function(x) {
  exgaus_cdf_range <- pexGAUS(20, mu = mu, sigma = sigma, nu = tau) - pexGAUS(1, mu = mu, sigma = sigma, nu = tau)
  dexGAUS(x, mu = mu, sigma = sigma, nu = tau) / exgaus_cdf_range
}

# Calculate max ratio of target to proposal
ratio <- f(x_grid)/g(x_grid)
M <- max(ratio)

# Efficient accept-reject sampler
accept_reject_exgaus <- function(n_samples) {
  samples <- numeric(n_samples)
  accepted <- 0
  exgaus_cdf_range <- pexGAUS(20, mu = mu, sigma = sigma, nu = tau) - pexGAUS(1, mu = mu, sigma = sigma, nu = tau)
  
  while (accepted < n_samples) {
    # Draw candidate from truncated exGAUS
    x_candidate <- rexGAUS(1, mu = mu, sigma = sigma, nu = tau)
    while (x_candidate < 1 || x_candidate > 20) {
      x_candidate <- rexGAUS(1, mu = mu, sigma = sigma, nu = tau)
    }
    # Acceptance probability
    accept_prob <- f(x_candidate)/(M * g(x_candidate))
    u <- runif(1)
    
    if (u <= accept_prob) {
      accepted <- accepted + 1
      samples[accepted] <- x_candidate
    }
  }
  return(samples)
}

set.seed(123)
samples_exgaus <- accept_reject_exgaus(1000)

Method 2: Metropolis-Hastings Sampling (Great for Complex Distributions) #

If accept-reject feels too slow, the Metropolis-Hastings algorithm (a type of MCMC) is a solid alternative. It uses a Markov chain to converge to your target distribution:

mh_sampler <- function(n_samples, init_val = 10, burn_in = 1000) {
  samples <- numeric(n_samples + burn_in)
  current <- init_val
  
  for (i in 1:(n_samples + burn_in)) {
    # Propose a new value from a normal distribution centered at current
    candidate <- rnorm(1, mean = current, sd = 2)
    # Reject candidates outside [1,20]
    if (candidate < 1 || candidate > 20) {
      samples[i] <- current
      next
    }
    # Calculate acceptance probability (symmetric proposal, so ratio simplifies)
    accept_prob <- min(1, f(candidate)/f(current))
    u <- runif(1)
    
    if (u <= accept_prob) {
      current <- candidate
    }
    samples[i] <- current
  }
  # Discard burn-in samples (chain needs time to converge to target)
  return(samples[(burn_in + 1):(burn_in + n_samples)])
}

set.seed(123)
samples_mh <- mh_sampler(1000)

# Validate
hist(samples_mh, freq = FALSE, main = "MH Samples vs Target PDF", 
     xlab = "x", col = "lightblue")
lines(x_grid, f(x_grid)/norm_const, col = "darkblue", lwd = 2)

Why Your First Attempt Failed #

Chances are you tried to sample directly from the exGAUS distribution and then modify it with pgamma—but pgamma is a CDF, not a valid transformation for generating random variables. You need to treat your custom expression as an unnormalized PDF and use proper Monte Carlo techniques to sample from it.

内容的提问来源于stack exchange，提问作者Alex