確率・統計

(18) 一般化線形モデル (Generalized Linear Model)

自然科学には様々な法則や関係式が存在します。例えば、バネの復元力はバネの伸び(または縮み)に比例し(「フックの法則 (Hooke's Law)」)、電圧は電流に比例します(「オームの法則 (Ohm's Law)」)。従って、その比例定数がわかっていれば、バネにものをぶら下げた時にどの程度バネが伸びるか、またどの程度の電圧にすれば必要な電流が得られるか、といったことを求めることができます。この比例定数は、あらかじめ実験で得られた測定値を元に求めたり、理論から計算することができます。自然科学だけではなく社会学・経済学・心理学などにおいても経験則などから得られた様々な関係式があり、それを元に値の推測などを行うことができます。
測定値を元に式を求める場合、測定誤差などにより完全に正しい結果を得ることはできません。その誤差はコントロールができないので、ある確率分布に従うと仮定して式を定義し、さらに得られた結果がどれだけ信頼できるのかを推定するための指標として利用します。このように定義した式の中で、値を求めるために使われる変数は「説明変数(独立変数)」といい、式から得られた値を「目的変数(従属変数)」といいます。独立変数に対する係数(比例定数)が従属変数に対する線形式で表され、誤差成分が正規分布に従うとした時の式は「線形重回帰モデル (Linear Multiple Regression Model)」と呼ばれ、その係数は最尤推定量として求めることができます。

今までに紹介したデータ分析法の中の「分散分析 (ANOVA)」や「共分散分析 (ANCOVA)」は、ダミー変数を利用した線形重回帰モデルの一種でした。このように、線形重回帰モデルは広く使われた手法ですが、これをさらに一般化した考え方として「一般化線形モデル (Generalized Linear Model)」が利用されるようになりました。今回は、この「一般化線形モデル」について紹介したいと思います。

1) 一般化線形モデル (Generalized Linear Model)

線形重回帰モデルは、次のような特徴を持っていました。

1. 従属変数 y_i は、独立変数 x_i に対する係数 α の線形式で表される
2. 従属変数 y_i は、期待値 μ_i、分散 σ² の正規分布 N( μ_i, σ² ) に従う

ここで x_i と α はベクトルであり、その要素数は独立変数の数を表します。また、i は独立変数ベクトル x_i そのものの番号を表しています。x_i = ( x_i1, x_i2, ... x_ip )、α = ( α₁, α₂, ... α_p ) とし、α₀ を定数項とすれば、

と表され、ε_i は正規分布 N( 0, σ² ) に従うと考えれば、y_i の期待値 μ_i = E[y_i] は

となります。x_i0 = 1 として x_i と α に定数項も追加すれば、

と非常にシンプルな形で表現することができるようになり、このとき y_i は正規分布 N( μ_i, σ² ) に従います。

独立変数ベクトルが N 個あれば、上記のような線形式も N 個作ることができます。これらを縦に並べると、

となり、μ = E[y] = ( E[y₁], E[y₂], ... E[y_N] )^T、X = ( x₁^T, x₂^T, ... x_N^T )^T とすれば

と表すことができます。X は N 行 p + 1 列の行列であり、μ は N 個、α は p + 1 個の要素数を持つベクトルです。「重回帰分析」「分散分析」「共分散分析」におけるモデル式は、全てこの式から表すことができるのでした。

測定などの誤差成分は線形式を使って表すことができないため、線形重回帰モデルでは正規分布に従うランダムな値としています。このとき、従属変数も正規分布に従うのは前述した通りです。しかし、一般的に全てのパラメータが正規分布に従うわけではなく、その他の確率密度関数を利用して推定した方がより精度の良い結果が得られる場合も充分に考えられます。また、従属変数と独立変数が線形関係だけで成り立つとは限らないので、さらに一般的なモデル式を採用することによって表現の幅を広げることを検討します。そこで、まずは上記線形モデル式の右辺にあたる線形式が単純に従属変数と一致するのではなく、単調増加(または短調減少)で微分可能なある関数 g(y) によって関連付けられると考え、

という形に拡張します。この関数 g(y) は「連結関数 (Link Function)」といい、g( μ_i ) = μ_i (恒等関数) ならば線形重回帰モデルと一致します。また、y_i が従う確率密度関数も正規分布から拡張して「指数型分布族 (Exponential Family of Distributions)」というより一般化したものが扱えるようにします。詳しくは後述しますが、正規分布は指数型分布族に含まれる確率分布になります。

2) 指数型分布族 (Exponential Family of Distributions)

確率密度関数が以下のような式で表現可能ならば、その確率分布は「(一母数)指数型分布族 (Single-parameter Exponential Family of Distributions)」に属するといいます。

h(y), T(y) は確率変数 y を変数とし、スカラー値を返す関数で、η(θ), A(θ) は確率密度関数が持つパラメータ θ を変数とし、スカラー値を返す関数です。g(θ) = e^-A(θ) とすれば、上式は

となり、逆に h(y) = e^B(y) とすれば

になります。各関数がスカラーを返す限り、y, θ はスカラーだけでなくベクトルであっても問題ありません。η(θ) は「自然パラメータ (Natural Parameter)」と呼ばれ、A(θ) を変数変換により A(η) で表すことが可能であれば、次のように記述することができます。

y がスカラーで T(y) = y のとき、f( y | η ) は「正準形 (Canonical Form)」であるといいます。多くのよく知られた確率分布は、正準形の指数型分布族に属しています。

パラメータ η が複数ある場合、その関数を η(θ) = ( η₁(θ), η₂(θ), ... η_k(θ) ) で表し、η(θ) と同数の y を変数とする関数を T(y) = ( T₁(y), T₂(y), ... T_k(y) ) とすれば

と表すことができます。これは「k-母数指数型分布族 (k-parameter Exponential Family of Distributions)」といいます。

前述の通り、正規分布は指数型分布族に属します。正規分布は

より、η₁ = μ / σ², η₂ = -1 / 2σ² とすれば

となるので

より、h(y) = 1、η = ( μ / σ², -1 / 2σ² )、T(y) = ( y, y² )、A(η) = -η₁² / 4η₂ - log( -η₂ / 2π ) / 2 としたとき、正規分布は 2-母数指数型分布族ということになります。しかし、通常はパラメータは一つとすることが多く、その他のパラメータは「局外パラメータ(Nuisance Parameter)」として既知の定数とみなします。正規分布の場合も、σ² を局外パラメータとすれば、

において η = μ / σ² としたとき、μ² / 2σ² = ( ησ )² / 2 より

となって、h(y) = exp( -y² / 2σ² )、η = μ / σ²、T(y) = y、A(η) = ( ησ )² / 2 + ( 1 / 2 ) log( 2πσ² ) の一母数指数型分布族として扱うことができます。

二項分布は下式で表される離散確率密度関数です。

π は一回の試行において成功する確率を表し、0 から 1 までの値をとるのでした。B_N( r | π ) は、N 回の試行で r 回成功する確率を示します。この式は以下のように変形することができます。

η = log( π / ( 1 - π ) ) のとき、π / ( 1 - π ) = e^η より

なので、二項分布は h(r) = _NC_r、T(r) = r、A(η) = Nlog( e^η + 1 ) の一母数指数型分布族です。

ポアソン分布は下式で表される離散確率密度関数です。

平均 λ 回発生する事象が r 回発生する確率はポアソン分布に従うと考えることができるのでした。この式を変形すると

なので、η = logλ とすれば λ = e^η となって、h(r) = 1 / r!、T(r) = r、A(η) = e^η の一母数指数型分布族であることが分かります。

指数分布は下式で表される連続確率密度関数です。

時間 y の間に発生した事象の回数がポアソン分布に従うとして、一回も発生しない確率を考えると指数分布になります。この式は

となるので、h(y) = 1、η = -λ、T(y) = y、A(η) = -log( -η ) の一母数指数型分布族となります。

次に、( 一母数 η、一確率変数 y の ) 指数型分布族が持つ性質を調べてみたいと思います。

に対して、左辺の θ に対する導関数を計算すると

になります。但し、微分と積分の式が入れ替え可能であることを前提としています。右辺はゼロになるので、η'(θ)E[T(y)] - A'(θ) = 0 より

となり、T(y) の期待値を求める公式が得られます。特に正準形の場合は

とすることができます。A(θ) の変数を η に変数変換して A(η) とした場合は

となるので、

と、非常に単純な式になります。次に、f(y|θ) の二階導関数を求めると、

と求めることができて、E[T(y)] = A'(θ) / η'(θ) より

を代入すると

となります。従って、

より、( ∂² / ∂θ² ) ∫{y∈Ω} f( y | θ ) dy = 0 であることから

となって、T(y) の分散を求める公式が得られます。これも A(θ) の変数を η に変数変換して A(η) とした場合、二階導関数は

になるので、E[T(y)] = A'(η) であることを利用して y による積分を求めると

という結果になります。これらの式を使い、指数型分布族に属する確率分布の期待値と分散を求めると次のようになります。

このように、指数型分布族は、A と η の二つの関数によって T(y) の (特に正準形なら y の) 期待値と分散を求めることができるという特徴があります。

3) スコア統計量 (Score Statistic)

互いに独立な確率変数 y = ( y₁, y₂, ... y_N ) が、同じ型の指数型分布族に従い、さらに確率分布のパラメータはそれぞれの y_i について異なるものとします。すなわち、

が成り立つとした時、y に対する尤度は、

になります。但し、θ = ( θ₁, θ₂, ... θ_N )、η = ( η₁, η₂, ... η_N ) で、θ_i, η_i が y_i に対するパラメータとなります。N 個の確率変数に対してパラメータの数も同数の N 個あるようなモデルは「飽和モデル(Satuated Model)」「最大モデル(Maximal Model)」「フルモデル(Full Model)」などとよばれます。この尤度の対数関数は

なので、∂l(θ|y) / ∂θ_i, ∂l(η|y) / ∂η_i は

となります。この式を「スコア統計量(Score Statistic)」と呼び U ( = U(θ_i) または U(η_i) ) で表します。U は確率変数 y_i に依存しているのでやはり確率変数であり、その期待値は

となります。但し、E[T(y_i)] = A'(θ_i) / η'(θ_i) であることを利用しています。分散 V[U] は、確率変数の一次結合に対する分散の公式から

となります。これを「フィッシャー情報量(Fisher Information)」といいます。

分散と期待値の関係式から

であり、E[U] = 0 なので

が成り立ちます。また、

より

となるので、最終的に

という結果が得られます。つまり、フィッシャー情報量はスコア統計量の二乗の期待値に等しく、またスコア統計量をパラメータで微分した値の期待値の符号を逆転した値とも等しいということになります。

