Up

ポアッソン分布を用いた多項式回帰モデル

 

ポアッソン分布を用いた多項式回帰モデルのベイズ分析は、このウェブサイトで扱っている。

ここでは、最尤法による分析例を示す。

まず、尤度比による最適な多項式の次数を求める。

次に、求めた次数の多項式によるポアッソン分布による回帰曲線を描く。

ファイルは、圧縮ファイルPoissonRegFiles.zipにまとめた。

 

 

 

多項式の次数の選択

 

 

警視庁のホームページ「各種交通事故発生状況（令和6年上半期）」からダウンロードしたPDFファイル「小学生の交通人身事故発生状況（令和6年上半期）更新日：2024年7月18日」に掲載されている月単位のデータを読み取り、CSVファイルDataMonth.csvとした（図A.1）。

図A.1　　月単位の事故件数（警視庁ホームページより）

 

図A.1のデータの棒グラフを描くと、図A.2のようになる。

図A.2　図A.1のデータの棒グラフ

 

月の事故件数がパラメータのポアッソン分布に従うとする。Spiegelhaltet (2019)は、殺人の件数にポアッソン分布を当てはめている。

独立変数の値１月から６月までに、１から６までの整数値を割り当て、パラメータが月の多項式で表されるとして、最尤法で多項式を求めると、図A.3のグラフを得た。

図A.3

 

まず、多項式の次数を求める。

月の事故件数がパラメータのポアッソン分布（岡本、2009、p.131）に従うとする。

さらに、がのの多項式で表されるとする。

このとき、尤度は、局所独立性の仮定の下で、次式で与えられる。

式（1）の値を最大にするパラメータ値を求める（最尤法）。

式（1）の最大化を極値探索法で行う場合、式（1）が確率の積であるので、アンダーフロー等のエラーの可能性がある。これを避けるために、対数をとって最大化を行う。

と置く。

ライブラリには、最大化ではなく最小化の関数が用意されているので、の最小化を行う。

最小化によりパラメータの値を求めるスクリプトをリストA.1のように作成した。

リストA.1では、極値探索が不適切な極小値に収束することを避けるために、から始め、次での解に基づいて、次での初期値を構成している。

スクリプトを実行すると、各の値に対するの値（で表す）が求められ、これらの値から

が求められる。この値は、帰無仮説の下で自由度１のカイ2乗分布に従う（岡本、2022、p.153）。

リストA.1の実行結果をまとめると、表A.1のようになる。

 

表A.1

ｐ値

1

22.51

 

 

2

21.78

1.452

0.228

3

20.03

3.516

0.061

4

19.74

0.566

0.452

5

18.99

1.504

0.220

 

多項式の次元数の変化に対するｐ値の変化のパターンから、３次の多項式が適切であると考えられる。

 

 

 

リストA.1　　多項式の次数を求めるためのスクリプト

 

"""

          ポアッソン分布による回帰モデル

          ポアッソン分布のパラメータを独立変数の多項式で表す

"""

import pandas as pd

import numpy as np

import scipy.stats as ss

import scipy.optimize as so

 

fin_nm = 'DataMonth.csv'

data_frm = pd.read_csv(fin_nm, encoding='shift-jis')

print(data_frm.keys().values)

data = data_frm.values[0]

print(data)

v_month = np.arange(1, 7, 1)    #    １月から６月までを、１から６までの整数で表す

print(v_month)

 

temp_b = np.array([data[0]])    #     極値探索の初期値設定に用いる

vLL = []                        #     極値探索で得た極値のリスト

for p in range(1,6):            #     p次多項式を用いる

    X= np.ones(6)

    for i in range(1,p+1):

        X = np.vstack((X, v_month**i))

    X = X.T

    print(X)

 

 

    def LL(b):                   #    対数尤度

        v = 0.0

        for i in range(6):

            mu = 0.0

            for j in range(p+1):

                mu += b[j]*X[i][j]

            if mu <= 1.0:           #    ポアッソン分布のパラメータmuを、１以上に保つ

                mu = 1.0

            v += ss.poisson.logpmf(data[i], mu)     #     対数尤度

        return v

 

 

    def min_func(b):      #    対数尤度にー１を掛ける。極値探索法が極小値を求めるものであるため

        return -LL(b)

 

    #

    #      初期値の設定。最初のp個の値は、前回の極小値を与えるものを継承

    #

    b0 = np.zeros(p+1)

    for k in range(p):

        b0[k] = temp_b[k]

    #       極小値を与えるパラメータを求める

    results = so.minimize(min_func, b0, method='Nelder-Mead',

                          options={'maxiter':10000})

    print(results)

    temp_b = results.x         #    極小値を与えるパラメータ値を

                               #    次の探索の初期値の設定に用いるために保存

 

    vLL.append(results.fun)    #    極小値をリストvLLに加える。

 

 

