Introduction

　ここでは，モデル分析を行う際にしばしば必要となる収束計算のR上での実装方法について，ニュートン・ラフソン（Newton-Raphson）法（以下，単に「ニュートン法」もしくは「NR法」と呼ぶ）を題材にして見ていきます．また，最適化計算の際によく用いられるoptim関数についても，その使い方を見ていきたいと思います．このWEBページの元となっているRmdファイルについては，筆者のWEBサイトに設置したこちらのRmdが入ったzipファイルをダウンロードしてください．
　このページの主眼はあくまで，収束計算のフルスクラッチでの実装方法を紹介することにあり，各種最適化アルゴリズムや，収束計算を要する各種統計手法の詳細には言及しません．加えて，筆者自身は数学は（算数の頃から）得意ではありませんので，以下で行う最適化理論・統計理論に関する解説には，数多くの誤りを含むものと想定されます．もし，誤りの箇所がありましたら，メールなどでお知らせくださると助かります．
　ページの作成にあたって，以下の文献を参考にしました．特に，川田(2008)は図を用いた最適化アルゴリズムの解説が分かりやすく，まずは2変数バージョンから各種手法を解説してくれるという点で，数学アレルギー持ちでも安心して読めます．また，蓑谷(2013)は，以下で紹介するポアソン回帰を含む一般化線形モデルに関する日本語の教科書として，最も網羅的かつ分かりやすいもののひとつだと思います．

　このページの数式の記述は，基本的にLaTeXベースでなされていますが，Rmarkdownの仕様上，ギリシャ文字を太字表記することができません．予めご了承ください．

Updates

数式・コードのタイポを修正しました（2020年6月7日）．
数式・コードのタイポを再度修正しました（2020年10月29日）．

NR Method

Basic Concept

　ニュートン法は，方程式系を数値計算によって解くための，反復計算に基づく求根アルゴリズム（$f(x)=0$を満たす根$x$を見つけるためのアルゴリズム）の1つです．以下では，その直観的な解釈について述べていきます．

図1：ニュートン法のイメージ

最初に， $f(x)$上のある点$A(\alpha,f(\alpha))$における$y=f(x)$の接線は次のように与えられます．

\[ y-f(\alpha) = f'(\alpha)(x - \alpha) \Rightarrow y=f(\alpha)+f'(\alpha)(x - \alpha) \tag{1} \]

ゆえに図1上での，$P_{1}$における$y=f(x)$の接線は次のように与えられます．

\[ y=f(x_{1})+f'(x_{1})(x-x_{1}) \tag{2} \]

この接線(2)と $x$軸との交点$(0,x_{2})$は次の式で求められます．

\[ 0=f(x_{1})+f'(x_{1})(x_{2}-x_{1}) \Rightarrow x_{2}=x_{1}-\frac{f(x_{1})}{f'(x_{1})} \tag{3} \]

また，$P_{2}$における$y=f(x)$の接線と$x$軸との交点$(0,x_{3})$も同様にして，

\[ 0=f(x_{2})+f'(x_{2})(x_{3}-x_{2}) \Rightarrow x_{3}=x_{2}-\frac{f(x_{2})}{f'(x_{2})} \tag{4} \]

と求められます．このような操作を繰り返すことにより得られる$x_{k}$は$f(x)=0$の解に漸近すると考えられます．これがニュートン法のざっくりとした考え方です．
　次に，テイラー展開によるニュートン法の解釈についても説明しておきます．$x=α$近傍における非線形関数$f(x)$のテイラー展開は以下のように与えられます．

