最尤法でフィッティング
pythonには便利なライブラリがいっぱいありますが，そのうちの一つにscipyがあります．
得られたデータ から,その確率密度分布を推定するときに,scipy.statsをよく使用します．
今回，自分が作成した関数を使用して最尤法でfittingをしたので備忘録としてメモります．
下記リンクを参考にしました．
https://austinrochford.com/posts/2015-03-03-mle-python-statsmodels.html
やり方としては，
GenericLikelihoodModelを
スーパークラスとして，新しくサブクラスを作成します．
例えば今回，任意の確率密度関数として
f(x;\beta,\alpha)=2\alpha^\beta x^{\beta-1} \frac{\exp({-\alpha x})}{\Gamma(\beta)}
というものを考えます．
単純にガンマ分布を2倍しただけです．
import numpy as np
from scipy import stats
import pandas as pd
from statsmodels.base.model import GenericLikelihoodModel
import math
import holoviews as hv
hv.extension('bokeh')
from bokeh.plotting import show
def newgamma_pdf(x, alpha, beta):
"""
新しく定義した関数です
"""
scale = 1/beta
return 2 * stats.gamma.pdf(x, alpha, loc=0, scale=scale)
class NewGammaFit(GenericLikelihoodModel):
"""
GenericLikelihoodModelを，
スーパークラスとして継承した
サブクラスを作成
"""
def __init__(self, endog, exog=None, **kwds):
"""
親クラスのコストラクタをオーバーライドします．
"""
if exog is None:
exog = np.zeros_like(endog)
super(GammaFitPMD, self).__init__(endog, exog, **kwds)
def nloglikeobs(self, params):
alpha = params[0]
beta = params[1]
return -np.log(newgamma_pdf(self.endog, alpha=alpha, beta=beta))
def fit(self, start_params=None, maxiter=10000, maxfun=5000, **kwds):
return super(GammaFitPMD, self).fit(start_params=start_params, maxiter=maxiter, maxfun=maxfun, **kwds)
これで準備OKです．
とりあえず普通のガンマ分布に従う乱数を発生させて，
それを今回作成した関数でfittingしてみます．
if __name__ == '__main__':
# とりあえず普通のガンマ分布に従う乱数を発生させる
data = stats.gamma.rvs(2, loc=0, scale = 1, size=1000)
#fittingさせる
model = NewGammaFit(data)
# 新しく作成した関数でfittingをする．fittingの際は初期値を設定する．
alpha, beta = 2, 1
results = model.fit(start_params=[alpha, beta])
print(results.summary)
X = np.linspace(0, 15, 1000)
hist = hv.Histogram(np.histogram(data, density=True))
scatter = hv.Scatter((X, pmd_pdf(X, 2, 1))).options(color="red")
plt = hv.render((scatter * hist).options(width=900, height=600))
show(plt)
これでできあがりです．
尚,グラフの出力は下記です．
scipy.curevefitとの違いは？
scipyには他にもcurve-fitという関数もあります．
こちらは最尤法ではなくて，最小二乗法でパラメータ推定をするようです．
最尤法と最小二乗法の違いは，
「実データとモデルとの誤差が正規分布するときは最小二乗法．そうじゃないときは最尤法」
と理解していますが，間違っていたら指摘して頂けると助かります．
パラメータ推定値の区間推定
ちなみに，
GenericLikelihoodModelクラスを継承したサブクラスでは，
色々と良い点があります．
上のコードで
print(result.summary)
とありますが，出力されるstd errが区間推定値です.
実はstd errの値が欲しくて，
GenericLikelihoodModelを継承したクラスを作成した，
という経緯もあります．
尚，実際の出力例が下です．
Current function value: 0.928056
Iterations: 45
Function evaluations: 85
GammaFitPMD Results
==============================================================================
Dep. Variable: y Log-Likelihood: -928.06
Model: GammaFitPMD AIC: 1856.
Method: Maximum Likelihood BIC: 1856.
Date: Wed, 08 Apr 2020
Time: 21:52:25
No. Observations: 1000
Df Residuals: 1000
Df Model: -1
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
const 5.2183 0.226 23.059 0.000 4.775 5.662
par0 1.7451 0.079 21.967 0.000 1.589 1.901
==============================================================================
最尤法でのパラメータ推定誤差
最尤法で求めた推定値のパラメーター誤差ですが，
N数が十分大きい時はクラメールラオの下限に近づくそうです．
さらにN数無限大の時は，誤差は正規分布に従うと考えちゃっても良さそうです．
(下記を参考にしています．)
http://www2.econ.osaka-u.ac.jp/~tanizaki/class/2015/micro_econome/03.pdf
実際にN数でどのように変わるのかシミュレーションしてみます．
ガンマ分布に従う乱数をN個発生させ，最尤法でfittingして，
パラメータα，βを求めます．
Nの数を100,5000,1000と変化させて，
その時のα，βの分布がどうなるかシミュレーションしました．
import numpy as np
from scipy import stats
from statsmodels.base.model import GenericLikelihoodModel
import holoviews as hv
hv.extension('bokeh')
from bokeh.plotting import show
def newgamma_pdf(x, alpha, beta):
scale = 1/beta
return stats.gamma.pdf(x, alpha, loc=10, scale=scale)
class GammaFit(GenericLikelihoodModel):
"""PMDの分布をgammafitするためのクラス"""
def __init__(self, endog, exog=None, **kwds):
if exog is None:
exog = np.zeros_like(endog)
super(GammaFit, self).__init__(endog, exog, **kwds)
def nloglikeobs(self, params):
alpha = params[0]
beta = params[1]
return -np.log(newgamma_pdf(self.endog, alpha=alpha, beta=beta))
def fit(self, start_params=None, maxiter=10000, maxfun=5000, **kwds):
return super(GammaFit, self).fit(start_params=start_params, maxiter=maxiter, maxfun=maxfun, **kwds)
if __name__ == '__main__':
size_list = [100, 500, 1000]
coeff_dic = {}
for size in size_list:
alpha_list = []
beta_list = []
for i in range(100):
data = stats.gamma.rvs(2, loc=10, scale=1, size=size)
print(data)
#data = data + stats.gamma.rvs(1, loc=0, scale=1, size=size)
model = GammaFit(data)
alpha, beta = 2, 1
results = model.fit(start_params=[alpha, beta])
alpha, beta = results.params[0], results.params[1]
X = np.linspace(0, 15, 1000)
alpha_list.append(alpha)
beta_list.append(beta)
coeff_dic[size] = [alpha_list, beta_list]
hist1 = hv.Histogram(np.histogram(coeff_dic[100][0]), label="N=100:{}".format(np.std(coeff_dic[100][0])), density=True)
hist2 = hv.Histogram(np.histogram(coeff_dic[500][0]), label="N=500:{}".format(np.std(coeff_dic[500][0])), density=True)
hist3 = hv.Histogram(np.histogram(coeff_dic[1000][0]), label="N=1000:{}".format(np.std(coeff_dic[1000][0])), density=True)
hist4 = hv.Histogram(np.histogram(coeff_dic[100][1]), label="N=100:{}".format(np.std(coeff_dic[100][1])), density=True)
hist5 = hv.Histogram(np.histogram(coeff_dic[500][1]), label="N=500:{}".format(np.std(coeff_dic[500][1])), density=True)
hist6 = hv.Histogram(np.histogram(coeff_dic[1000][1]), label="N=1000:{}".format(np.std(coeff_dic[1000][1])), density=True)
layout = hist1 + hist2 + hist3
layout2 = hist4 + hist5 + hist6
plt = hv.render(layout.options(width=900, height=600))
plt2 = hv.render(layout2.options(width=900, height=600))
show(plt)
show(plt2)
まずはαの結果です．
たしかにNが大きくなるにつれて，標準偏差は小さくなっています．
(グラフの上にある数字が標準偏差です)
N=100とN=1000を比べると,おおよそ√10程度違いますし，
クラメールラオの下限に従っていそうです．
続いてβの結果です．
βも同様にN数が大きくなるにつれて標準偏差は小さくなってます．
余談 ヒストグラムをfittingさせた時のN数の考え方は？
例えばbin数が10個のヒストグラムでfittingさせるとします．
そのヒストグラムを得るために何十個もデータがあったのに，
データ数としては10個になるのでしょうか？
unbinnedな方法を使うと，その限りではないそうです．
scipyのfittingもこの方法でやっているのかしら？
https://oku.edu.mie-u.ac.jp/~okumura/stat/141115.html