飽和モデルの場合、U は対数尤度関数の導関数そのものなので、U = 0 となるときの θ_i, η_i が最尤推定量であり、その値は、y_i が正準形の指数型分布族に従うと仮定した時

を満たす θ_i, η_i を求めることで得られます。しかし通常の場合、確率変数 y_i の期待値 μ_i = E[y_i] が連結関数 g(μ_i) によって

の形のモデル式で表されることを想定して θ_i や η_i ではなくパラメータ α の推定を行います。ここで、x_i, α の要素数を p としたとき p < N であり、飽和モデルのときより少数のパラメータで μ_i を表すわけです。u_j ≡ ∂l(θ|y) / ∂α_j は、l_i = η(θ_i)y_i - A(θ_i) + B(y_i) とすると l_i は θ_i を変数とする関数であり、μ_i = A'(θ_i) / η'(θ_i)、g(μ_i) = x_i^Tα より μ_i は θ_i, α を変数とする関数なので、

と表すことができて、個々の偏微分を求めると、

となり、これらを代入して

という結果が得られます。この結果からも u_j の期待値 E[u_j] は明らかにゼロであり、u_j と u_k の共分散は E[u_ju_k] から求められ

となります。ここで、y_i が互いに独立なら E[( y_i - μ_i )( y_l - μ_l )] = 0 ( i ≠ l ) であることから、

という結果が得られ、これは x_i ( i = 1, 2, ... N ) の j 番目と k 番目の積の和の形になっているので、

と表すことができます。

より、この対角行列を W で表し、j 番目のパラメータ α_j を係数とする独立変数からなるベクトルを ( x_1j, x_2j, ... x_Nj )^T = x'_j とすれば、上式は

と表されます。よって、これが j 行 k 列目の要素となるような行列を Φ としたとき

となります。ここで、Φ は p 行 p 列の正方行列で「フィッシャー情報行列 (Fisher information matrix ; FIM)」、X は x_i^T を i 行目 ( x'_j を j 列目 ) とする N 行 p 列の行列で「デザイン行列 (Design Matrix)」と呼ばれます。この最も単純な例は、パラメータ α の要素がただ一つのみで、g(μ_i) = μ_i = α₀ (定数) の場合で、このとき x_i = 1、W は 1 / V[y_i] を i 番目の対角成分とする対角行列で、X は全要素が 1 の N 次元ベクトルになります。スコア統計量は

と求められるので、u = 0 のときの α₀ (すなわち最尤推定量) a^₀ は

となり、

が u の分散を表します。

u_j は、y_i の線形結合の形で表されるので、y_i が正規分布に従うならば u_j も正確に正規分布に従います。y_i が正規分布以外の確率分布に従う場合はそうはなりませんが、y_i - μ_i は測定値と母集団平均の差分を表しており、u_j はその和の形となっていることから、中心極限定理より、N が十分に大きければ u_j は漸近的に正規分布に従います。u_j の平均はゼロで、u_j と u_k の共分散 E[u_ju_k] からなる共分散行列はフィッシャー情報行列 Φ そのものなので、u = ( u₁, u₂, ... u_p )^T としたとき、u は多変量正規分布 N( 0, Φ ) に漸近的に従うことになります。このことから、u^TΦ^-1u が自由度 p の χ²-分布に漸近的に従うことになり(補足 1)、その値を使ってモデル式の適合性を検定することができます。

μ_i は μ_i = A'(θ_i) / η'(θ_i) = A'(η_i) より θ_i や η_i に依存すると同時に g(μ_i) = x_i^Tα より α_j にも依存します。また、V[y_i] も V[y_i] = [ A⁽²⁾(θ_i)η'(θ_i) - A'(θ_i)η⁽²⁾(θ_i) ] / η'(θ_i)³ = A⁽²⁾(η_i) より θ_i や η_i に依存することから μ_i と関連性があり、さらに α_j に依存することになります。従って、X と y が与えられた時、y が互いに独立に一母数指数型分布族に従うなら u_j や Φ は N 個の変数 θ, η の関数になり、期待値の連結関数が p 個の係数 α によって Xα の形に表せるなら p 個の変数 α の関数になります。例えば、g(μ_i) = μ_i (恒等関数) で、y_i がそれぞれ独立に正規分布 N( μ_i, σ² ) に従うとした時の飽和モデル式は、

より

となるので、U(μ_i) = 0 となるときの μ_i (つまり μ_i の最尤推定量 m^_i ) は y_i と等しくなります。ここで、μ_i = α₀ (定数) というもっとも単純なモデルを考えると、y_i は正規分布 N( α₀, σ² ) に従い、x は全ての要素が 1 のベクトルなので、

より、α₀ の最尤推定量 a^₀ は m_y になります。但し、m_y は y_i の標本平均を表します。u₀ の期待値 E[u₀] は、E[m_y] = α₀ からゼロであり、分散 V[u₀] は

という結果になります。但し、y_i は互いに独立としているので、E[( y_i - α₀ )( y_l - α₀ )] = 0 ( i ≠ l ) であることを利用しています。u₀ は正規分布に従う確率変数 y_i の線形結合であり、やはり正規分布に従うことから、u₀ は N( 0, N / σ² ) に従うという結果が得られます。従って、

とすれば、z は標準正規分布 N( 0, 1 ) に従い、これを使った平均値の推定を行うことができます。これは、母集団の分散 σ² が既知の場合の平均値推定と全く同じ内容です (*3-1)。

μ_i = α₀ + α₁x_i1 というモデル式に対しては、x_i = ( x_i0, x_i1 )^T = ( 1, x_i1 )^T なので、

となるので、u₀ = u₁ = 0 を満たす時の解 α₀, α₁ は連立方程式

を解くことによって求めることができます。但し、m_x, m_y はそれぞれ x_i1, y_i の標本平均、m_xx, m_xy はそれぞれ x_i1², x_i1y_i の標本平均を表しています。これは、単回帰式の係数に対する最尤推定量を求めるときの方程式と全く同じものであり (*3-2)、その解は

になります。但し、s_xy は x_i1, y_i の標本共分散、V_x は x_i1 の標本分散をそれぞれ表します。

独立変数を増やして μ_i = x_i^Tα とします。但し、x_i の第一番目の要素は 1 とし、定数項を表します。このときも V[y_i] = σ²、g'(μ_i) = 1 なので W = ( 1 / σ² )E ( E は単位行列 ) となって、フィッシャー情報行列は

と非常に単純な形になります。また、u_j は

となり、u_j = 0 となるときの α が最尤推定量になります。この式は、「線形重回帰モデル」で得られた正規方程式と全く同じ形になり (*3-3)、

という結果が得られます。この正規方程式は、定数項となる一番目の要素を消去して式を一つ減らすと、(定数項を除いた場合の) X の標本共分散行列を V、x_i の j 番目の要素からなる N 次元ベクトル x'_j と y の標本共分散を要素とするベクトルを s_y としたとき

と表されるのでした。従って、V が逆行列を持てば、α の最尤推定量は

を解くことで求められます。正規分布のように、V[y_i] を定数として扱うことができるのであれば、線形重回帰モデルで得られるような形の正規方程式を解くことによって解析的に α を得ることができます。しかし、V[y_i] が y_i や θ_i、η_i によって変化すると式は複雑になり、さらに連結関数が恒等関数でない場合は線形式とはならず、簡単には解が得られなくなります。

*3-1) 「確率・統計 (8) 推定」の「4-1) 母集団の平均値の区間推定」に関する説明を参照

*3-2) 「確率・統計 (11) 二標本の解析 -1-」の「2) 相関係数と回帰曲線」に関する説明を参照

*3-3) 「確率・統計 (13) 回帰分析法」の「1) 線形重回帰モデル(Linear Multiple Regression Model)」に関する説明を参照

4) スコア法 (Method of Scoring)

一般に、α_j の最尤推定量 a^_j を求めるためには、p 個の正規方程式

を解く必要があります。ここで、m^_i は、a^ = ( a^₁, a^₂, ... a^_p ) として

を満たす N 個の変数であるため、未知数は p 個の a^_j であり、その値は一意に決めることができます。しかし、前節で述べたように、g( m^_i ) が恒等関数でないかぎり正規方程式は単純な線形式にはなりません。例えば、g( m^_i ) = log( m^_i ) ならば

より正規方程式は

となって、簡単に解を得ることはできなくなります。そこで、解析的に解くのではなく、数値計算アルゴリズムを使って数値解を求める方法を利用することになります。

「ニュートン - ラフソン法 (Newton-Raphson Method)」は、以下のような漸化式を使って f(x) = 0 となる x を求める手法でした (*4-1)。

但し、f(x) は微分可能で、その導関数が既知である必要があります。スコア統計量にこの漸化式を適用すると、

となります。例えば、r_i ( i = 1, 2, ... N ) が互いに独立で、同一のパラメータを持つポアソン分布

に従う場合、対数尤度関数は

なので、スコア統計量 U(η) は

であり、その導関数 U'(η) は

になります。よって、m_r = Σ_i{1→N}( r_i ) / N とすれば、漸化式は

になります。

スコア統計量 U と情報量 V[U] の間には、

の関係があるのでした。近似的に U' ≈ E[U'] = -V[U] が成り立つと仮定すれば、漸化式を

とすることができます。この方法を「スコア法 (Method of Scoring)」といいます。先ほどのポアソン分布の例では、

であり、E[r_i] = V[r_i] = e^η であることと、r_i が互いに独立であることを利用すると、

