在R中基于IID随机变量最小值数据拟合其原始正态分布

Hey there! Since you're new to R and working on fitting the original normal distribution from a dataset of IID minima, let's break this down into manageable steps that you can follow along with.

First, let's recap: each value in your dataset is the minimum of n independent normal variables. We know the CDF of the minimum is $F_{Min}(x)=1-(1-F(x))^n$, where $F(x)$ is the CDF of a single normal variable $N(\mu, \sigma^2)$. Our goal is to estimate $\mu$ and $\sigma$.

To make this concrete, let's start by simulating some example data (you can replace this with your actual dataset later):

# Set a random seed for reproducibility
set.seed(123)

# True parameters of the original normal distribution
true_mu <- 10
true_sigma <- 2

# Number of variables used to compute each minimum
n <- 5

# Generate 100 minimum values (replace this with your real data!)
min_data <- replicate(100, min(rnorm(n, mean = true_mu, sd = true_sigma)))

To fit the parameters, we'll use Maximum Likelihood Estimation (MLE). First, we need to write the negative log-likelihood function (since R's optimization functions minimize by default, it's easier to work with negative log-likelihood instead of maximizing log-likelihood).

The PDF of the minimum value is the derivative of its CDF:
$f_{Min}(x) = n \cdot (1 - \Phi\left(\frac{x-\mu}{\sigma}\right))^{n-1} \cdot \frac{1}{\sigma} \cdot \phi\left(\frac{x-\mu}{\sigma}\right)$
where $\Phi$ is the standard normal CDF, and $\phi$ is the standard normal PDF.

Here's how to code this in R:

neg_log_likelihood <- function(params, data, sample_size) {
  mu <- params[1]
  sigma <- params[2]
  
  # Ensure sigma is positive (standard deviation can't be zero or negative)
  if (sigma <= 0) return(Inf)
  
  # Calculate terms for the log-PDF
  z_scores <- (data - mu) / sigma
  cdf_term <- pnorm(z_scores)
  log_pdf <- log(sample_size) + (sample_size - 1)*log(1 - cdf_term) + 
    dnorm(z_scores, log = TRUE) - log(sigma)
  
  # Return the sum of negative log-PDF values
  return(-sum(log_pdf))
}

Now we'll use R's optim() function to find the values of $\mu$ and $\sigma$ that minimize the negative log-likelihood. We need to start with a reasonable initial guess for the parameters (we'll use the mean and standard deviation of the minimum data as a starting point):

# Initial parameter guesses
init_guess <- c(mean(min_data), sd(min_data))

# Run the optimization (we use L-BFGS-B to enforce sigma > 0)
fit_results <- optim(
  par = init_guess,
  fn = neg_log_likelihood,
  data = min_data,
  sample_size = n,
  method = "L-BFGS-B",
  lower = c(-Inf, 1e-6)  # Lower bound for sigma: very small positive number
)

# Check the estimated parameters
cat("Estimated mu:", fit_results$par[1], "\n")
cat("Estimated sigma:", fit_results$par[2], "\n")

For our simulated data, this should give you values close to the true $\mu=10$ and $\sigma=2$.

If you want a more user-friendly interface for MLE, the bbmle package is a great option. It handles parameter bounds and summary output more cleanly:

# Install the package if you haven't already
# install.packages("bbmle")
library(bbmle)

# Define the negative log-likelihood (slightly simplified for bbmle)
nll <- function(mu, sigma) {
  if (sigma <= 0) return(Inf)
  z_scores <- (min_data - mu)/sigma
  -sum(log(n) + (n-1)*log(1 - pnorm(z_scores)) + dnorm(z_scores, log=TRUE) - log(sigma))
}

# Fit the model
mle_fit <- mle2(
  nll,
  start = list(mu = mean(min_data), sigma = sd(min_data)),
  lower = list(mu = -Inf, sigma = 1e-6),
  method = "L-BFGS-B"
)

# View detailed results
summary(mle_fit)

The summary will give you the estimated parameters, standard errors, and other useful stats.

• You must know the sample size n: Each minimum in your dataset must come from exactly n original variables. If n is unknown, you'd need to estimate it alongside $\mu$ and $\sigma$, which adds complexity.
• Initial guesses matter: If your initial parameter guess is way off, the optimizer might converge to a local minimum. Try a few different starting points if your results seem off.
• Check convergence: Look at the convergence value in the optim() output—if it's 0, the optimization converged successfully. For bbmle, the summary will tell you if convergence was achieved.

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