R中多项式MLE的简单数值估计问题
Problem with simple numerical estimation for MLE of multinomial in R
我正在尝试建立多项式分布的简单数值 MLE 估计。
多项式有一个约束 - 所有单元格概率需要加起来为一。
通常具有此约束的方法是将其中一个概率重新表示为（1 - 其他概率之和）
然而，当我 运行 这样做时，我遇到了一个问题，因为在优化过程中，我可能有一个负值的对数。
有没有想过如何解决这个问题？我尝试使用另一个优化包（Rsolnp）并且它有效，但我试图让它与简单的默认 R optim 一起工作以避免 constrained/nonlinear 优化。
这是我的代码（我知道我可以通过分析得到这种特殊情况下的结果，但这是一个玩具示例，我的实际问题比这里更大）。
set.seed(1234)
test_data <- rmultinom(n = 1, size = 1000, prob = rep(1/4, 4))
N <- test_data
loglik_function <- function(theta){
output <- -1*(N[1]*log(theta[1]) + N[2]*log(theta[2]) + N[3]*log(theta[3]) + N[4]*log(1- sum(theta)))
return(output)
}
startval <- rep(0.1, 3)
my_optim <- optim(startval, loglik_function, lower = 0.0001, upper = 0.9999, method = "L-BFGS-B")
任何想法或帮助将不胜感激。谢谢
全神贯注：我知道你问过（约束）ML 估计，但是用贝叶斯方法 à la Stan/rstan 怎么样？如果不是 useful/missing 重点，我会删除它。
模型只有几行代码
library(rstan)
model_code <- "
data {
int<lower=1> K; // number of choices
int<lower=0> y[K]; // observed choices
}
parameters {
simplex[K] theta; // simplex of probabilities, one for every choice
}
model {
// Priors
theta ~ cauchy(0, 2.5); // weakly informative
// Likelihood
y ~ multinomial(theta);
}
generated quantities {
real ratio;
ratio = theta[1] / theta[2];
}
"
您可以看到使用 Stan 数据类型 simplex 在 theta 上实现单纯形约束是多么容易。在斯坦语中，simplex allows you to easily implement a probability (unit) simplex
其中 K 表示参数的数量（此处：选择）。
还要注意我们如何使用 generated quantities 代码块，根据参数（此处 theta[1] 和 theta[2]）计算派生量（此处 ratio）。由于我们可以访问所有参数的后验分布，因此计算派生量的分布是微不足道的。
然后我们将模型拟合到您的 test_data
fit <- stan(model_code = model_code, data = list(K = 4, y = test_data[, 1]))
并显示参数估计的摘要
summary(fit)$summary
# mean se_mean sd 2.5% 25%
#theta[1] 0.2379866 0.0002066858 0.01352791 0.2116417 0.2288498
#theta[2] 0.26 20013 0.0002208638 0.01365478 0.2358731 0.2526111
#theta[3] 0.2452539 0.0002101333 0.01344665 0.2196868 0.2361817
#theta[4] 0.2547582 0.0002110441 0.01375618 0.2277589 0.2458899
#ratio 0.9116350 0.0012555320 0.08050852 0.7639551 0.8545142
#lp__ -1392.6941655 0.0261794859 1.19050097 -1395.8297494 -1393.2406198
# 50% 75% 97.5% n_eff Rhat
#theta[1] 0.2381541 0.2472830 0.2645305 4283.904 0.9999816
#theta[2] 0.2615782 0.2710044 0.2898404 3822.257 1.0001742
#theta[3] 0.2448304 0.2543389 0.2722152 4094.852 1.0007501
#theta[4] 0.2545946 0.2638733 0.2822803 4248.632 0.9994449
#ratio 0.9078901 0.9648312 1.0764747 4111.764 0.9998184
#lp__ -1392.3914998 -1391.8199477 -1391.3274885 2067.937 1.0013440
以及显示 theta 参数的点估计和置信区间的图表
plot(fit, pars = "theta")
更新：使用 maxLik
的约束 ML 估计
您实际上可以使用 maxLik 库提供的方法实现受约束的 ML 估计。我发现有点"fiddly"，因为收敛似乎对起始值的变化和使用的优化方法很敏感。
对于它的价值，这是一个可重现的例子：
library(maxLik)
x <- test_data[, 1]
定义多项式分布的对数似然函数；我在此处包含了一个 if 语句，以防止 theta < 0 案例引发错误。
loglik <- function(theta, x)
if (all(theta > 0)) sum(dmultinom(x, prob = theta, log = TRUE)) else 0
我在这里使用 Nelder-Mead 优化方法来查找对数似然函数的最大值。这里重要的一点是 constraints 参数，它以等式 A theta + B = 0 的形式实现约束，详情和示例请参见 ?maxNM。
res <- maxNM(
loglik,
start = rep(0.25, length(x)),
constraints = list(
eqA = matrix(rep(1, length(x)), ncol = length(x)),
eqB = -1),
x = x)
我们可以检查结果
summary(res)
--------------------------------------------
Nelder-Mead maximization
Number of iterations: 111
Return code: 0
successful convergence
Function value: -10.34576
Estimates:
estimate gradient
[1,] 0.2380216 -0.014219040
[2,] 0.2620168 0.012664714
[3,] 0.2450181 0.002736670
[4,] 0.2550201 -0.002369234
Constrained optimization based on SUMT
Return code: 1
penalty close to zero
1 outer iterations, barrier value 5.868967e-09
--------------------------------------------
并确认估计总和确实等于 1（在准确范围内）
sum(res$estimate)
#[1] 1.000077
示例数据
set.seed(1234)
test_data <- rmultinom(n = 1, size = 1000, prob = rep(1/4, 4))
我正在尝试建立多项式分布的简单数值 MLE 估计。
多项式有一个约束 - 所有单元格概率需要加起来为一。
通常具有此约束的方法是将其中一个概率重新表示为（1 - 其他概率之和）
然而，当我 运行 这样做时，我遇到了一个问题，因为在优化过程中，我可能有一个负值的对数。
有没有想过如何解决这个问题？我尝试使用另一个优化包（Rsolnp）并且它有效，但我试图让它与简单的默认 R optim 一起工作以避免 constrained/nonlinear 优化。
这是我的代码（我知道我可以通过分析得到这种特殊情况下的结果，但这是一个玩具示例，我的实际问题比这里更大）。
set.seed(1234)
test_data <- rmultinom(n = 1, size = 1000, prob = rep(1/4, 4))
N <- test_data
loglik_function <- function(theta){
output <- -1*(N[1]*log(theta[1]) + N[2]*log(theta[2]) + N[3]*log(theta[3]) + N[4]*log(1- sum(theta)))
return(output)
}
startval <- rep(0.1, 3)
my_optim <- optim(startval, loglik_function, lower = 0.0001, upper = 0.9999, method = "L-BFGS-B")
任何想法或帮助将不胜感激。谢谢
全神贯注：我知道你问过（约束）ML 估计，但是用贝叶斯方法 à la Stan/rstan 怎么样？如果不是 useful/missing 重点，我会删除它。
模型只有几行代码
library(rstan) model_code <- " data { int<lower=1> K; // number of choices int<lower=0> y[K]; // observed choices } parameters { simplex[K] theta; // simplex of probabilities, one for every choice } model { // Priors theta ~ cauchy(0, 2.5); // weakly informative // Likelihood y ~ multinomial(theta); } generated quantities { real ratio; ratio = theta[1] / theta[2]; } "您可以看到使用 Stan 数据类型
simplex在theta上实现单纯形约束是多么容易。在斯坦语中，simplexallows you to easily implement a probability (unit) simplex其中 K 表示参数的数量（此处：选择）。
还要注意我们如何使用
generated quantities代码块，根据参数（此处theta[1]和theta[2]）计算派生量（此处ratio）。由于我们可以访问所有参数的后验分布，因此计算派生量的分布是微不足道的。然后我们将模型拟合到您的
test_datafit <- stan(model_code = model_code, data = list(K = 4, y = test_data[, 1]))并显示参数估计的摘要
summary(fit)$summary # mean se_mean sd 2.5% 25% #theta[1] 0.2379866 0.0002066858 0.01352791 0.2116417 0.2288498 #theta[2] 0.26 20013 0.0002208638 0.01365478 0.2358731 0.2526111 #theta[3] 0.2452539 0.0002101333 0.01344665 0.2196868 0.2361817 #theta[4] 0.2547582 0.0002110441 0.01375618 0.2277589 0.2458899 #ratio 0.9116350 0.0012555320 0.08050852 0.7639551 0.8545142 #lp__ -1392.6941655 0.0261794859 1.19050097 -1395.8297494 -1393.2406198 # 50% 75% 97.5% n_eff Rhat #theta[1] 0.2381541 0.2472830 0.2645305 4283.904 0.9999816 #theta[2] 0.2615782 0.2710044 0.2898404 3822.257 1.0001742 #theta[3] 0.2448304 0.2543389 0.2722152 4094.852 1.0007501 #theta[4] 0.2545946 0.2638733 0.2822803 4248.632 0.9994449 #ratio 0.9078901 0.9648312 1.0764747 4111.764 0.9998184 #lp__ -1392.3914998 -1391.8199477 -1391.3274885 2067.937 1.0013440以及显示
theta参数的点估计和置信区间的图表plot(fit, pars = "theta")
更新：使用 maxLik
的约束 ML 估计
您实际上可以使用 maxLik 库提供的方法实现受约束的 ML 估计。我发现有点"fiddly"，因为收敛似乎对起始值的变化和使用的优化方法很敏感。
对于它的价值，这是一个可重现的例子：
library(maxLik)
x <- test_data[, 1]
定义多项式分布的对数似然函数；我在此处包含了一个 if 语句，以防止 theta < 0 案例引发错误。
loglik <- function(theta, x)
if (all(theta > 0)) sum(dmultinom(x, prob = theta, log = TRUE)) else 0
我在这里使用 Nelder-Mead 优化方法来查找对数似然函数的最大值。这里重要的一点是 constraints 参数，它以等式 A theta + B = 0 的形式实现约束，详情和示例请参见 ?maxNM。
res <- maxNM(
loglik,
start = rep(0.25, length(x)),
constraints = list(
eqA = matrix(rep(1, length(x)), ncol = length(x)),
eqB = -1),
x = x)
我们可以检查结果
summary(res)
--------------------------------------------
Nelder-Mead maximization
Number of iterations: 111
Return code: 0
successful convergence
Function value: -10.34576
Estimates:
estimate gradient
[1,] 0.2380216 -0.014219040
[2,] 0.2620168 0.012664714
[3,] 0.2450181 0.002736670
[4,] 0.2550201 -0.002369234
Constrained optimization based on SUMT
Return code: 1
penalty close to zero
1 outer iterations, barrier value 5.868967e-09
--------------------------------------------
并确认估计总和确实等于 1（在准确范围内）
sum(res$estimate)
#[1] 1.000077
示例数据
set.seed(1234)
test_data <- rmultinom(n = 1, size = 1000, prob = rep(1/4, 4))