在R中无需qnorm()，用hit-or-miss模拟求满足P(L>x)=0.05的x值

Got it, let's break this down. We need to find the value x where only 5% of a normal distribution (mean=0, std=100) lies above x—that's the 95th percentile of the distribution. We'll use hit-or-miss Monte Carlo sampling to estimate this value, no qnorm() allowed. Here's how to do it:

Step 1: Build a Hit-or-Miss Probability Estimator #

First, we need a helper function that takes a candidate x and returns the proportion of samples from N(0,100) that are greater than x. This is the core of the hit-or-miss logic: each sample is either a "hit" (greater than x) or "miss" (not greater than x).

# Function to estimate P(L > x) via Monte Carlo sampling
estimate_upper_prob <- function(x, n_samples = 100000) {
  # Generate n_samples from the target normal distribution
  samples <- rnorm(n_samples, mean = 0, sd = 100)
  # Return the proportion of samples that are "hits" (greater than x)
  mean(samples > x)
}

Step 2: Use Binary Search to Narrow Down x #

We'll use binary search to efficiently find the right x. This works because the upper tail probability P(L > x) decreases as x increases—so we can adjust our search range based on whether our current estimate is above or below the target 0.05.

# Set our target and precision parameters
target_prob <- 0.05
tolerance <- 1e-4  # How precise we want our x estimate (0.01 error max)
n_samples <- 100000  # More samples = more stable estimates, but slower

# Initialize a wide search range (covers nearly all of N(0,100))
low <- 0
high <- 300

# Iterate until our range is smaller than the tolerance
while (high - low > tolerance) {
  mid <- (low + high) / 2
  current_prob <- estimate_upper_prob(mid, n_samples)
  
  if (current_prob > target_prob) {
    # Too many hits—we need a larger x to reduce the proportion of samples above it
    low <- mid
  } else {
    # Too few hits—we need a smaller x
    high <- mid
  }
  
  # Optional: Uncomment to track progress
  # cat("Current mid:", round(mid, 2), "| Estimated P(L>x):", round(current_prob, 4), "\n")
}

# Final estimated x value
estimated_x <- (low + high) / 2
cat("Estimated x where P(L > x) = 0.05:", round(estimated_x, 2), "\n")

What's Happening Here? #

• Hit-or-Miss Logic: Every sample we generate is a trial. If it's greater than our current candidate x, that's a hit—counted towards the 5% we're targeting. The mean(samples > x) gives us the empirical estimate of the upper tail probability.
• Binary Search: Starting with a range we know covers the true value (0 to 300, since 3 standard deviations covers nearly all normal distribution data), we split the range in half each iteration. We adjust the range based on whether our current midpoint gives a probability higher or lower than 0.05.
• Precision: The loop stops when our possible x range is smaller than 1e-4, meaning our estimate is accurate to within 0.01.

Optional Verification #

If you want to check how close we are to the true value (even though we can't use qnorm() for the solution, we can use it to confirm):

# True 95th percentile for comparison
true_x <- qnorm(0.95, mean = 0, sd = 100)
cat("True x (from qnorm):", round(true_x, 2), "\n")

You'll see our simulated estimate is almost identical to the true value (~164.49).

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