一般化線形モデル(GLM)の理論と実装③【ポアソン回帰編】|スクラッチ実装で理解する機械学習アルゴリズム #4

f:id:lib-arts:20190417225222p:plain

連載の経緯の詳細は#1でまとめましたが、本シリーズではあえてスクラッチ実装を元に機械学習アルゴリズムを実装していくことで、アルゴリズムの概要を掴んだり理論の流れを掴んだりできるようにできればと思います。

実装のほとんどが車輪の再発明に近くなりますが、ベーシックなアルゴリズムをあえて一から自分で追ってみるというのは引き出しを増やすという意味で非常に有意義です。
#2、#3では一般化線形モデル(GLM; Generalized Linear Model)の具体例の一つとしてロジスティック回帰について実装を通して理解を深める手助けとしました。
https://lib-arts.hatenablog.com/entry/ml_scratch2
https://lib-arts.hatenablog.com/entry/ml_scratch3
#4ではGLMの他の例としてポアソン回帰を取り扱います。
以下目次になります。
 
1. GLMの復習とポアソン回帰の概要
2. ポアソン回帰の例題と解き方の紹介
3. 実装で理解するポアソン回帰
4. まとめ

 

1. GLMの復習とポアソン回帰の概要
GLMは線形モデル(LM; Linear Model)の拡張で、具体的には二点における拡張です。
====
1. 目的変数yの分布
-> LMでは正規分布を仮定していたのに対し、GLMでは分布を指数型分布族に拡張します。

2. 目的変数の期待値と線形式部分の関係性
-> LMでは恒等関数を仮定していたのに対し、GLMではlogやexpなどの組み合わせでできるシンプルな関数(リンク関数)を用います。
====
#2では具体例としてロジスティック回帰を取り扱いましたが、今回はポアソン回帰を取り扱います。簡易化のために一次元で議論を進めますが、多次元でも同様の議論を行うことができます。
ポワソン回帰に用いられるポアソン分布は、計数(count)データや度数(frequency)データを取り扱うために利用されることが多いです。確率分布の数式としては、イベントの発生数をy、平均発生数\muをパラメータとして以下の数式で表せます。
f(y|\mu)=\frac{\mu^{y}exp(-\mu)}{y!} y=0,1,2..
ここで注意なのが、yが0以上の整数であることと、パラメータ\muの設定にあたって何らかの割合または率を利用することが多いことです。例えば商品の売り上げにおける店を訪れる数と購入する数の割合を利用したり、自動車事故件数の分析にあたって人口1,000人当たりの事故件数の割合を利用したりです。他にも単位時間あたりの発生率を利用することで、台風の発生数をモデリングすることなども可能になります。
またポアソン回帰ではg(\mu)=log(\mu)=β_{0}+β_{1}xというリンク関数を用います。この数式を変換することで、\mu=exp(β_{0}+β_{1}x)を得ることができ、\muを0よりも大きな値に設定することができます。(このようにリンク関数は問題の設定に応じて設定するものと考えておけば良いかと思います。)
ポアソン回帰の概要に関しては一通りまとめられたと思うので、2節では実際に例題とその解き方をご紹介します。

 

2. ポアソン回帰の例題と解き方の紹介
1節では基本的なポアソン回帰の概要についてまとめました。とはいえ話が抽象的でわかりにくいと思いますので、2節では実際の例題を元にポアソン回帰について取り扱います。下記のサンプルを元に回帰を行ってみます。

import numpy as np
import matplotlib.pyplot as plt

x = np.arange(1,8,0.1)
y = np.exp2(x) + np.random.randn(70) + 4

plt.scatter(x,y)
plt.show()

上記を実行すると下記のような結果が得られます。

f:id:lib-arts:20190519165752p:plain
上記に対して具体的なイメージをするなら、学習効果をイメージしていただければと思います。学習にあたっての単位時間が1増えるごとに、学習効果は指数関数状に増加すると巷に言われています。それらを数値化したとしましょう。(例題なので若干恣意的な作成となっています。)
これに対し通常の回帰分析を行うと不十分な結果になりそうだということが散布図を見ただけでもわかります。
問題の解法としてはポアソン回帰なので、リンク関数はg(μ)=log(μ)、目的変数の確率分布はポアソン分布に従うとします。これらの仮定の元での最尤法により、尤度関数L(μ)は下記のように計算できます。
L(μ)=\Pi_{k=1}^n\frac{μ_{k}^{y_{k}}exp(-μ_{k})}{y_{k}!}
このままだと掛け算が多くて計算しづらいのでL(μ)の対数をとります。
log(L(μ))=\Sigma_{k=1}^n(y_{k}log(μ_{k})-μ_{k}-logy_{k}!)
ここで、パラメータβの関数となっているのがμ_{k}のみのため、logy_{k}!の項についてはβについて偏微分するにあたっては定数(const.)として無視することができます。対数尤度をlとして、β_{0}β_{1}についてそれぞれ偏微分を計算します。
\frac{δl}{δβ_{0}}=\Sigma_{k=1}^n(y_{k}-exp(β_{0}+β_{1}x_{k}))
\frac{δl}{δβ_{1}}=\Sigma_{k=1}^n(y_{k}-exp(β_{0}+β_{1}x_{k}))x_{k}
上記の数式で計算した勾配をベースに勾配降下法(Gradient Descent)を用いて3節では学習(パラメータの推定)を行います。


3. 実装で理解するポアソン回帰
2節ではポアソン回帰に関しての例題の紹介と解法の提示を行いました。3節ではこれに対応した実装を取り扱っていければと思います。

beta = np.array([0.,0.])
for i in range(500000):
    exp_lindif_0 = -y + np.exp(beta[0]+beta[1]*x)
    exp_lindif_1 = x*(-y+np.exp(beta[0]+beta[1]*x))
    grad = np.array([np.sum(exp_lindif_0),np.sum(exp_lindif_1)])
    beta[0] = beta[0] - 0.00001*grad[0]
    beta[1] = beta[1] - 0.00001*grad[1]
    if (i+1)%100000==0:
        print("STEP_"+str(i+1))
        print(beta)
        print(grad)
        print("======")
print(beta)

上記を実行することによってβの値を求めることができます。実行結果は下記のようになります。

f:id:lib-arts:20190519171416p:plain
また、下記を実行するとグラフを描くことができます。

x_ = np.arange(1,9,0.05)
y_ = np.exp(beta[0]+beta[1]*x_)
plt.scatter(x_,y_,color="red")
plt.scatter(x,y)
plt.show()

f:id:lib-arts:20190519171449p:plain
上記が実行結果です。推定結果を元に予測した赤の点が元の青い点に重なっているように見えます。
また、推定結果を元の式と見比べると、np.exp(0.61036691)=1.8411067953634672で2に近いことから、概ね大元の数式を推定できているということがわかります。


4. まとめ
#4ではGLMの例としてポアソン回帰について取り扱いました。
#5では勾配法以外の最適化も取り扱ってみれればということで、MCMC(Markov Chain Monte Carlo)法による最適化について取り扱っていければと思います。