と求められるので、漸化式は変化はありません ( U'(η) = -Ne^η だったので、U' には確率変数 y_i がなく、E[U'] = U' は正確に成り立ち、計算をしなくても U' = -V[U] は成り立っていることがわかります )。

以下の表は、スコア法を使ってパラメータを求めた結果を示したものです。

m_r は標本平均を表し、その 1 / 2 を η の初期値として漸化式に代入して計算を行なっています。求めた値はやがて収束し、前回との差分は無視できるほど小さくなります。その結果は、m_r = 3 のとき η = 1.099、m_r = 5 のとき η = 1.609、そして m_r = 10 のとき η = 2.303 になります。ポアソン分布は、

とも表されるのでした (というよりもこちらの方が一般的な表し方です)。λ = e^η だったので、上記結果を元に λ を求めると、m_r = 3 のとき λ = 3.00、m_r = 5 のとき λ = 5.00、m_r = 10 のとき λ = 10.00 になり、λ は m_r と等しくなります。

パラメータが一つだけならば、正規方程式も一つのみであり、スコア法を使わなくても通常は簡単に U(η) = 0 の解を得ることができます。上記の例では

となって、スコア法での結果と一致します。スコア法が効果を発揮するのは、パラメータが複数になったときです。

ニュートン - ラフソン法は次のように考えることもできます。テイラー - マクローリン展開から、x が b に非常に近いときは

と近似することができます。f(x) = 0 のとき、両辺を f'(b) で割って

より

となって、ニュートン - ラフソン法での漸化式と一致します。同様に、多変量のテイラー - マクローリン展開 (補足 2) から、x = ( x₁, x₂, ... x_p )^T が b = ( b₁, b₂, ... b_p )^T に非常に近いときは、p 変数を持つ関数 f_j(x) は

と近似することができます。ここで、∇f_j(b) は ∂f_j(b) / ∂b_k ( k = 1, 2, ... p ) を要素とする p 次元のベクトルとします。

f_j(x) = 0 を満たす x の値はこの式だけで一意に決まるとは限りません。例えば二変数の場合、f_j(x) が xy-平面と交差すれば、交差した部分の曲線が f_j(x) = 0 を満たす x の軌跡となります。特に、上記近似式は一次式なので、交差した部分は ( b の近傍では ) 直線になります。もし、一点だけで接していれば ( b の近傍で xy-平面と重なれば ) 解は一意に決まり、全く接していない場合は ( b の近傍で xy-平面と平行な状態で ) 解が存在しないことになります。同様な式が p 個存在すれば、それらを使った連立方程式を解くことで、f_j(x) = 0 ( j = 1, 2, ... p ) を満たす解 x を得ることができます。この連立方程式は、

となるので、f_j(x), f_j(b) ( j = 1, 2, ... p ) を要素とする p 次元のベクトルをそれぞれ f(x), f(b)、j 行 k 列番目の要素を ∂f_j(b) / ∂b_k とする行列を H で表せば、上式は

になります。f(x) = 0 のとき、両辺に左側から H^-1 を掛ければ

という結果が得られ、これが多変量関数に対するニュートン - ラフソン法の漸化式になります。x を α = ( α₁, α₂, ... α_p )、f(x) を u(α) = ( u₁, u₂, ... u_p )、と置き換えれば、前に求めた α の値 α_(m-1) を使って

の漸化式から α_(m) を求めることができます。ここで、H_(m-1) の j 行 k 列目の要素は ∂u_j / ∂α_k = ∂²l / ∂α_j∂α_k になるので、これを計算すると

という結果が得られますが、∂u_j / ∂α_k ≈ E[∂u_j / ∂α_k] が成り立つと仮定すると、

となって、フィッシャー情報行列の符号を反転した -Φ の j 行 k 列目の要素と等しくなります。従って、漸化式は

と近似され、これが複数のパラメータを持つ場合のスコア法の漸化式となります。両辺に左側から Φ_(m-1) を掛けると

となります。Φ は、x_ij を i 行 j 列目 の要素とする N 行 p 列の行列を X、w_i = 1 / V[y_i]g'(μ_i)² を対角要素とする対角行列を W とした時

と表すことができたので、右辺は

であり、X^TW_(m-1)X の j 行 k 列目の要素は、j 番目のパラメータ α_j を係数とする独立変数からなるベクトルを x'_j = ( x_1j, x_2j, ... x_Nj )^T として

だったので、X^TW_(m-1)Xα_(m-1) の j 番目の要素は

となって、右辺の j 番目の要素は

で求められます。

とすれば

となって、漸化式は最終的に

となります。

テーラー - マクローリン展開を利用して、スコア統計量は以下の形に近似することができるのでした。

但し、a = ( a₁, a₂, ... a_p ) は α に非常に近い値をとるものとし、H は j 行 k 列番目の要素を ∂u_j(α) / ∂α_k とする行列です。∂u_j(α) / ∂α_k ≈ E[∂u_j(α) / ∂α_k] のとき、H はフィッシャー情報行列の各要素の符号を反転した -Φ を使って

になります。ここで、a は α の最尤推定量である (よって、a ≈ α である) と仮定すると、a は尤度を最大にする推定量なので、その対数の導関数である u(a) は 0 であり、上式は

と変形することができます。この結果より、Φ^-1 が定数なら、u(α) が漸近的に正規分布に従うことから a も正規分布に従うことになります。右辺の a - α の期待値 E[a - α] は

なので、

となって、a は α の不偏推定量になります。また、a の共分散行列 E[( a - α )( a - α )^T] は

と計算できますが、u(α)u(α)^T の j 行 k 列目の要素は u_ju_k なので E[u(α)u(α)^T] = Φ であり、さらに Φ^-1 は対称行列なので (Φ^-1)^T = Φ^-1 となって、

という結果が得られます。よって、a は多変量正規分布 N( α, Φ^-1 ) に漸近的に従い、

は自由度 p の χ²-分布に漸近的に従うことになります (補足 1)。

α の最尤推定量はスコア法を使って求めることができます。その値は多変量正規分布に漸近的に従う確率変数とみなすことができるので、それを利用して推定や検定を行うことができます。線形重回帰モデルでは、W の対角成分が全て従属変数 y_i の分散 σ² の逆数に等しいことから、フィッシャー情報行列は

であり、α の最尤推定量を求めるための漸化式は、x_iα_(m-1) = μ_i^(m-1) であることを利用すると、z_(m-1) の第 i 成分が

となることから

と求められ、線形重回帰モデルにおける正規方程式と一致します。

であり、定数項を j = 0 として係数 α = ( α₀, α₁, ... α_p ), a = ( a₀, a₁, ... a_p ) としたとき、

なので、

という結果が得られます。但し、j = 0 の項は定数なので、x_i,0 は全て 1 とします。この式の中で、Σ_j{0→p}( ( a_j - α_j )x_i,j ) の部分は

と変形できるので、最終的に

という結果になります。ここで、y^_i は α の最尤推定量を使って求めた y_i の予測値で、μ_i は y_i の期待値であり、真の値ともいうべきものになります。よって、線形重回帰モデルの場合、予測値と期待値の差の二乗和を分散で割った値は自由度 p + 1 の χ²-分布に従うという結果が得られます。特に、p = 0 のときは係数は α₀ のみであり、この値が y_i の期待値 μ になります。その最尤推定量は従属変数 y_i の標本平均 m_y に等しくなるので、上式は

となって、y_i が、同一の正規分布 N( μ, σ² ) から抽出した N 個の標本のとき、上式は自由度 1 の χ²-分布に従うので、先ほどの結果とも一致します(*4-2)。

今度は、スコア法を具体的に算出する手順を確認していきます。

において、両辺はともにパラメータ数 p を要素数とするベクトルで表されます。以下、一部を除いて添字の (m) や (m-1) は省略するので、左辺にある α だけが新たに得られる結果であり、その他は全て前回求めた α から算出されることに注意して下さい。上式の j 番目の要素は、左辺の場合 X^TWX = Φ であり、その j 行 k 列目の要素が x'_j^TWx'_k だったので、

となり、右辺の j 番目の要素は

となります。従って、左辺の α_k の係数 Σ_i{1→N}( x_ijx_ik / V[y_i]g'(μ_i)² ) を計算して係数行列を求め、右辺の値を計算して最終的に得られる連立方程式を解く操作を、得られた未知数 α_k が収束するまで繰り返すことによって解が得られます。

ここで問題になるのが V[y_i] と g'(μ_i) の値の計算方法です。前回求めた未知数 α_k^(m-1) を使い、μ_i は

の関係式から求められます。しかし、μ_i は g(μ_i) の逆関数 g^-1(ξ_i) から求める必要があるので、この逆関数が既知である必要があります。もしそうでないなら、

として、f(μ_i) = 0 となるときの μ_i を「二分法 (Bisection Method)」 (*4-3) や「ニュートン - ラフソン法 (Newton-Raphson method)」 (*4-1) で数値解として求める必要があります。μ_i が得られれば、この値を使って g'(μ_i) も求めることができます。但し、g(μ_i) の導関数が既知でない場合は、「差分法 (Finite Difference Method)」などを利用して近似解を計算することになります。V[y_i] と μ_i = E[y_i] はどちらも θ_i や η_i の関数なので、μ_i を元に θ_i や η_i を求め、それを代入して V[y_i] を求めます。η_i を例にとれば、μ_i = A'(η_i) だったので、A'(η_i) の逆関数が既知でなければ数値解法が必要になります。さらに、A'(η_i) が既知ではないときは、計算に差分法などを用いることになります。V[y_i] = A⁽²⁾(η_i) であり、今度は二階導関数が必要になるので、これも既知でなければ数値解法で求めることになります。

スコア法を使って回帰係数を求めるためのサンプル・プログラムを以下に示します。

/*
  DiffMethod : 差分法による導関数計算用関数

  F* f : 関数 f(x)
  double x : x の値
  double h : x の近傍を決めるための微小量

  計算には中心差分を使う

  f'(x) = [ f( x + h ) - f( x - h ) ] / 2h

  戻り値 : f'(x) の値
*/
template<class F>
double DiffMethod( F* f, double x, double h = DEFAULT_THRESHOLD )
{
  return( ( (*f)( x + h ) - (*f)( x - h ) ) / ( h * 2 ) );
}

/*
  Newton : ニュートン-ラフソン法により f(x) = y を満たす x を求める

  double y : y の値
  F* f( double ) : 関数 f(x)
  DF* df( double ) : 導関数 f'(x)
  double xInit : x の初期値
  double threshold : しきい値

  戻り値 : f(x) = y となる x の値 ( 計算ができなくなった場合は NaN を返す )
*/
template<class F, class DF>
double Newton( double y, F* f, DF* df, double xInit = 1, unsigned int maxCount = DEFAULT_MAX_COUNT, double threshold = DEFAULT_THRESHOLD )
{
  double x = xInit;

  for ( unsigned int i = 0 ; i < maxCount ; ++i ) {
    double curY = (*f)( x ) - y;
    double curDY = (*df)( x );
    if ( isnan( curY ) || isnan( curDY ) ) return( NAN );
    if ( curY == 0 ) return( x );
    if ( curDY == 0 ) return( NAN );

    double diff = curY / curDY;
    x -= diff;
    if ( fabs( diff ) < threshold )
      return( x );
  }

  return( NAN );
}

/*
  MemFunc : メンバ関数を変数 x のみを引数とする関数オブジェクトとして扱う

  Res : 関数の戻り値の型
  T : メンバ関数を実行する対象のインスタンス
  Arg : 引数の型
*/
template<class Res, class T, class Arg>
class MemFunc
{
  const T& t_; // 対象のインスタンス
  Res ( T::*f_ )( Arg ) const; // メンバ関数へのポインタ

public:

  /*
    コンストラクタ

    const T& t : 対象のインスタンス
    Res ( T::*f )( arg ) const : 対象のメンバ関数ポインタ
  */
  MemFunc( const T& t, Res ( T::*f )( Arg ) const )
    : t_( t ), f_( f ) {}

  /*
    operator() : 関数の実行

    Arg x : 対象の引数
  */
  Res operator()( Arg x )
  { return( ( t_.*f_ )( x ) ); }
};

/*
  ExpFamily_IF : 一母数指数型分布族用 I/F
*/
struct ExpFamily_IF
{
  virtual double a( double eta ) const = 0; // A(η)

  // A(η) の導関数 = 期待値 ( E[T(y)] = A'(η) )
  virtual double average( double eta ) const = 0;

  // A(η) の導関数の逆関数 ( η = A'^(-1)(E[T(y)]) )
  virtual double aveInv( double mu ) const = 0;

  // A(η) の二階導関数 = 分散 ( V[T(y)] = A''(η) )
  virtual double variance( double eta ) const = 0;

  // 属性を表す文字列
  virtual string ident() const = 0;
};

/*
  ExpFamily_Generic : 汎用一母数指数型分布族

  導関数・二階導関数の算出に DiffMethod を使用する
*/
class ExpFamily_Generic : public ExpFamily_IF
{
  static const double DEFAULT_H_ = 1E-6; // DiffMethod用の極少量Default値
  double h_; // DiffMethod用の極少量

public:

  ExpFamily_Generic( double h = DEFAULT_H_ ) : h_( h )
  {
    if ( h_ <= 0 ) {
      cerr << "Specified h [" << h_ << "] must be greater than zero.";
      cerr << " Changed to default value [" << DEFAULT_H_ << "]" << endl;
      h_ = DEFAULT_H_;
    }
  }

  // A(η) の導関数 = 期待値 ( E[T(y)] = A'(η) )
  virtual double average( double eta ) const
  {
    MemFunc<double,ExpFamily_IF,double> memFunc( *this, &ExpFamily_IF::a );
    return( DiffMethod( &memFunc, eta, h_ ) );
  }

  // A(η) の導関数の逆関数 ( η = A'^(-1)(E[T(y)]) )
  virtual double aveInv( double mu ) const
  {
    MemFunc<double,ExpFamily_IF,double> pdfExp( *this, &ExpFamily_IF::average );
    MemFunc<double,ExpFamily_IF,double> pdfVar( *this, &ExpFamily_IF::variance );
    return( Newton( mu, &pdfExp, &pdfVar ) );
  }

  // A(η) の二階導関数 = 分散 ( V[T(y)] = A''(η) )
  virtual double variance( double eta ) const
  {
    MemFunc<double,ExpFamily_IF,double> memFunc( *this, &ExpFamily_IF::average );
    return( DiffMethod( &memFunc, eta, h_ ) );
  }
};

/*
  ExpFamily_NormDist : 一母数指数型分布族(正規分布)
*/
class ExpFamily_NormDist : public ExpFamily_IF
{
  double var_; // 分散 σ^2

public:

  // A(η)
  virtual double a( double eta ) const
  { return( ( pow( eta, 2 ) * var_ + log( 2 * M_PI * var_ ) ) / 2 ); }

  // コンストラクタ
  // double sigma : 標準偏差(分散ではないことに注意)
  ExpFamily_NormDist( double sigma ) : var_( pow( sigma, 2 ) ) {}

  // A(η) の導関数 = 期待値 ( E[y] = A'(η) )
  virtual double average( double eta ) const
  { return( eta * var_ ); }

  // A(η) の導関数の逆関数 ( η = A'^(-1)(y) )
  virtual double aveInv( double mu ) const
  { return( mu / var_ ); }

  // A(η) の二階導関数 = 分散 ( V[y] = A''(η) )
  virtual double variance( double eta ) const
  { return( var_ ); }

  // 属性を表す文字列
  virtual string ident() const
  {
    const size_t VAR_SIZE = 9;
    char var[VAR_SIZE + 1];
    snprintf( var, sizeof( var ) / sizeof( var[0] ), "%8.3e", var_ );
    string s = var;
    return( "Normal Distribution (var = " + s + ")" );
  }
};

/*
  ExpFamily_Poisson : 一母数指数型分布族(ポアソン分布)
*/
struct ExpFamily_Poisson : public ExpFamily_IF
{
  // A(η)
  virtual double a( double eta ) const
  { return( exp( eta ) ); }

  // A(η) の導関数 = 期待値 ( E[y] = A'(η) )
  virtual double average( double eta ) const
  { return( exp( eta ) ); }

  // A(η) の導関数の逆関数 ( η = A'^(-1)(y) )
  virtual double aveInv( double mu ) const
  { return( log( mu ) ); }

  // A(η) の二階導関数 = 分散 ( V[y] = A''(η) )
  virtual double variance( double eta ) const
  { return( exp( eta ) ); }

  // 属性を表す文字列
  virtual string ident() const
  { return( "Poisson Distribution" ); }
};

/*
  LinkFunction_IF : 連結関数用 I/F
*/
class LinkFunction_IF
{
public:

  // 連結関数 g(x)
  virtual double operator()( double x ) const = 0;

  // 導関数 g'(x)
  virtual double df( double x ) const = 0;

  // 逆関数 g^-1(y)
  virtual double invf( double y ) const = 0;

  // 属性を表す文字列
  virtual string ident() const = 0;
};

/*
  LinkFunction_Generic : 汎用連結関数

  導関数・二階導関数の算出に DiffMethod を使用する
  逆関数の計算に Newton を使用する
*/
class LinkFunction_Generic : public LinkFunction_IF
{
  static const double DEFAULT_XINIT_ = 1; // xInit_ Default値
  static const unsigned int DEFAULT_MAX_COUNT_ = 1000; // maxCount_ Default値
  static const double DEFAULT_H_ = 1E-6; // h_ Default値
  static const double DEFAULT_THRESHOLD_ = 1E-6; // threshold_ Default値

  double xInit_; // Newton法での x の初期値
  unsigned int maxCount_; // Newton法での反復回数最大値
  double threshold_; // Newton法での収束判定しきい値
  double h_; // DiffMethod用の極少量

public:

  // コンストラクタ
  LinkFunction_Generic( double xInit = DEFAULT_XINIT_, unsigned int maxCount = DEFAULT_MAX_COUNT_,
                        double threshold = DEFAULT_THRESHOLD_, double h = DEFAULT_H_ )
    : xInit_( xInit ), maxCount_( maxCount ), threshold_( threshold ), h_( h )
  {
    if ( h_ <= 0 ) {
      cerr << "Specified h [" << h_ << "] must be greater than zero.";
      cerr << " Changed to default value [" << DEFAULT_H_ << "]" << endl;
      h_ = DEFAULT_H_;
    }
    if ( threshold_ <= 0 ) {
      cerr << "Specified threshold [" << threshold_ << "] must be greater than zero.";
      cerr << " Changed to default value [" << DEFAULT_THRESHOLD_ << "]" << endl;
      threshold_ = DEFAULT_THRESHOLD_;
    }
  }

  // 導関数 g'(x)
  virtual double df( double x ) const
  {
    MemFunc<double,LinkFunction_IF,double> memFunc( *this, &LinkFunction_IF::operator() );
    return( DiffMethod( &memFunc, x, h_ ) );
  }

  // 逆関数 g^-1(y)
  virtual double invf( double y ) const
  {
    MemFunc<double,LinkFunction_IF,double> memFunc( *this, &LinkFunction_IF::operator() );
    MemFunc<double,LinkFunction_IF,double> memDFunc( *this, &LinkFunction_IF::df );
    return( Newton( y, &memFunc, &memDFunc, xInit_, maxCount_, threshold_ ) );
  }
};

struct IFunc : public LinkFunction_IF
{
  // 連結関数 g(x)
  virtual double operator()( double x ) const
  { return( x ); }

  // 導関数 g'(x)
  virtual double df( double x ) const
  { return( 1 ); }

  // 逆関数 g^-1(y)
  virtual double invf( double y ) const
  { return( y ); }

  // 属性を表す文字列
  virtual string ident() const
  { return( "Identity Function" ); }
};

struct LogFunc : public LinkFunction_IF
{
  // 連結関数 g(x)
  virtual double operator()( double x ) const
  { return( log( x ) ); }

  // 導関数 g'(x)
  virtual double df( double x ) const
  { return( 1 / x ); }

  // 逆関数 g^-1(y)
  virtual double invf( double y ) const
  { return( exp( y ) ); }

  // 属性を表す文字列
  virtual string ident() const
  { return( "Logarithm Function" ); }
};

/*
  PrintVector : 可変長配列の要素表示

  const string& header : 最初に出力する文字列
  const vector<double>& vec : 出力対象の配列
*/
void PrintVector( const string& header, const vector<double>& vec )
{
  cout << header << "( ";
  vector<double>::size_type sz = vec.size();
  for ( vector<double>::size_type i = 1 ; i < sz ; ++i )
    cout << vec[i - 1] << ", ";
  cout << vec[sz - 1] << " )" << endl;
}

/*
  PrintMatrix : 二次元可変長配列(行列)の要素表示

  const string& header : 最初に出力する文字列
  const vector<double>& mat : 出力対象の行列
*/
void PrintMatrix( const string& header, const vector< vector<double> >& mat )
{
  vector< vector<double> >::size_type sz = mat.size();
  if ( sz == 0 ) return;
  PrintVector( header, mat[0] );

  string tab( header.length(), ' ' );
  for ( vector< vector<double> >::size_type i = 1 ; i < mat.size() ; ++i )
    PrintVector( tab, mat[i] );
}

/*
  PrintEquation : 方程式の表示

  a を係数とした方程式を表示する

  const string& header : 最初に出力する文字列
  const vector<double>& a : 係数
*/
void PrintEquation( const string& header, const vector<double>& a )
{
  vector<double>::size_type sz = a.size();

  cout << header;
  for ( vector<double>::size_type i = 1 ; i < sz ; ++i )
    cout << a[i - 1] << "x" << i - 1 << " + ";
  if ( sz > 0 )
    cout << a[sz - 1] << "x" << sz - 1 << endl;
}

/*
  ScoringMethod : スコア法による係数の推定

  const vector< vector<double> >& x : 独立変数(p個のパラメータのベクトルからなるn個のベクトル)
  const vector<double>& y : 従属変数
  vector<double>& a : 求めた係数
  const ExpFamily_IF& pdf : 指数型分布族
  const LinkFunction_IF& g : 連結関数
  bool verbose : 冗長モード(ON/OFF)
  unsigned int maxCount : 反復処理の最大回数
  double threshold : 収束条件(全係数が threshold 以下なら処理終了)

  戻り値 : 係数が得られた ... true ; データ異常・反復処理回数が最大値を超えた ... false
*/
bool ScoringMethod( const vector< vector<double> >& x, const vector<double>& y, vector<double>& a,
                    const ExpFamily_IF& pdf, const LinkFunction_IF& g,
                    bool verbose, unsigned int maxCount, double threshold )
{
  cout << "*** Scoring Method ***" << endl << endl;

  if ( &pdf == 0 ) {
    cerr << "Exponential family of distribution not defined." << endl;
    return( false );
  }
  cout << "Exponential Family of Distribution : " << pdf.ident() << endl;

  if ( &g == 0 ) {
    cerr << "Link function not defined." << endl;
    return( false );
  }
  cout << "Link Function : " << g.ident() << endl << endl;

  if ( &x == 0 ) {
    cerr << "x not defined." << endl;
    return( false );
  }
  if ( &y == 0 ) {
    cerr << "y not defined." << endl;
    return( false );
  }
  if ( &a == 0 ) {
    cerr << "Can't get coefficient. a is NULL." << endl;
    return( false );
  }

  unsigned int n = x.size(); // 独立変数ベクトル xi の数
  if ( n == 0 ) {
    cerr << "x has no data." << endl;
    return( false );
  }
  if ( y.size() != n ) {
    cerr << "x size (" << n << ") and y size (" << y.size() << ") not matched." << endl;
    return( false );
  }

  unsigned int p = x[0].size(); // 独立変数ベクトルの要素数
  for ( unsigned int i = 1 ; i < n ; ++i ) {
    if ( x[i].size() != p ) {
      cerr << "x[" << i << "] has different size data (" << x[i].size() << ")." << endl;
      return( false );
    }
  }

  cout << "N = " << n << "; p = " << p << endl << endl;

  if ( verbose ) {
    PrintMatrix( "x = ", x );
    cout << endl;
    PrintVector( "y = ", y );
    cout << endl;
  }

  vector<double> mu( y );  // y の期待値 (yで初期化)
  vector<double> dxi( n ); // g'(mu)
  vector<double> w( n );   // 対角行列の要素
  vector<double> ri( n );  // 右辺の x に掛ける係数

  a.assign( p, 0 ); // 係数の初期化

  LinearEquationSystem<double> s( p ); // 連立方程式計算用インスタンス

  bool isMatched;   // 収束したか
  unsigned int cnt; // 計算回数
  for ( cnt = 0 ; cnt < maxCount ; ++cnt ) {

    if ( verbose ) {
      cout << "----- cnt = " << cnt << " -----" << endl << endl;
      PrintVector( "mu = ", mu );
    }

    // g'(mu) と対角行列の計算
    for ( unsigned int i = 0 ; i < n ; ++i ) {
      dxi[i] = g.df( mu[i] );
      double eta = pdf.aveInv( mu[i] );
      w[i] = pdfVar( eta ) * pow( dxi[i], 2 );
    }
    if ( verbose ) {
      PrintVector( "g'(mu) = ", dxi );
      PrintVector( "w = ", w );
      cout << endl;
    }

    // 係数行列の計算
    for ( unsigned int j = 0 ; j < p ; ++j ) {
      for ( unsigned int k = j ; k < p ; ++k ) {
        s[j][k] = x[0][j] * x[0][k] / w[0];
        for ( unsigned int i = 1 ; i < n ; ++i )
          s[j][k] += x[i][j] * x[i][k] / w[i];
        s[k][j] = s[j][k];
      }
    }

    // 右辺の計算
    for ( unsigned int i = 0 ; i < n ; ++i )
      ri[i] = ( g( mu[i] ) + ( y[i] - mu[i] ) * dxi[i] ) / w[i];
    for ( unsigned int j = 0 ; j < p ; ++j ) {
      s.ans( j ) = x[0][j] * ri[0];
      for ( unsigned int i = 1 ; i < n ; ++i )
        s.ans( j ) += x[i][j] * ri[i];
    }

    if ( verbose ) {
      cout << "Equation System :" << endl;
      s.print();
      cout << endl;
    }

    // 連立方程式の計算
    if ( ! GaussianElimination( s ) ) {
      cerr << "Failed to calculate coefficients." << endl;
      return( false );
    }

    // 各係数が収束しているかを確認する
    isMatched = true;
    for ( unsigned int j = 0 ; j < p ; ++j ) {
      if ( fabs( a[j] - s.ans( j ) ) >= threshold )
        isMatched = false;
      a[j] = s.ans( j );
    }

    if ( verbose ) {
      PrintEquation( "Regression equation : y = ", a );
      cout << endl;
    }

    if ( isMatched ) break;

    for ( unsigned int i = 0 ; i < n ; ++i ) {
      double xi = a[0] * x[i][0];
      for ( unsigned int j = 1 ; j < p ; ++j )
        xi += a[j] * x[i][j];
      mu[i] = g.invf( xi );
    }
  }

  if ( cnt < maxCount ) {
    // 係数行列の再計算(=フィッシャー情報行列)
    for ( unsigned int j = 0 ; j < p ; ++j ) {
      for ( unsigned int k = j ; k < p ; ++k ) {
        s[j][k] = x[0][j] * x[0][k] /
          ( pdf.variance( pdf.aveInv( mu[0] ) ) * pow( g.df( mu[0] ), 2 ) );
        for ( unsigned int i = 1 ; i < n ; ++i )
          s[j][k] += x[i][j] * x[i][k] /
            ( pdf.variance( pdf.aveInv( mu[i] ) ) * pow( g.df( mu[i] ), 2 ) );
        s[k][j] = s[j][k];
      }
    }
    LinearEquationSystem<double> inv( 0 ); // 係数行列の逆行列
    Inverse( s, inv );

    PrintEquation( "Estimated regression equation : y = ", a );
    cout << "variance of a = ( " << inv[0][0];
    for ( unsigned int i = 1 ; i < inv.size() ; ++i )
      cout << ", " << inv[i][i];
    cout << " )" << endl << endl;
  } else {
    cout << "Failed to estimate regression coefficient" << endl << endl;
  }

  return( cnt < maxCount );
}

スコア法を利用するにあたり、導関数と逆関数の算出が必要になります。任意の関数に対して数値解を求めるために、差分法により微分係数を求めるための関数 DiffMethod と、ニュートン-ラフソン法を使って f(x) = y を満たす x を計算するための関数 Newton が用意されています。DiffMethod は中心差分と呼ばれる以下の計算式を使って近似解を求めています。

x を中心に微小量 h だけ正負側にシフトした時の f(x) の値を求め、その二点を結んだ直線の傾きを求めたい微分係数としています。理論的には h → 0 のときの極限が正確な値になりますが、できるだけ h を小さくした値で近似的に求めているわけです。

DiffMethod や Newton は、関数の型をテンプレート引数として渡せるようにして、関数ポインタ・関数オブジェクトのどちらでも利用できるようにしてます。しかし、今回は特にインスタンス内のメンバ関数を呼び出したいので、MemFunc というクラスを定義しています。このクラスは、インスタンス時にメンバ関数へのポインタと対象オブジェクトを受け取り、operator() を使って対象オブジェクトのメンバ関数を呼び出します。DiffMethod や Newton は func() の形で任意の関数の処理ができるようにしているのに対し、メンバ関数を呼び出すためには obj.memfunc() や p->memfunc() の形で呼び出す必要があるため、MemFunc オブジェクトの operator() で変換して、関数ポインタや関数オブジェクトと同じように扱うことができるようにしているわけです。

ExpFamily_IF は (一母数型) 指数型分布族を表すためのインターフェースです。スコア法では、A(η) とその導関数 (平均値)、導関数の逆関数 ( μ から η への変換 )、二階導関数 (分散) を求める必要があるため、それらをメンバ関数 a, average, aveInv, variance として定義しています。ExpFamily_Generic は DiffMethod を使って導関数と二階導関数の近似解を求めるように実装したもので、a を実装すればどのような関数でも値が求められるようにしています。また、正規分布とポアソン分布を表したクラス ExpFamily_NormDist, ExpFamily_Poisson もサンプルとして用意しています。

指数型分布族と同様に、連結関数に対してもインターフェース LinkFunction_IF を用意しています。連結関数では導関数と逆関数が必要となるので、汎用連結関数 LinkFunction_Generic では導関数に DiffMethod、逆関数に Newton を利用しています。また、その他のサンプルとして恒等関数 IFunc と対数関数 LogFunc を実装しています。

スコア法は ScoringMethod で実装されています。非常に長いプログラムですが、内容は先ほど説明した処理を順番に行なっているだけの簡単なものです。連立方程式は LinearEquationSystem クラスで定義して、解の計算は「ガウスの消去法(Gaussian Elimination)」を利用した関数 GaussianElimination を利用しています。全ての係数が収束するまで処理を繰り返し、収束するか、繰り返しの最大回数を超えたらループを抜けるようになっています。最後に求められた係数とその分散値を出力して処理を終了します。先に述べた通り、分散値はフィッシャー情報行列 Φ から得ることができますが、Φ は係数行列の逆行列と等しいため、収束後にもう一度係数行列を計算し、その逆行列を計算すれば Φ を求めることができます。逆行列の計算には、「数値演算法 (7) 連立方程式を解く -1-」の「3) 連立方程式による逆行列の計算」で作成したサンプル・プログラム Inverse を利用しています。なお、ここでは対角成分となる分散値のみを出力していますが、求めた行列の要素そのものを出力すれば、各係数どうしの共分散を得ることもできます。

「確率・統計 (13) 回帰分析法」の章で使ったサンプル・データを、サンプル・プログラムを使って処理してみます。データは次のようなものでした。

二教科 A, B のテストの点数を独立変数 x₁, x₂、全教科を総合評価するための総合テストの得点を従属変数 y として、総合テストの得点と二教科のテストの得点の間のモデル式を定義し、その回帰係数を計算します。線形重回帰モデルの場合、その結果は

でした。スコア法を使い、指数型分布族を正規分布、連結関数を恒等関数として計算した結果は次のようになります ( verbose = true にして冗長モードで出力しています )。但し、正規分布の分散値は y の標本分散 111.35 としています。

*** Scoring Method ***

Exponential Family of Distribution : Normal Distribution (var = 1.113e+02)
Link Function : Identity Function

N = 20; p = 3

x = ( 1, 90, 37 )
    ( 1, 52, 41 )
    ( 1, 80, 49 )
    ( 1, 77, 36 )
    ( 1, 84, 38 )
    ( 1, 86, 25 )
    ( 1, 98, 39 )
    ( 1, 96, 30 )
    ( 1, 52, 81 )
    ( 1, 56, 80 )
    ( 1, 84, 66 )
    ( 1, 89, 69 )
    ( 1, 51, 40 )
    ( 1, 87, 78 )
    ( 1, 68, 56 )
    ( 1, 60, 37 )
    ( 1, 76, 34 )
    ( 1, 57, 28 )
    ( 1, 84, 32 )
    ( 1, 86, 40 )

y = ( 70, 52, 78, 70, 78, 66, 84, 78, 64, 76, 77, 87, 60, 87, 71, 60, 59, 49, 66, 69 )

----- cnt = 0 -----

mu = ( 70, 52, 78, 70, 78, 66, 84, 78, 64, 76, 77, 87, 60, 87, 71, 60, 59, 49, 66, 69 )
g'(mu) = ( 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 )
w = ( 111.35, 111.35, 111.35, 111.35, 111.35, 111.35, 111.35, 111.35, 111.35, 111.35, 111.35, 111.35, 111.35, 111.35, 111.35, 111.35, 111.35, 111.35, 111.35, 111.35 )

Equation System :
(0.179614)x0 + (13.5878)x1 + (8.40593)x2 = 12.5819
(13.5878)x0 + (1069.54)x1 + (626.286)x2 = 971.612
(8.40593)x0 + (626.286)x1 + (449.645)x2 = 604.176

Regression equation : y = 10.3853x0 + 0.560639x1 + 0.36864x2

----- cnt = 1 -----

mu = ( 74.4825, 54.6528, 73.2998, 66.8255, 71.4873, 67.8163, 79.7049, 75.2659, 69.3984, 71.2723, 81.8092, 85.7184, 53.7235, 87.9148, 69.1526, 57.6633, 65.5276, 52.6636, 69.2755, 73.3459 )
g'(mu) = ( 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 )
w = ( 111.35, 111.35, 111.35, 111.35, 111.35, 111.35, 111.35, 111.35, 111.35, 111.35, 111.35, 111.35, 111.35, 111.35, 111.35, 111.35, 111.35, 111.35, 111.35, 111.35 )

Equation System :
(0.179614)x0 + (13.5878)x1 + (8.40593)x2 = 12.5819
(13.5878)x0 + (1069.54)x1 + (626.286)x2 = 971.612
(8.40593)x0 + (626.286)x1 + (449.645)x2 = 604.176

Regression equation : y = 10.3853x0 + 0.560639x1 + 0.36864x2

Estimated regression equation : y = 10.3853x0 + 0.560639x1 + 0.36864x2
variance of a = ( 219.571, 0.0250153, 0.0185105 )

計算結果は線形重回帰モデルと同じになることがわかります。連結関数は恒等関数なので g'(mu) は常に 1 となり、対角行列 W の対角成分も常に分散と同じ値 111.35 となっています。各係数の分散から標準誤差 ( S.E. ) を求めると、

となります。前述の通り、係数ベクトルは多変量正規分布に漸近的に従うので、その場合は信頼度 0.95 の区間が 平均 ± 1.96 x S.E. になることから、各係数の 95% 信頼区間は

になります。

次に、指数型分布族と連結関数に以下のような汎用クラスを使ってみます。

/*
  正規分布を ExpFamily_Generic で表現
*/
struct TestExp : public ExpFamily_Generic
{
  double var_;

  virtual double a( double eta ) const
  { return( ( pow( eta, 2 ) * var_ + log( 2 * M_PI * var_ ) ) / 2 ); }

  virtual string ident() const
  { return( "Test Function" ); }
};
TestExp testExp;
testExp.var_ = 111.35;

/*
  恒等関数を LinkFunction_Generic で表現
*/
struct TestLink : public LinkFunction_Generic
{
  virtual double operator()( double x ) const
  { return( x ); }

  virtual string ident() const
  { return( "Test Function" ); }
};
TestLink testLink;

同じデータで処理を行うと次のような結果になりました。

*** Scoring Method ***

Exponential Family of Distribution : Test Function
Link Function : Test Function

N = 20; p = 3

x = ( 1, 90, 37 )
    ( 1, 52, 41 )
    ( 1, 80, 49 )
    ( 1, 77, 36 )
    ( 1, 84, 38 )
    ( 1, 86, 25 )
    ( 1, 98, 39 )
    ( 1, 96, 30 )
    ( 1, 52, 81 )
    ( 1, 56, 80 )
    ( 1, 84, 66 )
    ( 1, 89, 69 )
    ( 1, 51, 40 )
    ( 1, 87, 78 )
    ( 1, 68, 56 )
    ( 1, 60, 37 )
    ( 1, 76, 34 )
    ( 1, 57, 28 )
    ( 1, 84, 32 )
    ( 1, 86, 40 )

y = ( 70, 52, 78, 70, 78, 66, 84, 78, 64, 76, 77, 87, 60, 87, 71, 60, 59, 49, 66, 69 )

----- cnt = 0 -----

mu = ( 70, 52, 78, 70, 78, 66, 84, 78, 64, 76, 77, 87, 60, 87, 71, 60, 59, 49, 66, 69 )
g'(mu) = ( 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 )
w = ( 111.351, 111.35, 111.347, 111.351, 111.347, 111.351, 111.349, 111.347, 111.351, 111.351, 111.349, 111.349, 111.349, 111.349, 111.351, 111.349, 111.351, 111.35, 111.351, 111.351 )

Equation System :
(0.179614)x0 + (13.5878)x1 + (8.40594)x2 = 12.582
(13.5878)x0 + (1069.54)x1 + (626.288)x2 = 971.615
(8.40594)x0 + (626.288)x1 + (449.646)x2 = 604.178

Regression equation : y = 10.3853x0 + 0.56064x1 + 0.36864x2

----- cnt = 1 -----

mu = ( 74.4825, 54.6528, 73.2998, 66.8256, 71.4873, 67.8163, 79.7049, 75.2659, 69.3984, 71.2723, 81.8093, 85.7184, 53.7235, 87.9148, 69.1526, 57.6633, 65.5277, 52.6637, 69.2755, 73.3459 )
g'(mu) = ( 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 )
w = ( 111.351, 111.351, 111.349, 111.353, 111.347, 111.351, 111.349, 111.351, 111.349, 111.351, 111.351, 111.349, 111.351, 111.351, 111.347, 111.35, 111.349, 111.35, 111.349, 111.349 )

Equation System :
(0.179614)x0 + (13.5878)x1 + (8.40593)x2 = 12.582
(13.5878)x0 + (1069.54)x1 + (626.287)x2 = 971.613
(8.40593)x0 + (626.287)x1 + (449.646)x2 = 604.176

Regression equation : y = 10.3852x0 + 0.56064x1 + 0.368641x2

Estimated regression equation : y = 10.3852x0 + 0.56064x1 + 0.368641x2
variance of a = ( 219.572, 0.0250154, 0.0185105 )

今度は W の対角成分で誤差が少し大きくなっていますが、ほぼ同じ結果が得られています。

*4-1) 「数値演算法 (5) 有理数と無理数の演算」の「6) ニュートン - ラフソン法」を参照

