1

我正在尝试研究专门针对大量协变量 p 和高维状态的最大似然估计的分布（这意味着 p/n，样本大小为 n，约为 1/5）。我正在生成数据，然后使用statsmodels.api.Logit将参数拟合到我的模型中。

问题是，这似乎只适用于低维状态（如 300 个协变量和 40000 个观察值）。具体来说，我知道已经达到最大迭代次数，对数似然是 inf 即已经发散，并且出现“奇异矩阵”错误。

我不知道如何解决这个问题。最初，当我仍在使用较小的值（比如 80 个协变量，4000 个观察值）时，我偶尔会遇到这个错误，我设置了最多 70 次迭代而不是 35 次。这似乎有帮助。

然而，它现在显然无济于事，因为我的对数似然函数正在发散。这不仅仅是在最大迭代次数内不收敛的问题。

很容易回答这些包根本不打算处理这些数字，但是已经有专门研究这种高维机制的论文，比如在这里使用p=800协变量和n=4000观察值。

当然，本文使用的是 R 而不是 python。不幸的是，我不知道 R。但是我应该认为 python 优化应该具有可比的“质量”？

我的问题：

R 是否比 python statsmodels 更适合处理这种高 p/n 机制中的数据？如果是这样，为什么以及可以使用 R 的技术来修改 python statsmodels 代码？

如何修改我的代码以适用于 p=800 和 n=4000 左右的数字？

4

1 回答 1

2

在您当前使用的代码中（来自其他几个问题），您隐式使用了 Newton-Raphson 方法。这是sm.Logit模型的默认设置。它计算并反转 Hessian 矩阵以加快估计速度，但这对于大型矩阵来说代价非常高——更不用说当矩阵接近奇异时经常导致数值不稳定，正如您已经看到的那样。这在相关的 Wikipedia条目中有简要说明。

您可以通过使用不同的求解器来解决这个问题，例如bfgs（或lbfgs），就像这样，

model = sm.Logit(y, X)
result = model.fit(method='bfgs')

这对我来说运行得很好，即使是n = 10000, p = 2000.

除了估计之外，更有问题的是，您用于生成样本的代码会导致数据受到很大程度的准可分离性影响，在这种情况下，整个 MLE 方法充其量是有问题的。您应该立即对此进行调查，因为它表明您的数据可能不像您希望的那样表现良好。准可分离性在这里得到了很好的解释。

于 2019-07-21T12:47:44.787 回答