\[ f(x)=f(\alpha)+\frac{1}{1!}f'(\alpha)(x-\alpha)+\frac{1}{2!}f''(\alpha)(x-\alpha)^2+\ldots \tag{5} \]

よって，$x=x_{k}$近傍における非線形関数$f(x)$のテイラー展開は，

\[ f(x)=f(x_{k})+\frac{1}{1!}f'(x_{k})(x-x_{k})+\frac{1}{2!}f''(x_{k})(x-x_{k})^2+\ldots \tag{6} \]

となり，1次の項まででの近似を行おうと考えれば，

\[ f(x)≈f(x_{k})+f'(x_{k})(x-x_{k}) \tag{7} \]

と表すことができます．したがってこれは，$y=f(x)$上の点$(x_{k},f(x_{k}))$における接線と同じであるということが分かります．よって，先ほど図1を用いて示したプロセスとのアナロジーを考えれば，以下に示す漸化式

\[ f(x_{k})+f'(x_{k})(x_{k+1}-x_{k})=0 \Rightarrow x_{k+1}=x_{k}-\frac{f(x_{k})}{f'(x_{k})} \tag{8} \]

を繰り返し計算することによって，$f(x)=0$の近似解を数値的に求めることが可能です．

Example on R: Cubic Equation

　ここでは，以下に示す簡単な1次元の非線形関数の求根を例に，実装方法を見てみます．

\[ f(x)=x^3+2x^2+3x+4 \tag{9} \]

まず，アルゴリズムを動かすにあたり，以下のように条件を設定します．
xは初期値，maxitは繰り返しの最大回数，tolは収束判定を行う為の閾値を表します．

x <- 10
maxit <- 50
tol <- 1e-08

以下では，ニュートン法の本体部分を示します．
1行目では，ループ変数（何回目の繰り返しかを示す変数）kを1からmaxit（今回は50）の範囲で動かす，すなわち，最大maxit回反復計算を行うことを宣言しています．
2行目では，式(9)で示した$f(x)$に対し，k番目のxを代入し計算した結果を，fというオブジェクトに代入しています．1番目のxは，上で設定した初期値10です．
3行目では，$f'(x)$に対しても同様に，k番目のxを代入し計算した結果を，gというオブジェクトに代入しています．
4行目では，反復計算によってxの値がどのように変化したかを計算終了後に表示する為に，xについてprint関数を記述しています．
5行目では，k番目のxをold.xというオブジェクトに格納しています．
6行目がいちばん重要な部分で，式(8)で示したニュートン法の漸化式にk番目のxを代入し，(k+1)番目のxを計算しています．
7行目が収束判定を行う部分で，k番目のxと(k+1)番目のxの差の絶対値が閾値tol未満であれば，解が収束したと判断して反復計算を抜けます．
8行目で，繰り返し回数kを表示します．

for(k in 1:maxit){
    f <- x^3 + 2*x^2 + 3*x +4
    g <- 3*x^2 + 4*x +3
    print(x)
    old.x <- x
    x <- old.x - f/g
    if(abs(old.x-x)<tol)break
}

## [1] 10
## [1] 6.402332
## [1] 3.977078
## [1] 2.312341
## [1] 1.110693
## [1] 0.1083721
## [1] -1.145657
## [1] -1.860904
## [1] -1.67571
## [1] -1.651029
## [1] -1.650629
## [1] -1.650629

## [1] 12

今回は，解の収束過程を見るために，敢えて初期値を最終的に求められる解から遠い場所におきました．式(9)の関数に実際に求められた解（赤丸）を重ねてプロットしてみます．

plot(function(x){x^3 + 2*x^2 + 3*x +4},xlim=c(-3,3),ylim=c(-15,30))
abline(h=0,v=0)
points(x,0,col="red")

解はうまいこと収束し，$f(x)=0$を満たす$x$を見つけることができたと判断してよさそうです．

NR Method for Multiple Parameters

　ここでは，複数パラメータの場合のニュートン法についても見ていきたいと思います．まず，点$\mathbf{x}=\mathbf{a}$近傍における$p$変数関数のテイラー展開を2次の項までで行うと，

\[ f(\mathbf{x})≈f(\mathbf{a})+\frac{1}{1!}\sum^p_{i=1}\frac{\partial f(\mathbf{a})}{\partial x_{i}}(x_{i}-a_i) + \frac{1}{2!}\sum^p_{i=1}\sum^p_{j=1}\frac{\partial^2 f(\mathbf{a})}{\partial x_{i}\partial x_{j}}(x_{i}-a_i)(x_{j}-a_j) \tag{10} \]

という風になります．ただし，微分については，

\[ \frac{\partial f(\mathbf{a})}{\partial x_{i}}= \frac{\partial f(\mathbf{x})}{\partial x_{i}}|_{\mathbf{x}=\mathbf{a}} \tag{11} \]

という意味の表記を導入しています．式(10)はこのままでは色々と扱いづらいので，行列形式に書き直します．

\[ f(\mathbf{x}) ≈ f(\mathbf{a}) + \left[ \begin{array}{ccc} \frac{\partial f(\mathbf{a})}{\partial x_{1}} & \ldots & \frac{\partial f(\mathbf{a})}{\partial x_{p}} \end{array} \right] \left[ \begin{array}{c} (x_1-a_1) \\ \ldots \\ (x_p-a_p) \end{array} \right] + \frac{1}{2} \left[ \begin{array}{ccc} (x_1-a_1) & \ldots & (x_p-a_p) \end{array} \right] \left[ \begin{array}{ccc} \frac{\partial^2 f(\mathbf{a})}{\partial x_{1} \partial x_{1}} & \ldots & \frac{\partial^2 f(\mathbf{a})}{\partial x_{1} \partial x_{p}} \\ \vdots & \ddots & \vdots \\ \frac{\partial^2 f(\mathbf{a})}{\partial x_{p} \partial x_{1}} & \ldots & \frac{\partial^2 f(\mathbf{a})}{\partial x_{p} \partial x_{p}} \end{array} \right] \left[ \begin{array}{c} (x_1-a_1) \\ \ldots \\ (x_p-a_p) \end{array} \right]\\ =f(\mathbf{a}) + \nabla f(\mathbf{a}) ' (\mathbf{x} - \mathbf{a}) + \frac{1}{2} (\mathbf{x} - \mathbf{a}) ' \nabla^2 f(\mathbf{a}) (\mathbf{x} - \mathbf{a}) \tag{12} \]

よって，点$\mathbf{x}=\mathbf{x}_{k}$近傍における2次の項までのテイラー展開も次のように表されます．

\[ f(\mathbf{x})≈f(\mathbf{x}_{k}) +\nabla f(\mathbf{x}_{k})'(\mathbf{x}-\mathbf{x}_{k}) +\frac{1}{2} (\mathbf{x} - \mathbf{x}_{k})'\nabla^2 f(\mathbf{x}_{k})(\mathbf{x}-\mathbf{x}_{k}) \tag{13} \]

1次の項までで近似を行おうとすれば，式(13)は更に以下のように近似されます．

\[ f(\mathbf{x})≈f(\mathbf{x}_{k})+\nabla f(\mathbf{x}_{k})'(\mathbf{x}-\mathbf{x}_{k}) \tag{14} \]

　ここで，以下のような連立方程式をニュートン法で解くことを考えましょう．

\[ \left \{ \begin{array}{c} f_{1}(\mathbf{x})=0 \\ \vdots \\ f_{p}(\mathbf{x})=0 \end{array} \right. \tag{15} \]

1変数関数の場合とのアナロジーを考えれば，$p$変数のニュートン法の場合でも同様，ある点 $\mathbf{x}=\mathbf{x}_k$近傍における$f_{1}(x),\ldots,f_{p}(x)$の1次近似関数

をゼロとした式，即ち以下に示す連立方程式

\[ \left \{ \begin{array}{c} f_{1}(\mathbf{x}_{k})+\nabla f_{1}(\mathbf{x}_{k})'(\mathbf{x}-\mathbf{x}_{k}) \\ \vdots \\ f_{p}(\mathbf{x}_{k})+\nabla f_{p}(\mathbf{x}_{k})'(\mathbf{x}-\mathbf{x}_{k}) \end{array} \right. \Rightarrow \left[ \begin{array}{c} f_{1}(\mathbf{x}_{k}) \\ \ldots \\ f_{p}(\mathbf{x}_{k}) \end{array} \right] + \left[ \begin{array}{ccc} \frac{\partial f_{1}(\mathbf{x}_{k})}{\partial x_{1}} & \ldots & \frac{\partial f_{1}(\mathbf{x}_{k})}{\partial x_{p}} \\ \vdots & \ddots & \vdots \\ \frac{\partial f_{p}(\mathbf{x}_{k})}{\partial x_{1}} & \ldots & \frac{\partial f_{p}(\mathbf{x}_{k})}{\partial x_{p}} \end{array} \right] (\mathbf{x}-\mathbf{x}_{k}) =\mathbf{0} \tag{17} \]

の解$\mathbf{x}=\mathbf{x}_{k+1}$を求める，という操作を反復的に行うことによって，連立方程式(15)の近似解を求めることが可能です．ここで，

\[ \mathbf{f}(\mathbf{x}_{k}) = \left[ \begin{array}{c} f_{1}(\mathbf{x}_{k}) \\ \ldots \\ f_{p}(\mathbf{x}_{k}) \end{array} \right] , \mathbf{J}(\mathbf{x}_{k}) = \frac{\partial \mathbf{f}(\mathbf{x}_{k})}{\partial \mathbf{x}_{k}} = \left[ \begin{array}{ccc} \frac{\partial f_{1}(\mathbf{x}_{k})}{\partial x_{1}} & \ldots & \frac{\partial f_{1}(\mathbf{x}_{k})}{\partial x_{p}} \\ \vdots & \ddots & \vdots \\ \frac{\partial f_{p}(\mathbf{x}_{k})}{\partial x_{1}} & \ldots & \frac{\partial f_{p}(\mathbf{x}_{k})}{\partial x_{p}} \end{array} \right] \tag{18} \]

と置き換えれば，$\mathbf{x}=\mathbf{x}_{k+1}$として，式(17)を次のように書き換えることが可能です．

\[ \mathbf{f}(\mathbf{x}_{k})+\mathbf{J}(\mathbf{x}_{k})(\mathbf{x}-\mathbf{x}_{k})=\mathbf{0} \Rightarrow \mathbf{x}_{k+1}=\mathbf{x}_{k}-[\mathbf{J}(\mathbf{x}_{k})]^{-1}\mathbf{f}(\mathbf{x}_{k}) \tag{19} \]

式(8)と比較すると$f(x)$の1階微分$f'(x)$の逆数が，$\mathbf{f}(\mathbf{x})$のヤコビアン$\mathbf{J}(\mathbf{x})$の逆行列に置き換わったものの，式の形は1変数の場合も$p$変数の場合も非常に似ていることが分かると思います．

Example on R: Count Data Regression (Theory)

　カウントデータとは，一般にある事象が決まった時間内に起こった回数を数え上げることで集めた統計（非負の整数）を指し，その発生頻度を調べ，分布関数を特定化し，それに基づいて回帰分析することをカウントデータ分析と呼びます．カウントデータの具体的な例としては，病院に行く回数や地震の発生件数，台風の上陸件数などが挙げられます．カウントデータの特徴としては，事象は稀にしか起こらず，多くの期間では事象が起こらない，いわゆるゼロ事象だということであるということが挙げられます．
　統計学上，稀にしか起こらない事象の発生確率はポアソン分布で表すことが多いので，こうした事象に関して回帰分析を行う際にも，事象発生のデータ生成過程にポアソン分布を仮定した「ポアソン回帰モデル」と呼ばれるモデルを用いた分析を行います．
まず，ポアソン分布は次のように定義されます．

\[ f(Y_{i}=y_{i})=\frac{{\rm exp}(-\mu_{i})\mu_{i}^{y_{i}}}{y_{i}!},\mu_{i}>0,y=0,1,2, \ldots \\ E[Y_{i}]=Var[Y_{i}]=\mu_{i} \tag{20} \]

ポアソン分布はこのように，未知のパラメータ$\mu_{i}$が決まれば分布形が決まる極めて簡便な分布です．このポアソン分布を用いたポアソン回帰モデルは次のように定義できます．

\[ f(Y_{i}=y_{i})=\frac{{\rm exp}(-\mu_{i})\mu_{i}^{y_{i}}}{y_{i}!},\mu_{i}={\rm exp}(\mathbf{x}_{i}'\mathbf{\beta}) \tag{21} \]

ただし，

\[ \mathbf{X} = \left[ \begin{array}{ccc} x_{11} & \ldots & x_{1p} \\ \vdots & \ddots & \vdots \\ x_{n1} & \ldots & x_{np} \end{array} \right] , \mathbf{x}_{i} = \left[ \begin{array}{c} x_{i1} \\ \vdots \\ x_{ip} \end{array} \right] , \mathbf{\beta} = \left[ \begin{array}{c} \beta_{1} \\ \vdots \\ \beta_{p} \end{array} \right], \mathbf{x}_{i}'\mathbf{\beta}=x_{i1}\beta_{1}+\ldots+x_{ip}\beta_{p} \tag{22} \]

であるとします．
ここでは，式(21)を以下で表す対数尤度関数に変換し，最尤推定する方法を考えます．

\[ {\rm ln}L(\mathbf{\beta})=\sum^n_{i=1}[y_{i}\mathbf{x}_{i}'\mathbf{\beta}-{\rm exp}(\mathbf{x}_{i}'\mathbf{\beta})-{\rm ln}y_{i}!] \tag{23} \]

式(23)を$\mathbf{\beta}$について1階微分した結果は次のように与えられます．

\[ \frac{\partial {\rm ln}L(\mathbf{\beta})}{\partial \mathbf{ \beta}}= \sum^n_{i=1}[y_{i}-{\rm exp}(\mathbf{x}_{i}'\mathbf{\beta})]\mathbf{x}_{i}= \mathbf{X}'[\mathbf{y}-{\rm exp}(\mathbf{X}'\mathbf{\beta})] \tag{24} \]

式(24)をゼロとした式，すなわち

\[ \frac{\partial {\rm ln}L(\mathbf{\beta})}{\partial \mathbf{\beta}}=\mathbf{0} \tag{25} \]

は，対数尤度関数$\rm{ln}L(\beta)$が$\beta$について最大化されるときの1階条件を表しています．ポアソン回帰に限らず，対数尤度関数の回帰係数についての1階微分は，一般的にグラディエント（gradient）と呼ばれることが多いです．式(25)を書き下せば

\[ \mathbf{g}(\mathbf{\beta})= \left[ \begin{array}{c} \frac{\partial {\rm ln}L(\mathbf{\beta})}{\partial \beta_{1}} \\ \vdots \\ \frac{\partial {\rm ln}L(\mathbf{\beta})}{\partial \beta_{p}} \end{array} \right] = \left[ \begin{array}{c} [y_{1}-{\rm exp}(\mathbf{x}_{1}'\mathbf{\beta})]x_{11} + \ldots + [y_{n}-{\rm exp}(\mathbf{x}_{n}'\mathbf{\beta})]x_{1n}\\ \vdots \\ [y_{1}-{\rm exp}(\mathbf{x}_{1}'\mathbf{\beta})]x_{p1} + \ldots + [y_{n}-{\rm exp}(\mathbf{x}_{n}'\mathbf{\beta})]x_{pn} \end{array} \right] = \left[ \begin{array}{c} 0 \\ \vdots \\ 0 \end{array} \right] \tag{26} \]

と表せますから，結局のところ，以下に示すような$p$変数の非線形連立方程式を解く必要があることが分かります．式(15)と比較を行ってみてください．

\[ \left \{ \begin{array}{c} \frac{\partial {\rm ln}L(\mathbf{\beta})}{\partial \beta_{1}}=0 \\ \vdots \\ \frac{\partial {\rm ln}L(\mathbf{\beta})}{\partial \beta_{1}}=0 \end{array} \right. \tag{27} \]

以下では，ニュートン法のアルゴリズムを用いて最尤推定を行っていきます．$\mathbf{\beta}=\mathbf{\beta}_{k}$とおくと，式(19)における$\mathbf{f}(\mathbf{x}_{k})$に対応するものが，今回の分析では，

\[ \mathbf{g}(\mathbf{\beta}_{k})= \left[ \begin{array}{c} \frac{\partial {\rm ln}L(\mathbf{\beta}_{k})}{\partial \beta_{1(k)}} \\ \vdots \\ \frac{\partial {\rm ln}L(\mathbf{\beta}_{k})}{\partial \beta_{p(k)}} \end{array} \right] \tag{28} \]

となっていることを確認してください．$(k)$は「k回目の計算で得られた」という意味です．次に，ニュートン法を実行する為に必要なヤコビアン$\mathbf{J}(\mathbf{x})$を求めていくわけですが，今回の分析において$\mathbf{J}(\mathbf{x})$に対応するものは，

\[ \mathbf{H}(\mathbf{\beta})= \frac{\partial \mathbf{g}(\mathbf{\beta})}{\partial \mathbf{\beta}}= \left[ \begin{array}{ccc} \frac{\partial^2 {\rm ln}L(\mathbf{\beta})}{\partial \beta_{1} \partial \beta_{1}} & \ldots & \frac{\partial^2 {\rm ln}L(\mathbf{\beta})}{\partial \beta_{1} \partial \beta_{p}} \\ \vdots & \ddots & \vdots \\ \frac{\partial^2 {\rm ln}L(\mathbf{\beta})}{\partial \beta_{p} \partial \beta_{1}} & \ldots & \frac{\partial^2 {\rm ln}L(\mathbf{\beta})}{\partial \beta_{p} \partial \beta_{p}} \end{array} \right] \tag{29} \]

になることが分かります．ポアソン回帰に限らず，対数尤度関数の回帰係数についての2階微分は，一般的にヘシアン（Hessian）と呼ばれることが多いです．同様に，$\mathbf{\beta}=\mathbf{\beta}_{k}$とおくと，式(19)における$\mathbf{J}(\mathbf{x}_{k})$に対応するものが，今回の分析では，

\[ \mathbf{H}(\mathbf{\beta}_{k})= \left[ \begin{array}{ccc} \frac{\partial^2 {\rm ln}L(\mathbf{\beta}_{k})}{\partial \beta_{1(k)} \partial \beta_{1(k)}} & \ldots & \frac{\partial^2 {\rm ln}L(\mathbf{\beta}_{k})}{\partial \beta_{1(k)} \partial \beta_{p(k)}} \\ \vdots & \ddots & \vdots \\ \frac{\partial^2 {\rm ln}L(\mathbf{\beta}_{k})}{\partial \beta_{p(k)} \partial \beta_{1(k)}} & \ldots & \frac{\partial^2 {\rm ln}L(\mathbf{\beta}_{k})}{\partial \beta_{p(k)} \partial \beta_{p(k)}} \end{array} \right] \tag{30} \]

になっていることに注意してください（「既に微分された」対数尤度のヘシアンなので，1階微分ではなく2階微分で表される形になっています）．よって，式(19)との対応関係から，ニュートン法での反復計算による回帰パラメータの推定は，次のように定式化されます．

\[ \mathbf{\beta}_{k+1}=\mathbf{\beta}_{k}-[\mathbf{H}(\mathbf{\beta}_{k})]^{-1}\mathbf{g}(\mathbf{\beta}_{k}) \tag{31} \]

　ただし，多変数関数に関し，ニュートン法を用いて解を求めることは，次の2つの理由から非現実的です．

$\mathbf{H}(\mathbf{\beta}_{k})$が正値定符号でない場合，$\mathbf{H}(\mathbf{\beta}_{k})$を繰り返し更新しながら最適化を行うと，鞍点（ある変数方向で見れば最大値だが，別の変数方向で見ると最小値を取るような点）がもたらされてしまう可能性があります．行列の符号と多変数関数の極大・極小に関する解説は，例えばこちらのサイトを参照してください．
煩雑な2階微分を解析的・数値的に計算する必要があります．

これらの問題を補正する為の方法が，スコア法（scoremethod）と呼ばれるものです．スコア法では，$\mathbf{H}(\mathbf{\beta}_{k})$を用いる代わりに，その期待値（符号は逆だが）であるフィッシャー情報行列$\mathbf{I}(\mathbf{\beta}_{k})$を用いた反復計算を行います．

\[ \mathbf{\beta}_{k+1}=\mathbf{\beta}_{k}+[\mathbf{I}(\mathbf{\beta}_{k})]^{-1}\mathbf{g}(\mathbf{\beta}_{k}) \tag{32} \]

ポアソン回帰の例では，情報行列$\mathbf{I}(\mathbf{\beta}_{k})$は次のように与えられます．

\[ \mathbf{I}(\mathbf{\beta})= - E[\mathbf{H}(\mathbf{\beta})]= -E\left[\frac{\partial^2 {\rm ln}L(\mathbf{\beta})}{\partial \beta \partial \beta'}\right] = \sum^n_{i=1}[{\rm exp}(\mathbf{x}_{i}'\mathbf{\beta})]\mathbf{x}_{i}\mathbf{x}_{i}'= \\ \left[ \begin{array}{ccc} x_{11} & \ldots & x_{n1} \\ \vdots & \ddots & \vdots \\ x_{1p} & \ldots & x_{np} \end{array} \right] \left[ \begin{array}{ccc} {\rm exp}(\mathbf{x}_{1}'\mathbf{\beta}) & \ldots & 0 \\ \vdots & \ddots & \vdots \\ 0 & \ldots & {\rm exp}(\mathbf{x}_{n}'\mathbf{\beta}) \end{array} \right] \left[ \begin{array}{ccc} x_{11} & \ldots & x_{1p} \\ \vdots & \ddots & \vdots \\ x_{n1} & \ldots & x_{np} \end{array} \right] \\ =\mathbf{X}'\mathbf{W}\mathbf{X} \tag{33} \]

Example on R: Count Data Regression (Empirics)

　ここでは例として，若年男性における，過去の懲役歴・雇用状態と1年あたり逮捕回数との関係を，ポアソン回帰を用いて調べてみます．分析に使うデータセットの原典は，Grogger (EI 1991)です．なお，この分析例は，Wooldridge (2020)の17章で示されている分析結果をレプリケートしたものです．
　データセットは，wooldridgeパッケージの中に含まれていますので，パッケージを読み込んだ後に，パッケージ内からデータセットを呼び出します．インストールがお済みでない方は，予めインストールをお願いします．

#crime1データが入っているwooldridgeパッケージを読み込み
library(wooldridge)
#crime1データを読み込み
data(crime1)
#データの先頭行を表示
head(crime1)

##   narr86 nfarr86 nparr86 pcnv avgsen tottime ptime86 qemp86 inc86 durat black
## 1      0       0       0 0.38   17.6    35.2      12      0   0.0     0     0
## 2      2       2       0 0.44    0.0     0.0       0      1   0.8     0     0
## 3      1       1       0 0.33   22.8    22.8       0      0   0.0    11     1
## 4      2       2       1 0.25    0.0     0.0       5      2   8.8     0     0
## 5      1       1       0 0.00    0.0     0.0       0      2   8.1     1     0
## 6      0       0       0 1.00    0.0     0.0       0      4  97.6     0     0
##   hispan born60 pcnvsq pt86sq    inc86sq
## 1      0      1 0.1444    144    0.00000
## 2      1      0 0.1936      0    0.64000
## 3      0      1 0.1089      0    0.00000
## 4      1      1 0.0625     25   77.44000
## 5      0      0 0.0000      0   65.61001
## 6      0      1 1.0000      0 9525.75977

用いる変数の定義は以下の通りです．

変数名	定義
narr86	逮捕された回数
pcnv	過去に有罪判決を受けた割合
avgsen	懲役の長さ
tottime	18歳以降の服役回数
ptime86	1986年に刑務所にいた月数
qemp86	4半期のうち従業状態だった期の数
inc86	合法的な所得[$100]
black	黒人ダミー
hispan	ヒスパニックダミー
born60	1960年に生まれダミー

最初に，1年当たり逮捕回数のヒストグラムをみてみましょう．普通の人は年に何度も逮捕されないので，0が極端に多い分布になっていることが分かります．

hist(crime1$narr86)

　早速，以下のコードを実行してデータの読み込みと各種初期設定を行いましょう．

y <- crime1$narr86
X <- cbind(crime1$pcnv,crime1$avgsen,crime1$tottime,crime1$ptime86,crime1$qemp86,crime1$inc86,crime1$black,crime1$hispan,crime1$born60)
colnames(X) <- c("pcnv","avgsen","tottime","ptime86","qemp86","inc86","black","hispan","born60")
maxit=50
tol=1e-08

1行目では，被説明変数として使うnarr86の列をyというオブジェクトに格納しています．
2行目では，cbind関数を用い，今回説明変数として用いる変数の列をひとまとめにして，Xというオブジェクトに代入を行っています．ただ，cbind関数を用いた場合には，どういうわけか列名の情報が落ちてしまうので，3行目でcolnames関数を用い，Xに改めて列名を定義しています．ここまでがデータの読み込みを表す部分になります．
4・5行目では，反復計算を行う上での初期値設定をいくつか行っています．
4行目では，反復計算の最大回数maxitを設定していて，今回はとりあえず50回にしておきました．
5行目では，推定値の収束判定を行うための閾値tolを設定していて，今回は1e-08，すなわち10^(-8)にしておきました．

　続いて，パラメータ推定を行うための各種設定を行っておきます．

y <- as.vector(y)
X <- as.matrix(X)
reg.X <- cbind(1, X)
beta <- rep(0, ncol(reg.X))

1行目では，as.vector関数を用いてyの型をベクトルに，2行目では，as.matrix関数を用いてXの型を行列に変換しました．この変換作業は大概の場合は必要ないのですが，たまにデータフレームやリストといった別の型が指定されている場合があるので，念のため行ったまでです．
3行目では，1列目の要素が全て1，2列目以降が先ほど読み込んだXになるような行列をreg.Xというオブジェクトに格納しました．Rでは，ベクトルの長さが違う場合は短い方のベクトルの要素が循環して使用される（今回は1が循環して使用されている）ので，このような操作が可能になります．
4行目では，今回反復計算によって求める回帰係数推定値の初期値ベクトルbetaを設定します．rep関数は1番目の引数に循環して使用される要素，2番目の引数に循環を何回行うかをしていします．今回は1番目の引数を0に設定しました．
2番目の引数にncol関数が登場していますが，ncol関数は行列の列数を取得する関数で，今回はreg.Xの列数（すなわち説明変数の数）と同じ10になります．
よって，全ての要素が0の10次元の初期値ベクトルが作成されることになります．

　ここからが，反復計算を行う部分のコードになります．

for(k in 1:maxit){
    g <- t(reg.X) %*% (y-exp(reg.X%*%beta))
    W <- diag(nrow(reg.X))
    diag(W) <- exp(reg.X%*%beta)
    I <- t(reg.X) %*% W %*% reg.X
    print(beta)
    old.beta <- beta
    beta <- old.beta + solve(I) %*% g
    if(sqrt(crossprod(beta - old.beta)) < tol)break
}

##  [1] 0 0 0 0 0 0 0 0 0 0
##                 [,1]
##         -0.423433967
## pcnv    -0.131886023
## avgsen  -0.011331585
## tottime  0.012069259
## ptime86 -0.040873459
## qemp86  -0.051309866
## inc86   -0.001461701
## black    0.327009689
## hispan   0.193809364
## born60  -0.022464972
##                 [,1]
##         -0.567827309
## pcnv    -0.301235529
## avgsen  -0.020150186
## tottime  0.021422760
## ptime86 -0.081727140
## qemp86  -0.072874340
## inc86   -0.004281888
## black    0.585251550
## hispan   0.406340432
## born60  -0.044233321
##                 [,1]
##         -0.596735929
## pcnv    -0.389677096
## avgsen  -0.023254354
## tottime  0.024192782
## ptime86 -0.097220695
## qemp86  -0.049704219
## inc86   -0.007118315
## black    0.656668488
## hispan   0.492183273
## born60  -0.050703298
##                [,1]
##         -0.59952307
## pcnv    -0.40134145
## avgsen  -0.02375351
## tottime  0.02448375
## ptime86 -0.09855834
## qemp86  -0.03864302
## inc86   -0.00802895
## black    0.66087303
## hispan   0.49975697
## born60  -0.05102785
##                 [,1]
##         -0.599588647
## pcnv    -0.401571002
## avgsen  -0.023772255
## tottime  0.024490352
## ptime86 -0.098558486
## qemp86  -0.038020462
## inc86   -0.008080557
## black    0.660837837
## hispan   0.499813376
## born60  -0.051028583
##                 [,1]
##         -0.599588795
## pcnv    -0.401571271
## avgsen  -0.023772299
## tottime  0.024490364
## ptime86 -0.098558447
## qemp86  -0.038018715
## inc86   -0.008080704
## black    0.660837581
## hispan   0.499813275
## born60  -0.051028583

1行目では，ループ変数（何回目の繰り返しかを示す変数）kを1からmaxit（今回は50）の範囲で動かす，すなわち，最大maxit回反復計算を行うことを宣言しています．
2行目では，式(24)で示されるポアソン回帰のグラディエントをコーディングしています．
t()は行列・ベクトルの転置，%*%は行列・ベクトルの積を表す演算子です．
3行目では，式(31)で示されるポアソン回帰のヘシアンを構成する行列のうち，対角行列$\mathbf{W}$を作成するため，単位行列を作成するための関数diagを用いて単位行列$\mathbf{W}$を作成しています（わかりにくくて申し訳ないですが）．diag関数の引数は，単位行列の次数です．$\mathbf{W}$の次数はサンプルサイズ，すなわちreg.Xの行数と同じなので，ncol関数でreg.Xの行数を取得し，引数として用いています．
4行目では，単位行列$\mathbf{W}$の対角項に${\rm exp}(\mathbf{x}_{i}'\mathbf{\beta})$を当てはめ，対角行列$\mathbf{W}$を作成しています．対角行列$\mathbf{W}$の作成が完了したので，5行目で情報行列の計算を行っています． 6行目では，forループを抜けた後に，各反復回で得られた回帰係数の推定値を表示するために，先ほど作成したベクトルbetaに対してprint関数を適用しています．
7行目では，k回目の計算で得られた回帰係数の推定値を，old.betaというオブジェクトに格納しました．これが式(32)における$\mathbf{\beta}_{k}$に対応します．
8行目では，式(32)をそのままコーディングしています．solve関数は行列の逆行列を計算する関数です．よって，betaは$\mathbf{\beta}_{k}$，solve(I)は$[\mathbf{I}(\mathbf{\beta}_{k})]^{-1}$，gは$\mathbf{g}(\mathbf{\beta}_{k})$に対応します．
9行目では，推定値の収束判定を行っています．収束判定の条件式は，

\[ \sqrt{(\mathbf{\beta}_{k}-\mathbf{\beta}_{k+1})'(\mathbf{\beta}_{k}-\mathbf{\beta}_{k+1})}<10^{-8} \tag{34} \]

としました．crossprod関数は，ベクトル・行列の内積を計算するための関数です．もしも収束条件を満たした場合には，breakでforループを抜けます．

　以下，各種統計量の計算及び表示の設定を行っていきます．

beta1 <- as.vector(beta)
var.beta <- solve(I)
dimnames(var.beta) <- NULL
stderr <- sqrt(diag(var.beta))
z.value <- beta1/stderr
coef <- cbind(beta1, stderr, z.value)
colnames(coef) <- c("coef", "s.e.", "z-value")
rownames(coef) <- append("(intercept)",colnames(X))
list(coefficient=coef, iteration=k)

## $coefficient
##                     coef       s.e.    z-value
## (intercept) -0.599588795 0.06725010 -8.9158052
## pcnv        -0.401571271 0.08497119 -4.7259698
## avgsen      -0.023772299 0.01994603 -1.1918308
## tottime      0.024490364 0.01475041  1.6603180
## ptime86     -0.098558447 0.02069464 -4.7625102
## qemp86      -0.038018715 0.02902421 -1.3098966
## inc86       -0.008080704 0.00104101 -7.7623727
## black        0.660837581 0.07383422  8.9502883
## hispan       0.499813275 0.07392671  6.7609296
## born60      -0.051028583 0.06405181 -0.7966767
## 
## $iteration
## [1] 7

1行目では，後ほどbetaを用いた計算を行う為，反復計算の結果得られたbetaの方をベクトルに変換し，beta1というオブジェクトに格納しています．
2行目では，推定された回帰係数の分散共分散行列を計算し，オブジェクトvar.betaに格納しています．
最尤推定された回帰係数の推定量$\hat{\mathbf{\beta}}$の漸近分布は，一般に次のように与えられます．

\[ \hat{\mathbf{\beta}} \sim N[\mathbf{\beta},[\mathbf{I}(\mathbf{\beta})]^{-1}] \tag{35} \]

よって，回帰係数の推定量$\hat{\mathbf{\beta}}$̂の分散共分散行列$Var(\hat{\mathbf{\beta}})$は，

\[ Var(\hat{\mathbf{\beta}})=[\mathbf{I}(\mathbf{\beta})]^{-1}=(\mathbf{X}'\mathbf{W}\mathbf{X})^{-1} \tag{36} \]

によって計算することができます．ただ，$\mathbf{W}$の中に含まれる${\mathbf{\beta}}$は通常未知ですので，反復計算の結果得られた回帰係数の推定値で置き換えて計算を行います．
2行目に計算された分散共分散行列には，どういうわけか行名と列名がついてしまっていて，そのままでは扱いにくいので，3行目ではdimnames関数で行・列名にアクセスし，行・列名をNULLに置き換えておきます．
4行目では，回帰係数推定値の標準誤差を計算しています．分散共分散行列は一般に，

\[ Var(\hat{\mathbf{\beta}})= \left[ \begin{array}{ccc} Var(\hat{\beta}_{1}) & \ldots & Cov(\hat{\beta}_{1},\hat{\beta}_{p}) \\ \vdots & \ddots & \vdots \\ Cov(\hat{\beta}_{p},\hat{\beta}_{1}) & \ldots & Var(\hat{\beta}_{p}) \end{array} \right] \tag{37} \]

ですから，その対角項の平方根を取ったものが，各回帰係数の推定値の標準誤差になります．diag関数を用いて分散共分散行列var.betaの対角項にアクセスし，それをstderrというオブジェクトに格納して，標準誤差ベクトルを作成しました．
5行目では，各回帰係数の推定値についてz値を計算し，それらを格納したz.valueというz値ベクトルを作成しました．

\[ z_{i}=\frac{\hat{\beta}_{i}}{\sqrt{Var(\hat{\beta}_{i})}} \tag{38} \]

6行目では，beta1・stderr・z.valueをまとめ，coefというオブジェクトに格納しています．
7・8行目ではcolnames関数を用いてcoefの列名を，rownames関数を用いて行名を指定しています．8行目で出てくるappend関数は2つのベクトルを結合する関数で，今回は ”intercept” という1行1列の文字型ベクトルに，colnames(X)，即ち最初に読み込んだXの列名（説明変数名）を表す文字列ベクトルを結合しています．
9行目では，list関数を用いてcoefと繰り返し回数kを同時に表示するための設定を行っています．ここまでで，回帰係数の推定結果は概ね出そろいました．

　最後に，Rにデフォルトで入っている，ポアソン回帰を実行するglm関数を実行して，答え合わせをしてみてください．

summary(glm(y~X, family=poisson()))

## 
## Call:
## glm(formula = y ~ X, family = poisson())
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.5731  -0.9076  -0.6651   0.2183   7.4609  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -0.599589   0.067250  -8.916  < 2e-16 ***
## Xpcnv       -0.401571   0.084971  -4.726 2.29e-06 ***
## Xavgsen     -0.023772   0.019946  -1.192   0.2333    
## Xtottime     0.024490   0.014750   1.660   0.0969 .  
## Xptime86    -0.098558   0.020695  -4.763 1.91e-06 ***
## Xqemp86     -0.038019   0.029024  -1.310   0.1902    
## Xinc86      -0.008081   0.001041  -7.762 8.34e-15 ***
## Xblack       0.660838   0.073834   8.950  < 2e-16 ***
## Xhispan      0.499813   0.073927   6.761 1.37e-11 ***
## Xborn60     -0.051029   0.064052  -0.797   0.4256    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 3208.5  on 2724  degrees of freedom
## Residual deviance: 2822.2  on 2715  degrees of freedom
## AIC: 4517.5
## 
## Number of Fisher Scoring iterations: 6

optim Function

　ここまでで，スコア法を用いた数値的な求根手法についてみてきましたが，以下では，Rにデフォルトで入っている汎用の最適化関数optimを用いた求根法についても説明しておきたいと思います．
optim関数は主に，解の初期値・最適化を行いたい関数を引数として，関数の最適化を実行する関数です．
ここでは，先ほどスコア法を用いて行ったポアソン回帰のパラメータ推定を，optim関数を用いて同様に解いてみます．その際，以下のようなコーディングを行います．

LL <- function(beta){sum(y*Design%*%beta-exp(Design%*%beta)-lfactorial(y))}
Design <- cbind(1,X)
k <- ncol(Design)
set.seed(123)
ML <- optim(par=rep(0,k),fn=LL,control=list(fnscale=-1),hessian=TRUE,method="SANN")
SE <- sqrt(diag(solve(-ML$hessian)))
cbind(coef=ML$par,s.e.=SE,z.value=ML$par/SE)

##               coef         s.e.   z.value
##  [1,] -0.633992449 0.0676589904 -9.370410
##  [2,] -0.250734746 0.0815081631 -3.076192
##  [3,] -0.103122465 0.0271599892 -3.796852
##  [4,]  0.067353911 0.0163106204  4.129451
##  [5,] -0.100555444 0.0231523961 -4.343198
##  [6,]  0.049562362 0.0283226351  1.749921
##  [7,] -0.007885675 0.0009474701 -8.322875
##  [8,]  0.346327773 0.0767868069  4.510251
##  [9,]  0.127227180 0.0765293292  1.662463
## [10,] -0.074793641 0.0646610512 -1.156703

1行目では，式(23)に示した対数尤度関数を回帰係数betaの関数として定義し，LLというオブジェクトに格納しています．Designは説明変数行列を意味します．
2行目では，分析に用いる説明変数行列を作成しています．スコア法を用いたパラメータ推定を行った際に用いたreg.Xと同じものです．
3行目では，optim関数で最尤法を実行する際に数値的に計算される，回帰係数betaの次数をncol関数で取得し，kに格納しました．
4行目では，optim関数実行に際して用いる乱数を指定します．この乱数を元に，optim関数内で用いる最適化アルゴリズムが動くためです．
5行目が実際に最尤推定を実行している部分になります．optim関数の引数について，ひとつずつ説明していきます．
parは最適化を行う際の初期値を意味します．今回は，最尤推定を行う際のbetaの初期値を代入しています．スコア法でパラメータ推定を行った際と同様，全て0とします．
fnは最適化を行う対象となる関数を意味します．今回は関数LLについて最適化を行いたいので，LLを代入しています．
control=list(fnscale=-1)は，「LLを最大化する」という意味のある種のおまじないです．この記述を行わなければLLは最小化されます．
hessian=TRUEは，「ヘシアン行列を表示する」という意味です．この記述を行わなければヘシアン行列は表示されません．
6行目では，推定された回帰係数の標準誤差を計算しています．手続きはスコア法で推定を行ったコーディングと同じです．比較してみてください．
7行目では，回帰係数の推定値・標準誤差・z値をひとまとめにして表示するような設定を行っています．

　このように，optim関数を用いた最尤推定の場合には，単に最大化したい対数尤度関数の関数形のみを引数として与えればよく，グラディエントも情報行列も明示的に引数として与えずに推定を実行できることが多いです（アルゴリズムによっては必要）．
ただし，スコア法の結果と見比べてもわかる通り，アルゴリズムの違いによって，求められる解が大きく変わってくることに注意が必要です．今回は，最適化アルゴリズムを焼きなまし法“SANN”で与えましたが，準ニュートン法“BFGS”で与えると，スコア法で得られた推定値と極めて近い推定値が得られます．

もしも最尤推定をフルスクラッチで書いたら

Keisuke Takano

4/17/2020