結論は対数尤度関数のへシアン(Fisherの情報量行列)の逆行列の対角成分なんだけどソースコード追っかけるのに割と時間を食ったのでメモを残す。
設定
- statsmodels ver0.12.2
- statsmodels.discrete.discrete_model.Logitを利用したロジスティック回帰分析について調べました
要旨
statsmodels.discrete.discrete_model.Logit
の戻り値のsummaryメソッドを実行した時に表示される回帰係数の標準誤差はLogit.hess()
で計算される対数尤度関数のへシアンの逆行列の対角成分。
本題
Logitのfitを実行すると返されるのはstatsmodels.discrete.discrete_model.BinaryResultsWrapper
これのsummaryメソッドを実行した時に表示される回帰係数の標準誤差がどう計算されているか知りたかった
BinaryReuslts.summary()内でさらに
Super(DiscreteResults)のsummaryメソッドを呼び出している
そこでさらにstatsmodels.iolib.summaryのadd_table_paramsメソッドを呼び出して、さらに同じクラス内のsummary_paramsメソッドにとぶ。
そこでようやくstd_errにresults.bseを代入しているのを発見。
もう一度BinaryResultsクラスに戻るがこの中でも、親クラス(DescreteResults)でも標準誤差を計算している様子はない。のでさらに親クラス(LikelihoodModelResults)まで戻ってようやくbse()というメソッドを発見・・・・長かった。
bse()の中身は同じクラス内のメソッドcov_params()を呼び出していて、、、、まだ旅が続く
cov_params()、if分岐が多くてミスってる可能性あるけど、基本設定(特にへシアンとかの計算に自作の関数与えてない)では多分self.normalized_cov_params * scaleを返しているだけ。
normalized_cov_paramsは__init__
内で引数で与えられているので、、、ようやくコード内で、BinaryResults(DescreteResults(LikelihoodModelResults))
にnormalized_cov_paramsを渡している箇所を見ればいいことにたどり着く。
元々BinaryResulutsはLogit.fit()の戻り値だったので、今度はLogit.fit()を辿るたびに出るが、、こっちはもはや長旅の後の下り坂みたいなもんで、大したことはない。
LogitのfitはSuper(BinaryModel)のfitを呼び出しているので、そっちをみにいくと、そこで定義されていないから、さらにその上のDiscreteModelを見に行く、、、、ここでもfitは親(base.LikelihoodModel)のfitを呼び出していて、、、LikelihoodModel.fit()の中で、LikelihoodModelResults(self,xopt,Hinv
とnormalized_cov_paramsにHinvを渡している箇所を発見!!!
Hinvは同fitメソッドの中で計算されているH = self.hessian(xopt)
の逆行列で、hessian()はLogitの中で定義されているからこれで行き着いた!!
BinaryResults.summary()で表示されるロジスティック回帰の回帰係数の標準誤差は対数尤度関数のへシアンの逆行列の対角成分!!!
TODO
へシアンの逆行列を標準誤差に使っているのはパラメタの推定量の分散の下限はクラメル・ラオの定理でFisherの情報量行列の逆行列が下限だと知られているから。最尤推定量の時その分散は下限に一致することを使っている。
そのうちこれについてもまとめたい。