R语言蒙特卡洛模拟问题：已知X方差估计参数p的实现方法

Hey there! Let's walk through how to estimate ( p ) using R's random number generators, and address your question about the polynom package.

先理清方差的理论基础 #

First, let's confirm the variance formula since you already have an analytical solution. Since ( Z ) (standard normal) and ( U ) (uniform (0,1)) are independent, the variance of ( X = pZ + (1-p)U ) is:

( \text{Var}(X) = p^2\text{Var}(Z) + (1-p)^2\text{Var}(U) )

We know:

• ( \text{Var}(Z) = 1 )
• ( \text{Var}(U) = \frac{1}{12} \approx 0.0833 )

Plugging in the given ( \text{Var}(X) = 0.4 ), we get the quadratic equation:
p² + (1-p)²/12 = 0.4

Your analytical result of ( p \approx 0.62 ) is spot-on (the exact positive root is ~0.6231).

用R估计p的两种方法 #

方法1：数值求解解析方程（快速准确） #

We can use base R's uniroot() function to solve the equation directly, no extra packages needed:

# Define the function that equals zero when Var(X)=0.4
variance_equation <- function(p) {
  p^2 + (1 - p)^2 / 12 - 0.4
}

# Solve over the interval [0,1] (since p is a weight)
solution <- uniroot(variance_equation, interval = c(0, 1))
estimated_p <- solution$root

cat("Analytic numerical solution: p ≈", round(estimated_p, 4), "\n")

This will output p ≈ 0.6231, matching your analytical result.

方法2：用随机数生成的模拟估计（符合题目要求） #

If you want to use R's random number generators as specified, we can use simulation and optimization to find ( p ):

set.seed(123) # Ensure reproducibility
n <- 1e6 # Use a large sample for precision

# Generate samples of Z and U
Z_samples <- rnorm(n)
U_samples <- runif(n)

# Define a loss function: squared difference between sample variance of X and 0.4
loss_function <- function(p) {
  X <- p * Z_samples + (1 - p) * U_samples
  (var(X) - 0.4)^2
}

# Minimize the loss function over [0,1]
simulation_result <- optimize(loss_function, interval = c(0, 1))
sim_estimated_p <- simulation_result$minimum

cat("Simulation-based estimate: p ≈", round(sim_estimated_p, 4), "\n")

Running this will give you an estimate very close to 0.6231 (e.g., ~0.6229 or 0.6232 depending on the random seed).

关于polynom包：完全不需要！ #

You don't need to install or use the polynom package here. This quadratic equation is simple enough to handle with base R functions like uniroot() or optimize(). The polynom package is designed for more complex polynomial operations (like creating polynomial objects, fitting polynomials to data, etc.), which isn't necessary for this problem.

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