*4-2) 「確率・統計 (6) 標本分布」の「1) χ²-分布(Chi-square Distribution)」にある「χ²-分布に対する性質」を参照

*4-3) 「検索アルゴリズム (2) 二分探索/二分探索木」の「1) 二分探索(Binary Search)」にある「二分法」を参照

5) 尤度比 (Likelihood-ratio)

飽和モデルは、N 個の測定値(独立変数)に対してそれと同数の N 個のパラメータを使ってモデル式を定義していました。各測定値は同じ確率分布を持ち互いに独立であると仮定し、各パラメータは一つの独立変数のみに関連付けられます。但し、繰り返し測定を行ったなどの理由で独立変数が重複する場合、パラメータの数は相異なる独立変数の数と同数とします。飽和モデルに対する対数尤度関数は、指数型分布族が正準形のとき次のように定義されるのでした。

ここで、m は飽和モデルにおける相異なる独立変数の数 ( = パラメータ数 ) を、n_i は i 番目の独立変数において重複した数をそれぞれ表します。∂l / ∂η_i は

なので、∂l / ∂η_i = 0 となるときの η_i は

となって、これを l(η|y) の式に代入することによって飽和モデルにおける最大対数尤度を計算することができます。

y_i の期待値 E[y_i] = μ_i が、連結関数 g(μ_i) によって、独立変数の線形結合で表されるとき、その式を g(μ_i) = x_i^Tα とし、パラメータ数を p とすれば、モデル式は

