statsmodels与pymc中的对数似然

Question

我正在尝试在Python中进行我的第一次最大似然估计。其中一个步骤需要我计算模型参数的可能性。我找到了一些示例数据，可以简单总结如下：

import numpy as np
import pandas as pd
life_test = pd.DataFrame(columns=['points', 'time'])
life_test['points'] = np.linspace(1,14,14)
life_test['time'] = np.concatenate((np.linspace(5,40,8), np.linspace(50,100,6)), axis=0)

如果我通过statsmodels.api运行一个简单的模型，我从结果总结中得到了一个值：-14.601。

import statsmodels.api as sm
endog=np.array(life_test['points'])
exog=np.array(life_test['time'])
exog = sm.add_constant(exog)
results = sm.OLS(endog, exog).fit()
results.summary()

在OLS的源代码中查看，似乎这是对对数似然的基本计算。

params = np.array(results.params)
nobs2=results.nobs/2.0 # decimal point is critical here!
-nobs2*np.log(2*np.pi)-nobs2*np.log(1.0/(2*nobs2) *\
    np.dot(np.transpose(endog - np.dot(exog, params)),\
    (endog - np.dot(exog,params)))) - nobs2

当我尝试用PyMC实现这个时，我得到了不同的结果。可能是我在计算位置（loc）和尺度（scale）时出错了。

import pymc.distributions as dist
mu = exog.mean()
sigma = exog.std()
dist.normal_like(exog, mu, 1/sigma**2)

在这里，我得到了一个值：-135.29。我觉得我一定是在计算我的尺度和位置值时出错了，但也可能是我实现中的其他错误。也许OLS使用的是除了普通对数似然以外的其他可能性？我对statsmodels、PyMC和最大似然估计都还很陌生。有没有人知道我在这里做错了什么？

Answer 1

你可以用下面的方式把 statsmodels 的结果和 sklearn 的结果进行比较：

>>> x=sklearn.linear_model.LinearRegression(fit_intercept=False).fit(exog,endog)
>>> x.coef_
array([ 1.45714286,  0.13428571])

这两个结果是可以相互比较的：

>>> sm.OLS(endog, exog).fit().params
array([ 1.45714286,  0.13428571])

结果是一致的。不过，另一方面，看起来你只是计算了把 gaussian 拟合到 exog 数据的可能性，这和 线性回归 是不一样的。

如果你想用 pymc 来重新创建 线性回归，你需要按照以下步骤进行：

• 定义你的模型的自由参数，并设置一些先验条件
• 将你的输入数据传入模型，并使用不同的自由参数值进行计算
• 最后，设置你的 Gaussian 可能性

所以，使用 pymc 的实现方式是：

life_test = pd.DataFrame(columns=['points', 'time'])
life_test['points'] = np.linspace(1,14,14)
life_test['time'] = np.concatenate((np.linspace(5,40,8), np.linspace(50,100,6)), axis=0)
endog=np.array(life_test['points'])
exog=np.array(life_test['time'])
alpha = pm.Normal('alpha', mu=0, tau=2)
beta = pm.Normal('beta', mu=0, tau=2)
sigma = pm.Uniform('sigma', lower=0, upper=1)
y_est = alpha + beta * exog
radon_like = pm.Normal('y', mu=y_est, tau=sigma, observed=True,value=endog)
model = dict(rand_like=radon_like,alpha=alpha,beta=beta,sigma=sigma)
S = pm.MCMC(model)
S.sample(iter=100000,burn=1000)
pm.Matplot.plot(S)

如果你按照下面的步骤计算对数可能性，你会得到使用 pm.normal_like 分布的接近结果：

>>> results = sm.OLS(endog, exog).fit()
>>> y_est = results.params[0] + results.params[1] * exog[:,1]
>>> pm.normal_like(endog, y_est, 1/np.sqrt(y_est.std()))
-19.348540432740464