如何在R语言的Rejection Sampling中控制样本量至指定数量？

Great question! Your current approach generates a fixed large set of samples upfront, which can waste computation if you only need a specific number of accepted draws. Let's adjust the code to stop sampling as soon as we hit your target number of accepted samples (e.g., 1000), and make it more efficient along the way.

Step-by-Step Solution #

First, we'll precompute key values to avoid redundant calculations, then use a loop that generates samples incrementally until we reach our target count.

Full Adjusted Code

# Define parameters for the truncated normal distribution
mu <- 2
sigma <- 2
truncation_min <- 1
target_sample_count <- 1000  # Your desired number of accepted samples
proposal_upper_bound <- 9

# Precompute the maximum density of the truncated normal distribution
# Since the original normal peaks at mu=2 (which is above our truncation min of 1),
# the peak of the truncated distribution remains at mu
max_target_density <- dnorm(mu, mu, sigma) / (1 - pnorm(truncation_min, mu, sigma))

# Initialize an empty vector to store accepted samples
accepted_samples <- c()

# Loop until we collect enough accepted samples
while (length(accepted_samples) < target_sample_count) {
  # Generate a single proposal from the uniform distribution [truncation_min, proposal_upper_bound]
  proposal <- runif(1, truncation_min, proposal_upper_bound)
  
  # Calculate the target density for this proposal
  target_density <- dnorm(proposal, mu, sigma) / (1 - pnorm(truncation_min, mu, sigma))
  
  # Compute acceptance probability and decide to accept/reject
  acceptance_probability <- target_density / max_target_density
  if (runif(1) < acceptance_probability) {
    accepted_samples <- c(accepted_samples, proposal)
  }
}

# Assign the final sample set to X as requested
X <- accepted_samples

# Optional: Verify the results with a histogram and theoretical curve
hist(X, freq = FALSE, col = "grey", breaks = 30, 
     xlim = c(truncation_min, proposal_upper_bound),
     main = paste("Truncated Normal Sample (", target_sample_count, " Accepted Draws)"),
     xlab = "Value")
curve(dnorm(x, mu, sigma)/(1-pnorm(truncation_min, mu, sigma)), 
      truncation_min, proposal_upper_bound, 
      add = TRUE, col = "red", lwd = 2)

Key Improvements #

1. Early Termination: The loop stops immediately once accepted_samples reaches target_sample_count, so you don't generate unnecessary extra samples.
2. Precomputed Max Density: Instead of calculating the maximum density from a large upfront sample, we compute it directly using the theoretical peak of the truncated distribution (since our truncation point is below the normal's mean, the peak stays at mu). This is more reliable and efficient.

Optional: Batch Sampling for Higher Efficiency #

If you're targeting a large number of samples, generating batches of proposals instead of single samples can reduce loop overhead:

# Batch-based approach (faster for large target counts)
accepted_samples <- c()
batch_size <- 200  # Adjust based on your needs

while (length(accepted_samples) < target_sample_count) {
  # Generate a batch of proposals
  proposals <- runif(batch_size, truncation_min, proposal_upper_bound)
  
  # Calculate target densities for the entire batch
  target_densities <- dnorm(proposals, mu, sigma) / (1 - pnorm(truncation_min, mu, sigma))
  
  # Determine which samples to accept
  acceptance_probs <- target_densities / max_target_density
  accepted_indices <- runif(batch_size) < acceptance_probs
  
  # Add accepted samples to our collection
  accepted_samples <- c(accepted_samples, proposals[accepted_indices])
}

# Trim to exactly the target count (in case we overshot)
accepted_samples <- accepted_samples[1:target_sample_count]
X <- accepted_samples

This batch method cuts down on the number of loop iterations, making the code run faster for larger target_sample_count values.

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