scipyを使って任意の関数で最尤法フィッティング

#最尤法でフィッティング
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