Python中的似然比检验

Question

我在Python 2.7中计算似然比检验时遇到了困难。

我有两个模型和相应的似然值。如果这些模型密切相关，我认为比较模型L2是否优于模型L1的规则是查看-2 * log(L2/L1)。

然后我想找到对应于-2 * log(L2/L1)的p值，并将其与L2优于L1的显著性联系起来。以下是我目前的进展：

import numpy as np
from scipy.stats import chisqprob

L1 = 467400. # log(likelihood) of my 1st fit
L2 = 467414. # log(likelihood) of my 2nd fit

LR = -2. * np.log(L2 / L1) # LR = -5.9905e-05

p = chisqprob(LR, 1) # L2 has 1 DoF more than L1

print 'p: %.30f' % p # p = 1.000000000000000000000000000000

five_sigma = 1 - scipy.special.erf(5 / np.sqrt(2.))                  :-)
print '5 sigma: %.30f' % five_sigma

five_sigma_check = 1 - 0.999999426696856                             :-(
print 'Check  : %.30f' % five_sigma_check

然而，我遇到了两个问题：

• 我的p值为1，但我预期它应该接近于0。
• 当我使用带有:-)标记的行上的公式来找到五西格玛时，例如，它与文献中引用的值不同，该行用:-(突出显示。我的five_sigma_check值取自这里。

请问有人能提供任何建议吗？我对Python和统计学相对较新。

谢谢。

- Sean Mooney

---

L1和L2是似然度还是对数似然度？如果它们是对数似然度，那么在计算LR时就不应该对它们取对数。 - C_Z_

L1和L2是似然对数的日志记录。我明白你的意思，对日志记录再次求对数没有意义... - Sean Mooney

Answer 1

使用以下公式计算给定对数似然的似然比：

from scipy.stats.distributions import chi2
def likelihood_ratio(llmin, llmax):
    return(2*(llmax-llmin))


LR = likelihood_ratio(L1,L2)


p = chi2.sf(LR, 1) # L2 has 1 DoF more than L1

print 'p: %.30f' % p 

# p: 0.000000121315450836607258011741

- C_Z_

---

chisqprob在scipy中已经不存在，或者对我来说根本没有导入。scipy的其他功能都可以正常使用（使用版本scipy-1.0.0）。是否有其他方法可以替代？ - CuriousDude

我认为新的等效函数是 scipy.stats.chi2.sf() - 或许你可以更新你的回答？ - Bede Constantinides

Answer 2

添加一些理论，以便有人能够理解：

在这里，我们有两个模型（假设是嵌套的），我们想比较这些模型在解释数据方面的有效性。从上图中，让我们现在使用Python 3实现LR测试：

from scipy.stats import chi2 
ll_0, ll_1 = 467400, 467414 # given, the log-likelihoods of the nested models m_0, m_1
# log likelihood for m_0 (H_0) must be <= log likelihood of m_1 (H_1)
Λ = -2 * (ll_0 - ll_1)
print(Λ)
# 28.0
df = 1 # given the difference in dof
# compute the p-value
pvalue = 1 - chi2(df).cdf(Λ) # since Λ follows χ2
print(pvalue)
# 1.2131545083660726e-07

我们可以绘制图表并清晰地看到，在α=0.05的情况下，我们可以拒绝零假设，支持模型m1。

α, df = 0.05, 1
x = np.linspace(0, 30, 1000)
plt.plot(x, chi2(df).pdf(x), label='χ2')
plt.axvline(chi2(df).ppf(1-α), color='red', label='α=0.05')
plt.scatter(Λ, 0, color='green', s=50, label='Λ')
plt.legend()
plt.title('χ2({}) distribution'.format(df), size=20)

- Sandipan Dey