になります。但し、N は独立変数の(重複も含めた)数を表しています。尤度が最大となるときの η, α をそれぞれ h, a とし、そのときの「尤度比 (Likelihood-ratio)」λ を

と定義すると、この値が大きければ、関心のあるモデル式が飽和モデルに比べてデータをうまく表せていない、すなわち当てはめ方がよくないことを意味します。実際には、尤度比の対数の二倍、

を使ってモデル式の有効性を判断します。この値 D を「対数尤度統計量(Log Likelihood (Ratio) Statistic)」といいます。

上式に対して、テイラー - マクローリン展開を利用した近似式を二次の項まで求めると、最尤推定量 a が変数 α に非常に近いときは

となります。但し、∇l(a) は ∂l(a) / ∂a_j ( j = 1, 2, ... ) からなるベクトルで、H は j 行 k 列の要素を ∂²l(a) / a_ja_k とする行列です。∇l(a) はスコア統計量からなるベクトルなので、a が最尤推定量ならばゼロベクトルになります。また、H はフィッシャー情報行列 Φ そのものなので、

となって、a の要素数 (すなわちモデル式のパラメータ数) を p とすれば、この値は自由度 p の χ²-分布に漸近的に従います (*5-1)。この式を使うと、対数尤度統計量 D は

と表すことができます。第一項は自由度 m の χ²-分布に従い、第二項は自由度 p の χ²-分布に従います。最後の項は、l(α) による当てはめがよいとき l(η) ≈ l(α) となり、無視することができるので、当てはめがよければ D は漸近的に自由度 m - p の χ²-分布に従うことになります。