for i in range(5):

    print(f'p = {i+1}:  {vLL[i]:.7f}')

    if i > 0:

        ch2 = 2 * (vLL[i-1]-vLL[i])      #   カイ２乗の値、岡本（2022）、p.153

        cum_p = ss.chi2.cdf(ch2, 1)

        print(f'        chi2 = {2 * (vLL[i-1]-vLL[i]):.5f}' +

              f'   p = {1-cum_p:.5f}')

 

 

 

 

回帰曲線の描画

 

上のセクションにおいて、多項式の次数をとすることが適切であるとの結論が得られたので、３次の多項式を用いて分析を行う。スクリプトをリストB.1のように作成した。

次多項式による回帰モデルの係数を、から順番に求めているが、これは次多項式のパラメータ値を次多項式における極値探索の初期値設定に用いて、不適切な極小値に収束することを避けるためである。

求められた３次多項式による回帰曲線は、図B.1のようになった。

図B.1（図A.3再掲）

 

３次多項式で与えられるポアッソン分布のパラメータに対して、ポアッソン分布の2.5パーセンタイル値と97.5パーセンタイル値を求め、95％区間とした。図中、橙色の縦の線分で表されている。データ値が全てこれらの95％区間内にあるので、モデルはデータに適合していると言える。事故は２月が最も少なく、５月に最大となる傾向が３次の多項式で表されている。

 

 

 

リストB.1　　３次多項式を用いたポアッソン分布による分析

 

import pandas as pd

import numpy as np

import scipy.stats as ss

import scipy.optimize as so

import matplotlib.pyplot as plt

import japanize_matplotlib

 

fin_nm = 'DataMonth.csv'

data_frm = pd.read_csv(fin_nm, encoding='shift-jis')

data = data_frm.values[0]

v_month = np.arange(1, 7, 1)

 

temp_b = np.array([data[0]])              #     極値探索の初期値

vLL = []

for p in range(1,4):

    #     独立変数の行列Xの設定：　岡本、2019、p.91

    X= np.ones(6)

    for i in range(1,p+1):

        X = np.vstack((X, v_month**i))

    X = X.T

 

    def LL(b):          #   対数尤度関数

        v = 0.0

        for i in range(6):

            mu = 0.0

            for j in range(p+1):

                mu += b[j]*X[i][j]

            if mu <= 1.0:

                mu = 1.0

            v += ss.poisson.logpmf(data[i], mu)

        return v

 

 

    def min_func(b):    #    対数尤度関数の符号を反転したもの

        return -LL(b)

 

 

    b0 = np.zeros(p+1)      #    パラメータの初期値

    for k in range(p):

        b0[k] = temp_b[k]   #    前回のパラメータ値を用いる

 

    results = so.minimize(min_func, b0, method='Nelder-Mead',

                          options={'maxiter':10000})

   

    temp_b = results.x        #    次回の初期値設定に用いるために保存

    vLL.append(results.fun)   #    求めた関数の最小値をリストvLLに格納

 

 

beta = temp_b       #   p=3におけるパラメータ値

 

X= np.ones(6)

p = 3

for i in range(1,p+1):

    X = np.vstack((X, v_month**i))

X = X.T

y_hat = X @ beta    #   回帰多項式の値、mu

 

p025 = []           #    2.5パーセンタイル値

p975 = []           #   97.5パーセンタイル値

for mu in y_hat:

    p025.append(ss.poisson.ppf(0.025, mu))

    p975.append(ss.poisson.ppf(0.975, mu))

   

               

plt.ylabel('件数', fontsize=16)

plt.title('小学生の交通人身事故', fontsize=16)

x_coord = data_frm.keys().values                     #    １月から６月までの値（文字列）

plt.bar(x_coord, data, alpha=0.5, label='Data')      #    データ値のプロット

#     多項式によるmuの予測値の折れ線グラフ

plt.plot(x_coord, y_hat, c='r', ls = '-', lw=3, label='Regression Line')

#     95%区間の描画

for i, cnf in enumerate(zip(p025, p975)):

    plt.plot([x_coord[i], x_coord[i]], [cnf[0], cnf[1]],

             'o', c='r', linestyle='-', lw=2)

plt.plot([], 'o', c='r', linestyle='-', lw=2, label='95% Interval')

plt.fill_between(x_coord, p025, p975, color='y', alpha=0.3)

plt.legend()

plt.tight_layout()

plt.show()

 

 

 

 

参考文献

岡本安晴（2009）統計学を学ぶための数学入門［下］、培風館

岡本安晴（2022）データを使いこなすための統計入門、22世紀アート

Spiegelhalter, D. (2019) The Art of Statistics. Penguin Random House.

 

 

Home