Rとpythonでパラメータの推計結果が違う

Question

前提・実現したいこと

Rとpythonでロジットモデルを作成しています。
しかし、RとPythonでは、パラメータの答えが違います。
何故かわかりません。
R:0.1864932 -0.5188502 -0.2187338 -2.3418396
python:0.25712579 -0.47485728 -0.2210386 -2.40078011

該当のソースコード

R:

データファイルの読み込み

Data <- read.csv("Test.csv",header=T)

データ数:Dataの行数を数える

hh <- nrow(Data)

Logit model の対数尤度関数の定義

fr <- function(x) {

パラメータの宣言

定数項

b1 <- x[1]
b2 <- x[2]

所要時間

d1 <- x[3]

女性ダミー

f1 <- x[4]

対数尤度のための変数を宣言€

LL = 0

効用の確定項の計算

            #定数項                        #所要時間              #女性ダミー

v1 <- b1matrix(1,nrow = hh,ncol=1) + d1( Data$所要時間1) + f1Data$女性ダミー
v2 <- b2matrix(1,nrow = hh,ncol=1) + d1*( Data$所要時間2)
v3 <- d1*( Data$所要時間3)

選択確率の計算

vehicle1 <- Data$選択1利用可能性 * exp(v1)
vehicle2 <- Data$選択2利用可能性 * exp(v2)
vehicle3 <- Data$選択3利用可能性 * exp(v3)

分母となる，各々のexp(V)の和をつくる

deno <- (vehicle1 + vehicle2 + vehicle3)

それぞれの選択確率を計算する

P1 <- vehicle1 / deno
P2 <- vehicle2 / deno
P3 <- vehicle3 / deno

選択確率が0になってしまった場合1にする(後で対数を取るため)

P1 <- (P1 != 0)*P1 + (P1 == 0)
P2 <- (P2 != 0)*P2 + (P2 == 0)
P3 <- (P3 != 0)*P3 + (P3 == 0)

実際の選択結果

C1 <- (Data$選択結果 == 1)
C2 <- (Data$選択結果 == 2)
C3 <- (Data$選択結果 == 3)

対数尤度の計算

LL <- colSums(C1log(P1) + C2log(P2) + C3*log(P3) )
}

対数尤度関数frの最大化#####

##パラメータの初期値を与える
b0 <- numeric(4)

##パラメータ値の最適化
res_BFGS <- optim(b0,fr, method = "BFGS", hessian = TRUE, control=list(fnscale=-1))

b <- res_BFGS$par
print("パラメータ")
print(b)

python:
from scipy.optimize import minimize
from numpy import log, exp, array, diagonal, sqrt, linalg
from pandas import read_csv
from numdifftools import Hessian

データファイル読み込み

c = read_csv('Test.csv')

対数尤度関数の定義

def func(Beta):
# 効用関数パラメータ
ASC_Car, ASC_Bus, B_TIME, B_FemaleC = Beta
# 選択結果判定指標
Choice = array([(c.Choice == 1), (c.Choice == 2), (c.Choice == 3)])
# 効用関数
V = array([
ASC_Car + B_TIME * c.TimeC + B_FemaleC * c.Female,
ASC_Bus + B_TIME * c.TimeB,
B_TIME * c.TimeR
])
# ロジットモデルの選択確率と尤度
PD = (array([c.AvailC, c.AvailB, c.AvailR]) * exp(V)).sum(axis = 0)
P = exp(V) / PD
IndLL = (Choice * log(P))
return -IndLL.sum()

パラメータの初期値を与える

Beta0 = [0,0,0,0]

最尤推定の実行

output = minimize(func, Beta0, method='BFGS', tol=1e-10, options={'disp': False })

パラメータの標準偏差・t値の計算

optBeta = output.x
fHessian = Hessian(func)
stdev = sqrt(diagonal(linalg.inv(fHessian(optBeta))))

結果の出力 （表記を修正）

print('Estimated Parameter:', optBeta)

Comments

[nandymak]

今のままでは読めないのでコードを「```」で囲ってください。RとPythonは分けてください。

プレビューで可読性があるか確認をお願いします。