两个二项分布混合模型的极大似然函数最优性

Got it, let's walk through how to estimate the parameters (n, p, a, b) from your mixed binomial distribution using R's JDEoptim package. Here's a practical, step-by-step implementation:

First, make sure you have the package installed. If not, run this:

install.packages("JDEoptim")
library(JDEoptim)

If you don't have your real dataset ready yet, let's simulate some data from the mixed binomial distribution to test our method. We'll use known true parameters so we can check if our estimation works:

# True parameters
true_n <- 10
true_p <- 0.3
true_a <- 0.2
true_b <- 0.7

# Generate 1000 observations
set.seed(123) # For reproducibility
n_obs <- 1000
# First, draw which component each observation comes from
component <- rbinom(n_obs, 1, true_p)
# Generate data from each component
x <- ifelse(component == 1, 
            rbinom(n_obs, true_n, true_a),
            rbinom(n_obs, true_n, true_b))

Since JDEoptim minimizes functions by default, we'll define the negative log-likelihood (maximizing the log-likelihood is equivalent to minimizing its negative).

Note: We'll treat (n) as a continuous value during optimization first, then round it to the nearest integer at the end (since (n) must be a positive integer for the binomial distribution). We also add small buffer values to (p, a, b) to avoid log(0) errors.

neg_log_likelihood <- function(params, x) {
  n <- params[1]
  p <- params[2]
  a <- params[3]
  b <- params[4]
  
  # Enforce valid parameter ranges (avoid numerical issues)
  if (n < 1 || p <= 0 || p >= 1 || a <= 0 || a >= 1 || b <= 0 || b >= 1) {
    return(Inf) # Penalize invalid parameter combinations
  }
  
  # Calculate log-likelihood for each observation
  ll <- log(p * dbinom(x, size = n, prob = a) + (1 - p) * dbinom(x, size = n, prob = b))
  
  # Return negative sum of log-likelihoods (since we minimize)
  return(-sum(ll))
}

Next, set reasonable bounds for each parameter:

• (n): Must be at least the maximum value in your dataset (you can't have more successes than trials). We'll set lower bound to 1, upper bound to max(x) + 10 (gives some room for estimation).
• (p, a, b): Must be between 0 and 1. We use 0.001 and 0.999 instead of strict 0/1 to avoid numerical instability.

Then call JDEoptim:

# Set parameter bounds: c(n_lower, p_lower, a_lower, b_lower), c(n_upper, p_upper, a_upper, b_upper)
lower_bounds <- c(1, 0.001, 0.001, 0.001)
upper_bounds <- c(max(x) + 10, 0.999, 0.999, 0.999)

# Run optimization
optim_result <- JDEoptim(lower = lower_bounds,
                         upper = upper_bounds,
                         fn = neg_log_likelihood,
                         x = x) # Pass our dataset as an argument

Pull out the estimated parameters and compare them to the true values (if using simulated data):

# Extract optimized parameters
est_params <- optim_result$par
# Round n to nearest integer
est_n <- round(est_params[1])
est_p <- est_params[2]
est_a <- est_params[3]
est_b <- est_params[4]

# Print results
cat("Estimated Parameters:\n")
cat(sprintf("n: %d (True: %d)\n", est_n, true_n))
cat(sprintf("p: %.4f (True: %.4f)\n", est_p, true_p))
cat(sprintf("a: %.4f (True: %.4f)\n", est_a, true_a))
cat(sprintf("b: %.4f (True: %.4f)\n", est_b, true_b))

When you run this with the simulated data, you should get estimates very close to the true values. For your real dataset, just replace the simulated x with your actual data.

Quick Notes for Your Real Dataset #

• If your dataset has a large maximum value, adjust the upper bound for (n) accordingly.
• If the optimization converges slowly, tweak JDEoptim's control parameters (like maxiter or popsize)—check the package documentation for details.
• Since (n) is discrete, do a final check around the rounded estimated (n) (e.g., test (n = est_n -1, est_n, est_n +1)) to find the exact integer that gives the highest log-likelihood.

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