1

Python で最初の最尤推定を実行しようとしています。この手順の 1 つで、モデル パラメーターの可能性を計算する必要があります。ここに要約できるサンプルデータをいくつか見つけました。

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 で実装しようとすると、異なる結果が得られます。私の側での場所とスケールの計算ミスかもしれません。

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

ここでは -135.29 という値を取得します。scale と loc の値を誤って計算しているに違いないと感じていますが、実装の他のエラーである可能性があります。おそらく、OLS は通常の対数尤度以外の他の尤度を使用していますか? 私は統計モデル、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])

結果は一貫しています。一方、 とは異なるデータにagaussianをフィッティングする可能性を計算したようです。exoglinear-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 に答える