1

我正在尝试在 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 运行一个简单的模型。我从 results.summary() 中得到 -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 的值。我觉得我一定是误算了我的 scale 和 loc 值，但这可能是我的实现中的一些其他错误。除了正常的对数似然之外，也许 OLS 正在使用其他一些可能性？一般来说，我对 statsmodels、PyMC 和 MLE 还是很陌生。有谁知道我在这里做错了什么？

4

1 回答 1

3

您可以使用以下方法比较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])

结果是一致的。另一方面，您似乎只是计算了将 a 拟合gaussian到与exog不同的数据的可能性linear-reqression。

为了重新创建linear regression，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
于 2014-08-25T12:45:56.833 回答