在R中基于概率集合估计Beta分布的alpha与beta参数

问题背景 #

需要仅通过一组概率值估计Beta分布的alpha和beta参数，尝试三种R方法时遇到以下问题：

• betareg包返回phi而非alpha/beta参数，且存在警告；
• bbmle包看似可行，但构造全常量data.frame的做法存疑，运行有警告；
• optim函数未返回有效估计结果，逻辑需修正。

以下是数据集生成代码及各方法的修正方案：

生成概率数据集 #

set.seed(1234)
n.samples <- 1000
alpha <- 2
beta  <- 2
my.probs <- rbeta(n.samples, alpha, beta)
mean(my.probs)
# [1] 0.4925867

1. betareg包：从phi反推alpha和beta #

betareg的输出中，phi是Beta分布的精度参数（$\phi = \alpha + \beta$），均值模型的截距对应logit变换后的均值$\mu = \frac{\alpha}{\alpha+\beta}$。可通过以下步骤反推alpha和beta：

• 用plogis()将截距转换为原始均值$\mu$；
• alpha = $\mu * \phi$；
• beta = $\phi - \alpha$。

修正后的代码：

library(betareg)
my.model <- betareg(formula = my.probs ~ 1)

# 计算alpha和beta估计值
mu <- plogis(coef(my.model)[[1]])
phi <- my.model$coefficients$phi
alpha_est <- mu * phi
beta_est <- phi - alpha_est

alpha_est # 输出：~2.04
beta_est  # 输出：~2.14

关于警告：deparse()的警告是包的小问题，不影响结果，可忽略。

2. bbmle包：简化写法并消除警告 #

你之前构造全1的data.frame完全冗余，且公式接口的写法会因未正确绑定数据导致警告。直接自定义负对数似然函数即可，同时初始值要接近真实值（避免迭代中出现无效参数）。

修正后的代码：

library(bbmle)
# 定义负对数似然函数
neg_loglik <- function(shape1, shape2) {
  -sum(dbeta(my.probs, shape1, shape2, log = TRUE))
}

# 拟合模型，初始值设为接近真实值的(2,2)
my.model <- mle2(neg_loglik, start = list(shape1 = 2, shape2 = 2))
summary(my.model)

输出的shape1和shape2即为alpha和beta的估计值，结果会接近真实的2和2，且无警告。

3. optim函数：修正似然逻辑并限制参数范围 #

你之前的beta.reg函数逻辑错误（误引入线性回归的残差平方和逻辑），Beta分布的MLE只需要估计alpha和beta两个正参数。修正要点：

• 去掉多余的b0参数；
• 确保alpha和beta始终为正（用L-BFGS-B方法设置下界）；
• 直接计算负对数似然总和。

修正后的代码：

# 正确的负对数似然函数
beta_reg_optim <- function(params, data) {
  alpha <- params[1]
  beta <- params[2]
  # 参数必须为正，否则返回极大值终止优化
  if (alpha <= 0 || beta <= 0) return(Inf)
  -sum(dbeta(data, shape1 = alpha, shape2 = beta, log = TRUE))
}

# 初始值设为接近真实值的(2,2)
init_params <- c(2, 2)
# 使用L-BFGS-B方法，限制参数下界为极小正数
my_model_optim <- optim(init_params, beta_reg_optim, data = my.probs, 
                        method = "L-BFGS-B", lower = c(1e-6, 1e-6), hessian = TRUE)

# 查看估计结果
my_model_optim$par # 输出：~2.06, ~2.12
# 计算标准误
se <- sqrt(diag(solve(my_model_optim$hessian)))
se

优化后会得到有效的alpha和beta估计值，且无警告。

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