R中Dirichlet-Multinomial分布MLE估计参数卡在边界问题排查

• 目标：构建Dirichlet-Multinomial分布的参数最大似然估计器
• 现象：单独校验自定义密度函数ddirichlet.multinom()时结果看似正常，但将其传入stats4::mle()执行优化时，算法始终将参数推到设定的取值下界，调整起始值也无法得到合理估计结果
• 初步怀疑方向：自定义负对数似然函数ll()逻辑错误，或是mle()参数配置有误，多次调试未定位根因

复现代码 #

library(tidyverse)

# 用已知参数的Dirichlet分布生成测试数据
ex_data <- 
  gtools::rdirichlet(500, c(37, 5, 13, 120)) %>%
  as_tibble(.name_repair = "universal") %>%
  rename_with(~str_replace(.x, "...", "x")) %>%
  add_column(n = round(rnorm(500, 750, 20))) %>%
  mutate(across(starts_with("x"), ~round(.x * n))) %>%
  select(-n)

ex_data
#> # A tibble: 500 × 4
#>       x1    x2    x3    x4
#>    <dbl> <dbl> <dbl> <dbl>
#>  1   153    23    81   498
#>  2   156    18    65   500
#>  3   208     6    40   479
#>  4   142    25    73   524
#>  5   181     9    54   497
#>  6   139    21    40   512
#>  7   184    24    49   495
#>  8   137    10    39   544
#>  9   153     9    44   519
#> 10   186    21    42   464
#> # … with 490 more rows

# 自定义Dirichlet-Multinomial密度函数
ddirichlet.multinom <- function(x, alpha, log = FALSE) {
  
  a0 <- sum(alpha)
  n <- sum(x)
  const <- lgamma(a0) + lgamma(n + 1) - lgamma(n + a0)
  
  pr <- vector(length = length(x))
  
  for (i in 1:length(x)) {
    pr[i] <- lgamma(x[i] + alpha[i]) - lgamma(alpha[i]) - lgamma(x[i] + 1)
  }
  
  prob <- const * prod(pr)
  
  if (log) prob else exp(prob)
  
}

# 负对数似然函数
ll <- function(a1, a2, a3, a4) {
  
  x1 <- ex_data$x1
  x2 <- ex_data$x2
  x3 <- ex_data$x3
  x4 <- ex_data$x4
  
  log_likelihood <- vector("double", length = length(x1))
  
  for (i in 1:length(x1)) {
    log_likelihood[i] <- 
      ddirichlet.multinom(
        c(x1[i], x2[i], x3[i], x4[i]),
        c(a1, a2, a3, a4),
        log = TRUE
      )
  }
  
  -sum(log_likelihood)
  
}

# 执行最大似然估计
stats4::mle(
  ll,
  method = "L-BFGS-B",
  start = c(37, 5, 13, 120),
  lower = c(1.0001, 1.0001, 1.0001, 1.0001)
)
#> 
#> Call:
#> stats4::mle(minuslogl = ll, start = c(37, 5, 13, 120), method = "L-BFGS-B", 
#>     lower = c(1.0001, 1.0001, 1.0001, 1.0001))
#> 
#> Coefficients:
#>     a1     a2     a3     a4 
#> 1.0001 1.0001 1.0001 1.0001

问题出在自定义密度函数ddirichlet.multinom()的概率合并逻辑：

• 代码中const和pr向量的所有计算都用了lgamma()函数，输出的全是对数尺度下的数值
• 对数尺度下计算联合概率，应该用加法（对应原始尺度的乘法），但原代码错误使用了const * prod(pr)做乘法+求积，完全不符合对数运算规则，返回的对数似然值是完全错误的
• 之前的交叉校验没有覆盖对数密度返回值的正确性校验，才误以为密度函数工作正常

将密度函数中合并概率的错误行替换为对数尺度正确的加和计算即可：

# 修复后的密度函数
ddirichlet.multinom <- function(x, alpha, log = FALSE) {
  
  a0 <- sum(alpha)
  n <- sum(x)
  const <- lgamma(a0) + lgamma(n + 1) - lgamma(n + a0)
  
  pr <- vector(length = length(x))
  
  for (i in 1:length(x)) {
    pr[i] <- lgamma(x[i] + alpha[i]) - lgamma(alpha[i]) - lgamma(x[i] + 1)
  }
  
  # 对数尺度下用加和合并，不是乘法
  prob <- const + sum(pr)
  
  if (log) prob else exp(prob)
}

修复后重新运行mle()，即可得到接近真实生成参数c(37,5,13,120)的估计结果，不会再出现参数被推到下界的问题。

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