飽和モデルの最大尤度は、他のモデル式の最大尤度よりも必ず大きくなります。なぜなら、飽和モデルは他のモデルよりも変数が多く、尤度をできるだけ大きくするための自由度も最も高いからです。わかりやすい例では、分散が既知の正規分布に従う独立な N 個の確率変数 y_i ( i = 1, 2, ... N ) に対し、飽和モデルは

なので、μ_i = y_i のとき尤度は最大値 1 / ( 2πσ² )^N/2 になるのに対し、パラメータ数が一つのみの場合は

で、最大値は μ = Σ_i{1→N}( y_i ) / N のときで、全ての y_i が等しくない限り 1 / ( 2πσ² )^N/2 より小さくなります。これは、変数が少ない分、尤度が最大となるようにパラメータを調整する自由度が小さいためです。このことから、l(η) - l(α) は必ず正値をとり、当てはまりが良くなければその値は大きくなっていきます。従って、D が χ²-分布に従うことを利用して、推定や検定を行うことができます。この検定法は「尤比度検定 (Likelihood-ratio Test)」と呼ばれています。

y_i が互いに独立に正規分布 N( μ_i, σ² ) に従う場合を考えると、飽和モデルにおける対数尤度関数は

なので、l(μ) が最大値となるのは μ_i = y_i のときで

