Pythonのstatsmodels.OLSとRのlmの対数尤度が合わない件

対数尤度が合わないとAICやBICといった値も必然的に合わなくなると思われるので、けっこう困った話になる。
今回はずれるという事象だけをメモしておく。気が向いたらなぜずれるのかを検証してみるかもしれない。

回帰式は例によってこれ。

Pythonで作ったデータをstatsmodels.OLSに突っ込むのとあわせて、csvファイルに出力して、Rで読んでlmに食わせる。
まずはPythonのコード。

# -*- coding: utf-8 -*-
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns
import statsmodels.api as sm

#サンプリング数
N=1000

#データ生成
np.random.seed(0)
x=np.random.normal(1,1,N)
eps=np.random.normal(0,1,N)
y=2+1.5*x+eps

np.savetxt("xy.csv", np.vstack((x,y)).T, delimiter=",")

#OLS(statsmodels)
X=sm.add_constant(x, prepend=True) #1だけが詰まった列を左端に追加
model = sm.OLS(y,X)
results = model.fit()

次にRのコード。

dir = "C:\\Users\\owner\\data\\"
data_file = "xy.csv"
data = read.csv(paste(dir,data_file,sep=""),header=F,as.is=T)

x=as.matrix(data[,1])
y=as.matrix(data[,2])

result=lm(y~x)
summary(result)
logLik(result)
AIC(result)
BIC(result)

Pythonの実行結果は下記の通り。

OLS① - WHAT' CHA GONNA DO FOR ME?
でやったときと使ったデータが同じなので、当然ながら完全に同じ結果が得られている。

In [225]: X=sm.add_constant(x, prepend=True) #1だけが詰まった列を左端に追加
     ...: model = sm.OLS(y,X)
     ...: results = model.fit()
     ...: print(results.summary())
                            OLS Regression Results                            
==============================================================================
Dep. Variable:                      y   R-squared:                       0.692
Model:                            OLS   Adj. R-squared:                  0.692
Method:                 Least Squares   F-statistic:                     2241.
Date:                Thu, 12 Feb 2015   Prob (F-statistic):          2.28e-257
Time:                        01:56:08   Log-Likelihood:                -1386.1
No. Observations:                1000   AIC:                             2776.
Df Residuals:                     998   BIC:                             2786.
Df Model:                           1                                         
==============================================================================
                 coef    std err          t      P>|t|      [95.0% Conf. Int.]
------------------------------------------------------------------------------
const          2.0431      0.043     47.942      0.000         1.959     2.127
x1             1.4691      0.031     47.339      0.000         1.408     1.530
==============================================================================
Omnibus:                        1.076   Durbin-Watson:                   2.074
Prob(Omnibus):                  0.584   Jarque-Bera (JB):                1.155
Skew:                           0.059   Prob(JB):                        0.561
Kurtosis:                       2.883   Cond. No.                         2.53
==============================================================================

次にRの実行結果がこちら。Pythonと同じデータを使って回帰式を推定しているため、パラメータ等は完全に一致している。ただ、対数尤度とAIC、BICはPythonのstatsmodels.OLSと比べると微妙にずれている…

> result=lm(y~x)
> summary(result)

Call:
lm(formula = y ~ x)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.99269 -0.67483  0.01292  0.62140  3.15789 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  2.04312    0.04262   47.94   <2e-16 ***
x            1.46909    0.03103   47.34   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.9686 on 998 degrees of freedom
Multiple R-squared:  0.6919,	Adjusted R-squared:  0.6916 
F-statistic:  2241 on 1 and 998 DF,  p-value: < 2.2e-16

> logLik(result)
'log Lik.' -1386.084 (df=3)
> AIC(result)
[1] 2778.168
> BIC(result)
[1] 2792.891

次のステップとしては、どちらが正しいか確認するため、自分で対数尤度を計算してみることか。
Rならdnormですぐやれるけど、Pythonで正規分布の密度関数計算するのはググらないとわからない…
なんとなく、不偏分散を計算する時の自由度の違いとかが原因な気がする…