基于R语言optim()实现自定义分布的MLE估计方法问询

问题背景 #

我们可以编写通用函数生成指定分布的随机变量：

random_variates <- function(n, distribution, ...){
  distribution(n, ...)
}

set.seed(123)

# 负二项分布（2个参数）
random_variates(n = 5, distribution = rnbinom, mu = 5, size = 0.5)
# [1]  1  3  0 34  3

# 泊松分布（1个参数）
random_variates(n = 5, distribution = rpois, lambda = 5)
# [1] 9 8 6 7 1

但使用optim()进行极大似然估计（MLE）时，每个分布都需要编写专属的对数似然函数：

# 负二项分布的对数似然函数
nb_ll <- function(params, data){
  -sum(log(dnbinom(data, mu = params[1], size = params[2])))
}

# 示例数据
z_vals <- c(5, 0, 0, 3, 0, 0, 7, 0, 0, 0, 1, 0, 0, 0, 0, 1, 2, 0, 0, 
            3, 2, 0, 0, 1, 2, 0, 0, 2, 3, 0, 0, 0, 2, 1, 0, 2, 0, 0)

# MLE估计结果
optim(c(0.01, 0.01), nb_ll, data = z_vals)[["par"]]
# [1] 0.973937 0.476989

现在希望编写一个通用的似然函数，无需为每个分布单独实现，类似以下无法运行的代码：

# 这段代码无法运行
any_ll <- function(params, distribution, data, ...){
  -sum(log(distribution(data, ...)))
}
optim(c(0.01), any_ll, data = z_vals, distribution = dpois)
optim(c(0.01, 0.01), any_ll, data = z_vals, distribution = dnbinom)

请问在R语言中是否有实现该需求的方法？

解决方案 #

当然可以实现，核心是通过动态参数映射解决不同分布参数数量、名称不一致的问题，以下是两种可行方法：

方法1：指定参数名的通用对数似然函数 #

通过显式指定参数名，将params向量与分布的参数绑定，结合do.call动态调用分布函数：

any_ll <- function(params, distribution, data, param_names, ...) {
  # 将参数向量与参数名绑定为列表
  param_list <- setNames(as.list(params), param_names)
  # 合并数据、待估参数与固定参数（如log=TRUE）
  all_args <- c(list(x = data), param_list, list(...))
  # 返回负对数似然值（适配optim的最小化逻辑）
  -sum(do.call(distribution, all_args))
}

使用示例

# 泊松分布MLE估计：指定参数名为lambda
optim(c(0.01), any_ll, 
      data = z_vals, 
      distribution = dpois, 
      param_names = "lambda",
      log = TRUE)
# $par
# [1] 0.9736842

# 负二项分布MLE估计：指定参数名为mu和size
optim(c(0.01, 0.01), any_ll, 
      data = z_vals, 
      distribution = dnbinom, 
      param_names = c("mu", "size"),
      log = TRUE)
# $par
# [1] 0.973937 0.476989

方法2：基于参数顺序的简化版 #

如果熟悉目标分布的参数顺序，可以直接按位置传递params，无需指定参数名：

any_ll_simple <- function(params, distribution, data, ...) {
  # 按顺序构造参数列表：数据 -> 待估参数 -> 固定参数
  all_args <- c(list(x = data), as.list(params), list(...))
  -sum(do.call(distribution, all_args))
}

使用示例

# 泊松分布：dpois参数顺序为x, lambda, ...
optim(c(0.01), any_ll_simple, 
      data = z_vals, 
      distribution = dpois,
      log = TRUE)

# 负二项分布：需按dnbinom的参数顺序传递（size在前，mu在后）
optim(c(0.01, 0.01), any_ll_simple, 
      data = z_vals, 
      distribution = dnbinom,
      log = TRUE)

关键说明 #

• do.call是实现动态函数调用的核心，能灵活处理不同分布的参数差异。
• 直接指定log=TRUE避免额外的log()调用，同时保证数值稳定性。
• optim()会将第一个参数params作为优化变量，其余参数通过...传递给似然函数，需确保这些参数正确传入。

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