になります。連結関数が恒等関数で、g(μ_i) = μ_i = x_i^Tα のときは

であり、α の最尤推定量を a とすれば

です。よって、両者の対数尤度統計量は

であり、x_i^Ta は予測値そのものなのでこれを y^_i とすれば

という結果になります。これは、観測値と予測値の残差の二乗和を分散で割った結果 (残差の平方和) であり、自由度 N - p - 1 の χ²-分布に従うのでした(*5-2)。飽和モデルのパラメータ数は N であり、線形重回帰モデルの場合は(定数項の添字をゼロとした時) p + 1 なので、D が従う確率分布に一致します。また、このとき、D は「漸近的に」ではなく「正確に」自由度 N - p - 1 の χ²-分布に従うことになります。

ここで、p = 0 の場合を考えると、これは μ_i = α₀ すなわち定数であることを意味し、そのときの対数尤度の最大値は、α₀ の最尤推定量が a₀ = m_y すなわち y_i の標本平均だったので

となり、l(μ)_max との間の対数尤度統計量を D₀ とすると

となります。この値は観測値の ( 平均を原点とした時の ) 平方和を分散で割った値 (観測値の平方和) であり、自由度 N - 1 の χ²-分布に従います(*5-2)。p = 0 のときのモデル式は当てはめの最も悪いものであると考えられます。最も単純でパラメータ数も最小であることから、このモデルは「最小モデル (Minimal Model)」と呼ばれています。D₀ は飽和モデルと最小モデルの尤度比なので、最も値が大きくなります。D₀ と D の差が大きければ、D の値は相対的に小さいことになり、当てはめも良くなると考えられます。D₀ と D の差と D₀ との比率を R² として

とすると、分散 σ² を打ち消すことができて、「回帰分析法」で紹介した「決定係数 (Coefficient Of Determination)」「寄与率 (Contribution Ratio)」と全く同じ内容になります (*5-3)。ここで、分子部分は「観測値の平方和」と「残差の平方和」の差であり、この値は予測値とその標本平均との差の平方和 (予測値の平方和) と等しくなります。すなわち、予測値の標本平均を m^_y として

が成り立ちます。なお、R² の平方根は「重相関係数 (Multiple Correlation Coefficient)」と呼ばれます (*5-3)。

また、自由度 M, N の χ²-分布に従う互いに独立な確率変数 t, u に対して

としたとき、F は自由度 ( M, N ) の F-分布に従うため、D₀ - D = Σ_i{1→N}( ( y^_i - m^_y )² ) / σ² が自由度 p の χ²-分布に従うことを利用して(*5-4)、

を使って推定や検定を行うことができます。これらは、y_i に共通な分散 σ² が未知の場合に有効な推定・検定方法となるのでした。この結果から、尤度比検定は、重回帰分析において得られた「決定係数」や「分散分析表」の概念を一般化したものであると考えることができます。

最小モデルに対する対数尤度関数の一般式は次のようになります。

したがって、∂l / ∂η は

で、∂l / ∂η = 0 のとき

なので、これを l(η|y) に代入することで最小モデルに対する最大尤度関数が得られます。

一般的な尤比度検定の使い方は、あるモデル M₁ と、その中から不要と思われるパラメータを除外した、より単純なモデル M₀ を比較するというものです。M₀, M₁ それぞれの独立変数の数を q, p としたとき、M₀ は不要なパラメータが除外されているので p > q が成り立ちます。モデル M₀, M₁ の最大対数尤度をそれぞれ l(a₀|y), l(a₁|y) としたとき ( つまり、a₀, a₁ は最尤推定量を表しています )、両者の対数尤度統計量 D は

と表すことができます。相異なる独立変数ベクトルの数を m とする飽和モデルの最大対数尤度を l(h|y)、飽和モデルと M₀, M₁ の対数尤度統計量を D₀, D₁ とすれば、

であり、M₀, M₁ のデータへの当てはめがよいほど D₀, D₁ の値は小さくなり、それらはそれぞれ自由度 m - q, m - p の χ²-分布に従います。D₀, D₁ によってどちらも当てはめがよいと確認できれば、その差 D₀ - D₁ を求めると

となって、この値は自由度 p - q の χ²-分布に従うので、D の信頼度が高ければ、より単純な M₀ を最適なモデルとして選ぶことができます。前述した「回帰分析法」の中では、D₀ に最小モデルを使い、回帰係数によりパラメータ数を少なくしたモデルとの対数尤度統計量により検定を行なっていることになります。

*5-1) 詳細については、この章の「4) スコア法 (Method of Scoring)」を参照

*5-2) 「確率・統計 (13) 回帰分析法」の「2) 重相関係数 (Multiple Correlation Coefficient)」にある「分散分析表」の説明を参照

*5-3) 「確率・統計 (13) 回帰分析法」の「2) 重相関係数 (Multiple Correlation Coefficient)」参照

*5-4) 「確率・統計 (13) 回帰分析法」の「補足2) 予測値の標本分散」参照

6) 残差統計量 (Residuals)

線形重回帰モデル y_i = x_i^Tα + ε_i では、ε_i はモデル式では説明できない誤差成分であり、N( 0, σ² ) に従う確率変数として扱っていました。y_i の期待値 μ_i = E[y_i] が、α の最尤推定量 a によって当てはめ値 y^_i = x_i^Ta として得られたとき、y_i - y^_i を「残差 (Residual)」といいます。誤差成分 ε_i の標準偏差 σ は通常、未知の値ですが、これが推定量 σ^ として得られたとしたとき、残差を σ^ で割った値

は「標準化残差 (Standardized Residual)」または「スチューデント化残差 (Studentized Residual)」と呼ばれます。仮定したモデルがデータをよく表現しているのなら、残差は独立変数や従属変数とは無関係であり、さらには他のパラメータや測定順番などによっても相関は発生しないはずです。従って、横軸に独立変数や従属変数、測定順番、さらには他の関連するパラメータなどをとり、縦軸を標準化残差としたグラフを描いて関連性がないかをチェックすることによって、モデル式による当てはめのよさを調べることができます。標準化残差は利用する確率分布によって異なり、例えばポアソン分布の場合は E[y_i] = V[y_i] なので、この推定値を y^_i としたとき、標準化残差は

になります。

線形重回帰モデルでの残差 y_i - y^_i を成分とする残差ベクトル e = y - y^ = ( y₁ - y^₁, y₂ - y^₂, ... )^T は、正規方程式

より

であることを利用すると

となります。但し、E は「単位行列 (Identity Matrix)」とします。このとき、e の期待値 E[e] は

であり、共分散行列 V[e] は

を計算することで求めることができます。E[yy^T] は

と変形ができます。但し、ε は ε_i からなる残差ベクトルで、互いに独立であり (従って、E[εε^T] の非対角成分は全てゼロになります)、正規分布 N( 0, σ² ) に従うとします (従って、E[ε] = 0, E[ε^T] = 0^T です)。この式を代入すれば、

となります。X^TX は対称行列であり、その逆行列も対称行列なので(*6-1)、E - X(X^TX)^-1X^T も対称行列です。従って、

となり、E - X(X^TX)^-1X^T はべき等行列になります ( 上記の計算過程から X(X^TX)^-1X^T もべき等行列であることがわかります )。H = X(X^TX)^-1X^T とすれば、

という結果が得られ、その対角成分が y_i - y^_i の分散となることから、H の第 i 番目の対角成分を h_i とすれば

となります。H は「射影行列 (Projection Matrix)」または「ハット行列 (Hat Matrix)」と呼ばれます。

σ は相変わらず未知であるため、その推定量を σ^ とすれば、線形重回帰モデルでの標準化残差 r_i は

で計算することができます。h_i は「てこ比 (Leverage)」といい、y^ = Hy より、てこ比の大きな成分は当てはめ値に対して大きな影響を及ぼすことになります。h_i は必ず正値になり、その和は独立変数ベクトルの成分の数 (以下、これを p とします) に等しいため、実測値の総数を N としたとき、てこ比が p / N の 2 から 3 倍以上であれば注意が必要だとされています (補足 3)。標準化残差とてこ比を組み合わせた指標として、「DEFFITS」と「クックの距離 (Cook's Distance)」があり、それぞれ

で表されます。これらが大きな値を示す場合、その実測値の影響力が大きいと判断することができます。

*6-1) 対称行列の逆行列がやはり対称行列になることの証明は「固有値問題 (2) カルーネン・レーベ展開」の「4) 主成分分析」にて最後の方で紹介しています。

今回は、一般化線形モデルの概念を、主に重回帰分析と比較した形で紹介しました。このモデルを利用した具体的な手法については、次の章で紹介したいと思います。

補足1) 多変量正規分布とカイ二乗分布

x = ( x₁, x₂, ... x_p )^T を、正規分布に従う p 個の確率変数からなるベクトルとします。x の共分散行列を V で表すとその r 行 c 列めの要素は E[( x_r - μ_r )( x_c - μ_c )] であり、さらに x の平均ベクトルを μ = ( μ₁, μ₂, ... μ_p )^T とすれば、x は多変量正規分布 N( μ, V ) に従うことになります。V は対称行列であり、その固有値からなる対角行列を D、固有ベクトルからなる直交行列を Q としたとき V = QDQ^T または V^-1 = QD^-1Q^T が成り立つのでした。従って、

と変換することができます。但し、平均がゼロになるように座標をシフトした上で Q^T によって直交変換したベクトル Q^T( x - μ ) を x' としています。このときのベクトル x' は直交行列 Q^T によって回転・鏡映した状態であることを意味するので(主軸変換)、ベクトルどうしの角度や大きさは変化していません。x' = ( x'₁, x'₂, ... x'_p )^T、D の対角成分(すなわち V の固有値)を λ_j ( j = 1, 2, ... p ) とすれば、

