ガウス・ニュートン法

非線形モデルにたいして最小二乗法によりデータのあてはめをおこなう方法の一つに、 ガウス・ニュートン法がある。 ここで説明変数 x から 目的変数 y を得る式は非線形であるとする。

説明変数 `x_i` と 目的変数 `y_i` のデータの組が `n` 個与えられているとする。また、 与えられたモデル関数が、m 個のパラメータ `alpha_1, alpha_2, cdots, alpha_n` で表されているとする。以降、これらを `bbalpha` で表す。 目的変数の実測値と、説明変数によりモデルで予測された値の差を残差と呼ぶ。すなわち、 `i` 番めのデータの組に関する残差 `r_i` は次の式で与えられる：
`r_i(bbalpha) = y_i - f(x_i; bbalpha)`。
データがもっともよく当てはまった状態とは、次の式：
`S(bbalpha) = r_1^2 + r_2^2 + cdots + r_n^2`
が最小になった状態と考えられる。そこで、`bbalpha` を変化させて、 もっとも`S(bbalpha)` が小さくなるような反復計算方法を与えるのがガウス・ニュートン法である。 反復の第 k 回におけるパラメータの組を `bbalpha^((k))` と記すことにする。

`bbalpha^((k+1)) = bbalpha^((k)) + (J_g^T J_g)^-1 J_g^T bbr(bbalpha^((k)))`

ここで `J_g` は関数 g に関するヤコビ行列と呼ばれ、 n 行 m 列の行列で、次の形で表される。
`J = (((del g_1(bbalpha))/(del alpha_1),(del g_1(bbalpha))/(del alpha_2),cdots,(del g_1(bbalpha))/(del alpha_m)), ((del g_2(bbalpha))/(del alpha_1),(del g_2(bbalpha))/(del alpha_2),cdots,(del g_2(bbalpha))/(del alpha_m )), (vdots, vdots,,vdots), ((del g_n(bbalpha))/(del alpha_1),(del g_n(bbalpha))/(del alpha_2), cdots, (del g_n(bbalpha))/(del alpha_m)))`

非線形モデルとヤコビ行列の例

ミカエリス・メンテン式

ミカエリス・メンテン式とは、酵素の反応速度v に関する式で、次で与えられる：
`v = (Vmax[S])/(K_m + [S])`
ここで、[S]は基質濃度、Vmax は基質濃度が無限大のときの反応速度である。また、Km はミカエリス・メンテン定数である。 [S]とvは実測可能であるので、Vmax と Km をパラメータ推定により求めることができる。 `alpha_1 = Vmax, alpha_2 = Km `として、 ヤコビ行列の i 番めの行は次で表される：

`((del g_i)/(del alpha_1), (del g_i)/(del alpha_2)) = (x_i / (alpha_2 + x_i), -(alpha_1 x_i) / (alpha_2 + x_i)^2)`

ローレンツ分布

分光学において、スペクトル線の形状を記述するために用いられるのがローレンツ分布である。 電磁波の強度 I を角周波数 x によって表すローレンツ分布は次の式で与えられる：
`I = (h gamma^2 )/ ((x-M)^2 + gamma^2) + b + ax`
ここで `M` は共鳴各周波数、`gamma` は半値幅、`h` はピーク強度、`b + ax` はバックグラウンドの強度である。
参考：奥村晴彦　 [改訂新版]C言語による標準アルゴリズム事典 ポリトープ法の項参照。
`M` を精度良く求めることが、物質の正確な同定につながる。このときのヤコビ行列は、 `M, gamma, h, b, a` をそれぞれ `alpha_1, alpha_2, alpha_3, alpha_4, alpha_5` とおくと、 ヤコビ行列の i 番めの行は次で表される：

`((del g_i)/(del alpha_1), (del g_i)/(del alpha_2), (del g_i)/(del alpha_3), (del g_i)/(del alpha_4), (del g_i)/(del alpha_5)) = ( (2 (x_i -alpha_1)alpha_3 alpha_2^2)/((x_i-alpha_1)^2+alpha_2^2)^2, (2alpha_2 alpha_3)/((x_i-alpha_1)^2 +alpha_2^2) - (2alpha_2^3 alpha_3)/((x_i - alpha_1)^2 + alpha_2^2)^2, alpha_2^2/((x_i - alpha_1)^2 + alpha_2^2), 1, x_i)`

以上の式は誤りである可能性が高い。

血糖曲線

佐藤 總夫：自然の数理と社会の数理　II から引用する。正常者の糖の代謝と糖尿病者の糖の代謝に関して、 アッカーマン（Ackerman）らは血糖値 G の時間的変化に関するモデル（血糖曲線）を次のように立てた。
`G(t) = R exp (-alpha t) sin omega t + G_0`
このときのヤコビ行列は、 `R, alpha, omega, G_0` をそれぞれ `alpha_1, alpha_2, alpha_3, alpha_4` とおくと、 ヤコビ行列の i 番めの行は次で表される：

`((del g_i)/(del alpha_1), (del g_i)/(del alpha_2), (del g_i)/(del alpha_3), (del g_i)/(del alpha_4)) = ( exp (-alpha_2 t) sin alpha_3 t , - alpha_1 t exp (-alpha_2 t) sin alpha_3 t , alpha_1 t exp (-alpha_2 t) cos alpha_3 t , 1)`

以上の式は誤りである可能性が高い。

計算例

では、アッカーマンによる血糖曲線の例で確かめてみる。 前述の佐藤の著書では実測値は数値としては記されていない。視察によって読み取ると、 次の通り推察される。

一方、`bbalpha^((0))` は、上記数値を視察により曲線で内挿し、その特徴的な値から求めてみた。 次の値がもっともらしいと考えられる：

以上から、残差とヤコビ行列を求めることができる。

数式

数式記述は MathJax を用いている。

まりんきょ学問所 ＞ JavaScript 手習い ＞ ガウス・ニュートン法