となります。ところで、多変量正規分布 N( 0, V ) は

なので、|V| = |D| であることを利用して

と式を変更することができます。但し、|det(J)| は変数変換によって生じるヤコビアンの絶対値で、その i 行 j 列の要素は ∂x'_i / ∂x_j になります。Q の i 列目の固有ベクトルを u_i = ( u_i,1, u_i,2, ... u_i,p ) としたとき、この値は

となるので、ヤコビ行列は直交行列 Q と等しくなり、その行列式は ±1 なので (*N1-1)、最終的に N( 0, V ) は

となります。これは、平均ゼロ、分散 λ_j の正規分布 N( 0, λ_j ) の同時分布を意味し、このとき x'_j は互いに独立な状態になります。標準正規分布 N( 0, 1 ) に従う互いに独立な確率変数の二乗和は、χ²-分布に従うことから、x'_j / √λ_j の二乗和

は自由度 p のχ²-分布に従い、

より ( x - μ )^TV^-1( x - μ ) が自由度 p の χ²-分布に従うことになります。但し、V は逆行列を持つことを前提としているので、固有値にゼロが含まれている場合はこの式は成り立たなくなります。

特に一変量の正規分布に対しては、上式は

となって、( x - μ ) / σ が標準正規分布に従うことから、その二乗は自由度 1 の χ²-分布に従います。

*N1-1) 「確率・統計 (5) 正規分布」の「補足4) 行列の積の行列式」の中で、|AB| = |A||B| が成り立つことを証明しています。これを使えば、Q^TQ = E (単位行列) より Q の行列式は ±1 になります。

補足2) 多変量のテイラー - マクローリン展開 (Multivariate Taylor-Maclaurin Expansion)

多変量のテイラー - マクローリン展開を証明するためには、多変数関数をいったん一変数関数に置き換える必要があります。それにはまず、多変数関数の合成関数の微分に関する公式を証明します。

一変数関数 f(x) が微分可能で、x は微分可能な関数 x = x(t) が成り立つとした時、微分の定義

より

とおけば ε(h) → 0 ( h → 0 ) であり、上式を変形した式

が成り立ちます。この式 (1) を使って合成関数を変形すると、

と表せるので、δ = hx'(t) + hε₁(h) としたとき

となりますが、h → 0 のとき ε₁(h) → 0 であり、従って δ → 0 なので、ε₂(δ) → 0 も成り立ち

となります。よって、(1) 式と見比べると f'(x)x'(t) = ( df / dx )( dx / dt ) が f( x(t+h) ) の導関数であり、従って

が成り立ちます。(1) 式を見ると、x を a に固定して h を変数とすれば、hε(h) が無視できるとき f( a + h ) は h の一次式であり、その傾きは f'(a) になります。hε(h) を無視して、x = a + h ( つまり h = x - a ) とすれば、(1) 式は

となりますが、これは一変数のテイラー-マクローリン展開を一次項まで表した式と同一になります。これらの式は、f(x) を微小な領域において直線と見立てたときの近似式を表していると考えることができます。これを二変数の場合に拡張してみると、(1) 式は次のように表すことができます。

ここで、f_x( x, y ) と f_y( x, y ) はそれぞれ f の x, y による偏微分です。k = 0 ならば (2) 式は

であり、h = 0 ならば

となって、x, y を固定して一変数関数として考えればこの式は (1) と同等になります。また、δ = ( h, k ), ε = ( ε_x(h), ε_y(k) ) としたとき、「コーシー = シュワルツの不等式」より | ( δ, ε ) | ≤ ||δ||・||ε|| が成り立つことから

であり、ε( h, k ) = ( ε_x(h)² + ε_y(k)² )^1/2 とすれば、ε_x(h) → 0, ε_y(k) → 0 のとき ε( h, k ) → 0 であり、(2) 式は

と表すこともできます。(2) 式において hε_x(h), kε_y(k) が無視できるとすれば、f( x + h, y + k ) は平面の方程式を表し、y を固定した時の x の増分は h、x を固定した時の y の増分は k になります。従って、この式は f( x, y ) を微小な領域において平面見立てたときの近似式を表していると考えることができます。さらに一般に、多変数 x = ( x₁, x₂, ... ) の関数 f(x) における (2) 式は

と表すことができます。但し、

であり、

を意味します。

もう一度、二変数の場合に戻って、いよいよ合成関数の微分について考えます。二変数関数 f( x, y ) の変数 x, y が x( u, v ), y( u, v ) としてやはり二変数 u, v によって表されるとした時、(2) 式を使って

と表すことができます (引数の一部は省略しています)。よって、

と変形することができますが、h → 0 のとき ε_u → 0、k → 0 のとき ε_v → 0 で、h → 0, k → 0 であれば Δx → 0, Δy → 0 なので ε_x → 0, ε_y → 0 が成り立ち、

が成り立ちます。よって、(2) 式と見比べると

が成り立つことが証明されます。一般に、多変数関数 f(x) の変数 x が、他の変数 u によって x_i = x_i(u) で表されるとき、

が成り立つことも、同様なやり方で証明できます。

二変数関数 f( x, y ) に対し、x(t) = x + ht, y(t) = y + kt として、φ(t) を

とします。テイラー-マクローリンの公式から、

であり、合成関数の微分の公式から

となります。但し、f_x( x + ht, y + kt ) = ( ∂ / ∂x )f( x + ht, y + kt ), f_y( x + ht, y + kt ) = ( ∂ / ∂y )f( x + ht, y + kt ) であり、それぞれ f を x, y で偏微分した結果に相当します。このとき、上式を

と表すことがあります。このときの ∂ / ∂x, ∂ / ∂y は「微分作用素 (Differential Operator)」と呼ばれます。

今度は二階微分を求めてみると、

となります。但し、f_xx( x + ht, y + kt ), f_yy( x + ht, y + kt ) はそれぞれ x, y による二階偏導関数、f_xy( x + ht, y + kt ) は x, y の順番で偏微分を、f_yx( x + ht, y + kt ) は y, x の順で偏微分を行った二階偏導関数です。f_xy と f_yx は必ずしも等しくないことに注意して下さい。微分作用素の部分だけ括るようにすれば、

と表すことができます。この結果から、

が成り立つことが予想され、これが成り立つと仮定した時、

なので、帰納法により実際に成り立つことが証明できます。従って、

となり、これが二変数でのテイラー-マクローリン展開になります。同様の考え方によって、p 変数のテイラー-マクローリン展開は

となります。

補足 3) てこ比の性質

「てこ比 (Leverage)」は、射影行列 H = X(X^TX)^-1X^T の対角成分を指すのでした。H は対称行列かつべき等行列であり、E - H ( E は「単位行列 (Identity Matrix)」) も対称行列かつべき等行列です。任意の対称行列 A の i 番めの列ベクトルを a_i としたとき、これは i 番めの行ベクトルでもあります。よって、A² の i 番めの対角成分は || a_i ||² で必ずゼロ以上の値をとり、ゼロになるのは a_i = 0 の場合のみです。A がべき等行列ならば、A² = A なので、A が対称かつべき等ならその対角成分は必ずゼロ以上になります。従って、i 番めのてこ比を h_i としたとき、H が対称行列かつべき等行列であることから h_i ≥ 0 であり、E - H が対称行列かつべき等行列であることから 1 - h_i ≥ 0 が成り立ちます。よって、

であることが示されます。

A = BF となる行列 F が存在するならば、A の列ベクトルは全て B の列ベクトルの線形結合で表すことができます (*n3-1)。B の列ベクトルの中で線形独立なものの数は B の階数であり、その線形結合を列ベクトルとする A の階数は、B の階数より減ることはありますが大きくなることは決してありません (階数が減る簡単な例は、F = 0 のときで、この場合は A もゼロ行列なので、B の階数にかかわらず A の階数はゼロになります)。従って、

すなわち、任意の行列を掛けた場合は階数が小さくなる場合があることを意味します。これを利用すると、射影行列の定義から rank(H) ≤ rank(X) であることがすぐにわかります。しかし、

より rank(X) = rank(HX) ≤ rank(H) でもあり、両者が成り立つためには rank(H) = rank(X) でなければなりません。また、H はべき等行列なので、対角成分の和「トレース(Trace)」は階数と等しく、rank(H) = tr(H) が成り立ちます。tr(H) はてこ比の総和を意味するので、結局

という結果が得られます。X の行ベクトルが各独立変数のベクトルを表し、その成分数が p であるとします。また、そのベクトルが N 個あったとすれば、X は N 行 p 列の行列になります。N > p のとき、X^TX が正則なら、rank(X^TX) = p であり、従って rank(X) = p です。よって、Σ_i{1→N}( h_i ) = p であり、h_i の平均は p / N 程度となります。h_i > 0 なので、p / N の 2 から 3 倍以上のてこ比があるということは、他のてこ比は相対的に小さいことになり、特定の独立変数における影響度が大きいことを意味するので注意が必要になります。

英語の "Leverage" は「(目的を達成するための)力、行動力、影響力、勢力」を意味するとともに「てこの作用、てこ装置、てこ比」の意味もあるようです。てこは、小さな力から大きな力を生み出すおなじみの道具のことですが、意味としては「影響力」という訳の方が合っているような気がします。なぜこのような和訳になったのかは不明です。

*n3-1) 「(12) 二標本の解析 -2-」の「補足2) べき等行列の階数」参照。

◆◇◆更新履歴◆◇◆

◎ 「4) スコア法」において、漸化式の右辺の計算式が誤っていたため修正しました (2013-11-30)

連結関数が恒等関数でないと、計算結果がおかしくなります。サンプル・プログラムも併せて修正しています。

◎ 「スコア法」のサンプル・プログラムを少し見直しました (2013-12-15)

1. 一母数指数型分布族用 I/F の ExpFamily_IF にメンバ関数 aveInv を追加
2. スコア法での eta の計算を、Newton 法から aveInv の利用に変更

η の計算には μ = A'(η) の逆関数 η = A'^-1(μ) を使って求めることができます。変更前は、A'(η) とその導関数 A''(η) を使ってニュートン法から η を得ていましたが、逆関数が既知の場合はそれを使った方が効率も精度もよくなります。ニュートン法による解の計算は、汎用一母数指数型分布族クラスの ExpFamily_Generic で実装しています。

◎ 「6) 残差統計量」を追加しました (2013-12-25)

◎ 「4) スコア法」のサンプル・プログラムを少し見直しました (2014-03-09)

◎ 「4) スコア法」の例にて信頼区間の計算が誤っていたので修正しました (2019-04-27)

前に戻る

タイトルに戻る