本稿は，筆者の修士学位論文の中から，非線形最小二乗法（非線形回帰分析；ニュートン法，パターン法，ガウス−ニュートン法）の解法の解説を行っている部分を抜き出して，一般的な非線形最小二乗法の解説としたものである．；keywords = 非線形,最小二乗法,ガウス法,パターン法,ガウス−ニュートン法,回帰分析
   このホームページのトップへ

　

　

非線形最小二乗法

（ニュートン法，パターン法，ガウス−ニュートン法）

　

神山　卓也

　

　

　本稿は私の修士学位論文の局部を抜粋したものであり、非線形最小二乗法（ＮＬＬＳ：ニュートン法，パターン法，ガウス−ニュートン法），その実用の過程で用いる制約付最小二乗法（ＣＬＳ）と各種統計量算出外挿テスト用データ比較，などについて解説を行うものである．非線形最小二乗法の計算用ツールとしては，私が開発した行列数値計算・計量経済分析用統計計算用のＣ言語関数ライブラリー「エクスカリバー (Excalibur)」をご参照いただきたい．ここでは，Ｃ関数ソース・コーディングののlzh形式アーカイブ・ファイルをダウンロード可能である．

　

　

１．非線形最小二乗法とは何か

　

　一般に，ｎ次元列ベクトル空間

　

　

からｍ次元ベクトル空間Ｖ_ｍへの写像ｆが以下の２式を充たすとき，ｆを線形写像という．すなわち，Ｖ_ｎの任意の元ｘ，ｙと任意の数λに対して

　

2)　f(ｘ＋ｙ) ＝ f(ｘ) ＋ f(ｙ)

3)　f(λｘ) ＝ λf(ｘ)

　

が成立することである．最も分かりやすく申せば，一次の式が線形の代表例であると考えれば良い．ここで，パラメタに対して線形であるとは，被説明変数をy_ｔ，説明変数をx_１ｔ，x_２ｔ，パラメタをa_１ｔ，a_２ｔ，a_３ｔ，攪乱項をu_ｔとすると，例えば

　

4)　y_ｔ ＝ a_１ｔ ＋ a_２ｔx_１ｔ ＋ a_３ｔx_２ｔ ＋ u_ｔ

　

のようになることである．ここで注意すべき点は，

　

5)　y_ｔ ＝ a_１ｔ ＋ a_２ｔx_１ｔ／x_２ｔ ＋ a_３ｔx_２ｔ^２ ＋ u_ｔ

　

のような場合，説明変数が線形でなくてもパラメタに対しては線形たり得るが，例えば，

　

6)　y_ｔ ＝ a_１ｔ ＋ （a_２ｔ／a_３ｔ）x_１ｔ ＋ u_ｔ

　

のようなモデルは，変数に対しては線形であるもののパラメタに対しては線形ではないことである．したがって，例えばミクロ経済学の消費者行動基礎理論の実証分析を行う際に、構造方程式としてベルヌイ−ラプラス型に効用指標関数を特定化した場合に誘導型方程式として得られる需要関数は

　

　

という形であり，パラメタに関して線形であるから，線形の重回帰分析を用いてパラメタの値を計測することが可能である．一方，例えば効用指標関数を二次式に特定化した場合の需要関数は

　

8)　q_１　＝｛(a_２２p_１ − a_１２p_２)Ｙ｝／｛a_２２p_１^２−2a_１２p_１p_２＋a_１１p_２^２｝

　　　　　＋　｛a_２p_１p_２ − a_１p_２^２｝／｛a_２２p_１^２−2a_１２p_１p_２＋a_１１p_２^２｝＋ u_ｔ

　

という形になるのでパラメタに対して線形ではない（ちなみに，余談であるが説明変数に対しても線形ではない）．

　

　重回帰分析が使える場合には，

　

9)　β_ｐ＝ (Ｘ'Ｘ)^−１･Ｘ'ｙ

　

の式を用いれば，Ｘ'Ｘが特異行列にならない限り，回帰係数を単純・機械的に求めることができる（ただし，ここではＸ'は転置行列である）．しかし，非線形の回帰分析の際には，線形回帰における9)式のように１回の計算で答えを得られるような式はないので，「何らかの方法」を用いてパラメタの値を動かしながらその度に残差平方和を計算することを繰り返し，その残差平方和が一定値以下になったら収束したものと判断する，という方法を採らざるを得ない．これが非線形回帰推定，非線形最小二乗法，あるいは非線形回帰分析と呼ばれるものである．その非線形最小二乗法の方法の内，ニュートン法，パターン法，ガウス−ニュートン法の３種類方法に基づいたのプログラムを先述のごとく筆者は開発し，(1)で述べた各モデルを３種類の方法で測定し，収束したものの中から最も当てはまりの良いものを選ぶことにした．この節では，実験計画の一環として，この３種類の非線形最小二乗法について解説を行う．

　

　

２．ニュートン法

　

　ニュートン法（Newton method）とは，数値計算プログラムの方程式の解法として有名なニュートン法の考え方を，最小二乗推定に応用するものである．最小二乗法の考え方は，目的関数である残差平方和の極小値を実現するパラメタを真の値の推定値とするものである．このとき，パラメタで残差平方和を偏微分すると0に等しくなるから，これをパラメタ個数分の本数の方程式に置き換えて数値計算のニュートン法で解こうというものである．ニュートン法の概略を，岩田暁一・黒田昌裕「最適値プログラム − FORTRAN �Uによる計算プログラム・シリーズ (6) −」（『三田商学研究』11巻３号，1968年８月）pp.107-112から引用・要約して説明する．

　

　区間［a,b］において，関数f(x)の第２次導関数がf''(x)＞0で，f(a)＞0,f(b)＜0を仮定する．いま，aを定点，a_１を任意の区間［a,b］における点とし，f(a_１)の値をテイラー展開によって近似すると，

　

10) 　f(a_１)＝f(a)＋(a_１−a)f'(a)／1！＋(a_１−a)^２f''(a)／2！

　　　　　　 ＋･･･＋(a_１−a)^ｎ−１f^{（ｎ−１）}(a)／(n-1)!＋(a_１−a)^ｎf^（ｎ）(ξ)／n!

ただし，ξ＝a＋θ(a_１−a)，0＜θ＜1

　

になる．10)式において，ε＝a_１−aとし，εを十分小さくとれば２次以降の高次項を近似的に無視可能と考えることができ，10)式は，

　

11)　 f(a_１)＝f(a)＋εf'(a)

　

と近似できる．このとき，方程式f(ｘ)＝0の根はf(a_１)＝0となるようなa_１であるから，ステップ幅εはf'(a)≠0と仮定すれば，

　

12)　 f(a_１)＝0＝f(a)＋εf'(a)

　　　　　∴ εf'(a)＝−f(a)

　　　　　∴　　　ε＝−f(a)／f'(a)

　

の計算で求めることができる．したがって，

　

13)　　a_１＝a＋ε＝a−f(a)／f'(a)

　

となるa_１を計算し，同様の計算を収束するまで繰り返すと，f(ｘ)＝0の根を得ることができる．問題の仮定からf''(x)＞0だからはf'(x)は単調増加するので，根は１つである．この考え方がニュートン法の基本であるが，これだけでは単なる方程式の数値計算法であって，最小二乗法にはまだ到達していない．そこで，これを最小二乗法に拡張する．一般に，構造方程式を

　

14)　 y＝y(x_１,x_２,･･･,x_ｐ;a_１,a_２,･･･,a_ｎ)

　

とする．ここで，x_１,x_２,･･･,x_ｐは説明変数であり，a_１,a_２,･･･,a_ｎはパラメタである．ここで、我々が行うべきことは，目的関数としての残差平方和

　

　

の値を最小にするようなパラメタa_１,a_２,･･･,a_ｎを求めることである．ただし，ここで，添え字pは，当該（被説明）変数予測値である．このことは目的関数Ｓの極値を求めることにほかならないから，結局，

　

16)　　∂S／∂a_１ ＝ 0，　∂S／∂a_２ ＝ 0，･･････，∂S／∂a_ｎ ＝ 0

　

のように，ｎ本の連立方程式を解く問題になる．ここで，16)式の偏導関数（複数）が12)式の関数fに相当するものであるので，f'になるのは15)式の２階の偏導関数である点が，紛らわしいので注意を要する点である．偏導関数が２階になるので，テイラー展開を２次元に拡張する場合を考える．いま，

　

17)　 df(x,y)＝hf_ｘ(x,y)＋kf_ｙ(x,y)

　　　d^２f(x,y)＝h^２f_ｘｘ(x,y)＋2hkf_ｘｙ(x,y)＋k^２f_ｙｙ(x,y)

　

とおいた場合，２次元に拡張されたテイラー展開は，

　

18)　 f(x＋h,y＋k)＝f(x,y)＋df(x,y)＋(1/2)･d^２f(x,y)＋･･･

　　　　　　　　　　＋（1／(n-1)!）d^ｎ−１f(x,y)＋（1／n!）d^ｎf(x+θh,y+θk)

　

である．18)式において，左記のニュートン近似法にしたがってh,kを十分に小さくとって２次以降の項を無視することにより，h,kの近似式を得ることができる．この18)式を用いて，ｎ個のパラメタに関するｎ本の偏導関数式をそれぞれテイラー展開すると，

　

を導くことができる．ここで，εa_１,･･･,εa_ｎについて解いたものが収束幅になる．この収束幅をパラメタの初期値a_１,a_２,･･･,a_ｎに加えて，新たな19)式を再び解く．これを繰り返すと，非線形最小二乗推定値に収束する．収束の幅としては，或る微係数∂S/∂a_ｉの最下限ε^＊を予め定めておいて，各微係数の値がε^＊を下回ったら収束点とすれば良い．

　

　問題となるのは１次と２次の微係数の与え方である．筆者が開発した統計計算プログラム「エクスカリバー」の中でのニュートン法プログラムは，岩田暁一・黒田昌裕「最適値プログラム − FORTRAN �Uによる計算プログラム・シリーズ (6) −」（『三田商学研究』11巻３号，1968年８月）のフォートラン・プログラムをＣ言語に移植したものであるので，微係数の与え方も岩田・黒田［1968］にしたがい，数学上の微分の定義にしたがっての数値計算で微係数を求める方法を採用した．つまり，

　

20)　　∂S／∂a_ｉ ＝ ｛S(a_ｉ＋δt_ｉ，a_ｊ)−S(a_ｉ，a_ｊ)｝／δt_ｉ

　

で１次の微係数を求め，

　

21)　∂^２S／∂a_ｉ∂a_ｊ

　　　 ＝ ｛S(a_ｉ+δt_ｉ,a_ｊ+δt_ｊ)−S(a_ｉ+δt_ｉ,a_ｊ)−｛S(a_ｉ,a_ｊ+δt_ｊ)−S(a_ｉ,a_ｊ)｝｝

　　　　　／｛δt_ｉ･δt_ｊ｝

　　　 ＝ ｛S(a_ｉ+δt_ｉ,a_ｊ+δt_ｊ)−S(a_ｉ+δt_ｉ,a_ｊ)−S(a_ｉ,a_ｊ+δt_ｊ)＋S(a_ｉ,a_ｊ)｝

　　　　　／｛δt_ｉ･δt_ｊ｝

　

で２次の微係数を求めるというものである．このような計算方法の利点としては，２階の微分が可能でありさえすれば，どのような関数系の計測を行う場合にも，関数の１次微係数も２次微係数も求めずに，関数の原型だけを用いたプログラムを組めることである（もしも，20)式や21)式の方法をとらなければ，例えばdlog|x|／dx＝1/xのごとく，常に導関数と２階の導関数の形状を把握していないと計測ができなくなってしまう．）．一方，20)式や21)式を用いたこのような計算方法の欠点としては，演算精度によっては丸めの誤差を無視できなくなってしまうことである．ニュートン法の収束条件は２次微係数に逆行列が存在すること，すなわち，２次微係数が非特異行列（non-singular matrix）であることであるが，筆者の修士学位論文における消費関数の実証分析での実測経験上でも，他の方法で収束する場合でも，ニュートン法プログラムでは途中で２次微係数行列が特異行列（singular matrix）になるエラーが発生してしまう確率が高かった．このことは，ニュートン法自体の問題ではなく，21)式の計算を行う際の演算の精度の問題が大きいのではないか，と考えている．

　

　

３．パターン法

　

　パターン法（pattern search method）は，目的関数y＝g(θ)の値をθのいろいろな値について直接計算計算しながらg(θ)の最小点を見出していく方法である直接探索法　　　（direct search method）の一種である．直接探索法の長所は，目的関数g(θ)が微分可能でない場合にも適用できるという点である．ここでは，岩田暁一『計量経済学』（有斐閣［経済学叢書］，1982）pp.235-239から引用紹介して，パターン法の説明を行う．

　

　いま，計測したい関数を

　

22) 　y_ｐ＝y（Ｘ,θ）

　

とする．ただし，Ｘは変数の行列で，θは求めるべきパラメタのp次元の行列である．このとき，真のyの値をyとおくと，最小二乗法において最小化したい目的関数は

　

　

である．このg(θ)が，θ*で最小点★をとるものとする．いま，初期点▲のベクトルθ_１を

　

　

とし，これを第１基地点（base point）と呼ぶ．一方，ベクトル

　

　

を設定すると，

　

26)　　b_１ ＋ δ_ｉ

　

は第１気地点をθi軸に沿ってδだけ移動させたものである．いま，t_１０＝b_１を中心にして

目的関数をもっと小さくする点の探索を開始する．ここで，点t_１０＋δ_１での目的関数

g(t_１０＋δ_１)がg(t_１０)（＝g(b_１)）よりも小さく，すなわち，g(t_１０＋δ_１)＜g(t_１０)であればt_１１＝t_１０＋δ_１とする．g(t_１０＋δ_１)≧g(t_１０)であれば点t_１０−δ_１で関数g(t_１０−δ_１)と

g(b_１)とを比べて，g(t_１０)＞g(t_１０−δ_１)であればt_１１＝t_１０−δ_１とし，

g(t_１０)≦g(t_１０−δ_１)であればt_１１＝t_１０とする．i＝2以上の点においても，i＝pになるまで同様の探索を行う．以上を式の形で簡潔にまとめると，

　　　

　　　 t_１ｉ＝ ｔ_{１，ｉ−１}＋δ_ｉ，　　if　g(ｔ_{１，ｉ−１}＋δ_ｉ)＜g(ｔ_{１，ｉ−１})

27)　　　　＝ ｔ_{１，ｉ−１}−δ_ｉ，　　if　g(ｔ_{１，ｉ−１}＋δ_ｉ)≧g(ｔ_{１，ｉ−１})＞g(ｔ_{１，ｉ−１}−δ_ｉ)

　　　　　 ＝ ｔ_{１，ｉ−１}，　　if　g(ｔ_{１，ｉ−１})≦Min［g(ｔ_{１，ｉ−１}＋δ_ｉ)，g(ｔ_{１，ｉ−１}−δ_ｉ)］，

　　　　　　　　　　　　　　　　i＝1,2,･･･,p

　　　

になる．i＝pになるまで探索が終わったら，b_２＝t_１ｐとして，点t_１ｐを第２回探索の基地点とする．これを，パターンの拡張という．そして，今度は，第１基地点b_１から第２基地点b_２方向に，両基地点間の倍の距離をとった点

　

28)　　t_２０＝ b_１＋ 2(t_１ｐ−b_１) ＝ b_１＋ 2(b_２−b_１) ＝ 2b_２− b_１

　

を中心にして，新たな探索を行う．以降，同様にして，第ｋ回における通常の探索

　

　

　　　 t_ｋｉ＝ ｔ_{ｋ，ｉ−１}＋δ_ｉ，　　if　g(ｔ_{ｋ，ｉ−１}＋δ_ｉ)＜g(ｔ_{ｋ，ｉ−１})

29)　　　　＝ ｔ_{ｋ，ｉ−１}−δ_ｉ，　　if　g(ｔ_{ｋ，ｉ−１}＋δ_ｉ)≧g(ｔ_{ｋ，ｉ−１})＞g(ｔ_{ｋ，ｉ−１}−δ_ｉ)

　　　　　 ＝ ｔ_{ｋ，ｉ−１}，　　if　g(ｔ_{ｋ，ｉ−１})≦Min［g(ｔ_{ｋ，ｉ−１}＋δ_ｉ)，g(ｔ_{ｋ，ｉ−１}−δ_ｉ)］，

　　　　　　　　　　　　　　　　i＝1,2,･･･,p

　

を行った後に，第ｋ回終了後のパターンの拡張は以下のようになる．

　

30)　　t_{ｋ＋１，０}＝ b_ｋ＋ 2(t_{ｋ＋１，ｐ}−b_ｋ) ＝ b_ｋ＋ 2(b_ｋ＋１−b_ｋ) ＝ 2b_ｋ＋１− b_ｋ

　

この方法は，拡張時に次の探索点にいくまでのベクトルの長さを２倍とすることにより，探索の速くできること，すなわち収束を速くできことに特徴がある．しかしながら，ベクトルの長さを２倍にすることによって行き過ぎた探索点に到達してしまい，最小点★に未だ到達していないにも関わらず，29)式で示す探索を行うどの点ｔ_{ｋ，ｉ−１}±δ_ｉでも，目的関数g(ｔ_{ｋ，ｉ−１}±δ_ｉ)の値が点ｔ_{ｋ，ｉ−１}のものg(ｔ_{ｋ，ｉ−１})よりも小さくならなくなってしまう危険性がある．このときは，点b_ｋにまで戻り，歩幅b_ｋにδを１／２に縮小して，点b_ｋの周りで再度探索を行う．すなわち，

　

31)　　δ_ｉ^{（ｋ＋１）} ＝ (1/2)･δ_ｉ^ｋ；　　　i＝1,2,･･･,p

　

としてから，再度探索を行うのである．もしも，この縮小した歩幅のもとでもg(ｔ_{ｋ，ｉ−１})の値よりもg(ｔ_{ｋ，ｉ−１}±δ_ｉ)の値が小さくならなければ，さらに31)の短縮を繰り返す．そして，各δ_ｉについて予め設けておいた許容限界ε_ｉを，いずれかのｉで

　

32)　 δ_ｉ^{（ｋ＋１）} ＝ (1/2)･ε_ｉ；　　　　i=1,2,･･･,p

　

が成立したら定常点に達したとみなすものである．

　

　筆者が自分の修士学位論文の研究において二次式効用指標関数に基づく様々な需要関数を計測してみた所感としては，確かにパターン法はニュートン法よりも収束し易いが，他方法でも収束した場合と比べての相対的な当てはまり具合は良くなかった．また，解とは全然異なる値に収束してしまう場合もあった．これは，目的関数の等高線の底が複数ある場合，最も深い底が別にあっても近場の底にはまってしまう危険性が高いことを意味するのではないかと考える．

　

　

４．ガウス−ニュートン法

　

　最後に，ガウス−ニュートン法（Gauss-Newton Method）について，Takeshi Amemiya, Advanced Econometrics, Harvard University Press,1985，pp.139-141と1995年度の筆者出身社会人向け夜間共学大学院・東洋英和女学院大学大学院修士課程社会科学研究科の授業「経済数学」受講録（担当：雨宮健教授［当時］）から引用説明を行う．ガウス−ニュートン法も非線形最小二乗法であるから，最小化ならしめたい目的関数は，ニュートン法やパターン法と同様に残差平方和

　

　

である．いま，βの推測値を添え字pをつけてβ_ｐと表すものとする．33)式のf_ｔ(β)を初期推測値β_ｐ１の周りでテイラー展開すると，

　

　

　

になる．ここで，観測期間をｎ期，計測対象となる式の本数をｌ本，未知の計測したいパラメタをｍ個とすると，34)式の最右辺第２項左側の

　

　

は，ｌ×ｎ行ｍ列のサイズの行列となる．一方，最右辺第２項右側の（β −β_ｐ１）は，　ｍ行１列のサイズの行列（ベクトル）となる．ここで，34)式において，予め値を計算可能な項を左辺に移項すると，

　

　

となる．ここで，左辺は全て計算可能で線形であり，ｌ×ｎ行１列のサイズの行列（ベクトル）となる．右辺は，未知のパラメタβに対して線形であるｌ×ｎ本の式である．したがって，36)式の未知のパラメタβは線形の最小二乗法を用いて求めることができる．この新たな推測値β_ｐ２でβ_ｐ１を置き換えて，収束するまで36)式の回帰計算を繰り返す．収束の判定条件としては，予め設けておいた許容限界εを用いて，

　

　

となれば良い．

　

　ここまでが教科書に載っているガウス−ニュートン法の説明であるが，筆者はさらに計算上の工夫を一つ加えた．それは，線形の最小二乗法の方法である．36)式右辺は線形になるので線形の回帰分析が可能になる，というのがガウス−ニュートン法の特徴であるが，唯一の使い辛さとしては，この右辺は定数項のない線形方程式となる点が上げられる．少なくとも，右辺を行列表記したときに第１列が全て１になるなどということは，偶然の女神のいたずらでもない限りあり得ず，先験的には期待し難い．したがって，そのまま単純回帰や重回帰を適用することは困難であると筆者は考えた．当初は，各行毎に第１列の値で同一行の左右両辺を割ることで，重回帰分析が可能な形に無理矢理持っていってしまおうかとも考えたが，それでは残差の分散均一性がなくなってしまう．そこで，定数項が0のときの回帰分析，すなわち原点を通る回帰分析を行えるようにするべく，制約付最小二乗法（Constaraind　Least Squares；ＣＬＳ）を行うことにした．そこで，岩田暁一『計量経済学』（有斐閣［経済学叢書］，1982）から引用して，制約付最小二乗法について説明する．

　

　Ａをr×(p＋1)の一定値の行列，ｂをr×1の一定値のベクトルとし，Ａの位はρ(Ａ)＝rとする．Ａ，ｂの要素は全て既知であるとしたとき，βが

　

38)　 Ａβ ＝ ｂ

　

という線形制約式を満たしているとする．これが制約条件を表す式であり，ここで，例えば，筆者が行いたい制約条件であるところの，定数項が0である（＝原点を通る回帰）という条件は，［1 0 ･･･ 0］β＝0になる．

　

　この38)式を制約条件として，Ｓ＝(ｙ−ｘβ)'･(ｙ−ｘβ)を最小にする制約付最小二乗法を考える．λをラグランジュの未定乗数を要素とするｒ×１のベクトルとして，

　

39)　　Ｑ ＝ (ｙ−Ｘβ)'(ｙ−Ｘβ) − 2(Ａβ−ｂ)λ

　

をβについて偏微分して０とすると，

　

40)　 ∂Ｑ／∂β　＝ − 2Ｘ'ｙ ＋ 2Ｘ'Ｘβ − 2Ａ’λ ＝ ０

　

となる．ここで，40)式の解をβ_ｐで表すと，

　

41)　 β_ｐ ＝ (Ｘ’Ｘ)^−１Ｘ'ｙ ＋ (Ｘ'Ｘ)^−１･Ａ'･λ ＝ β_ｏｌｓ＋ (Ｘ’Ｘ)^−１Ａ'λ

　

ただし，41)式において，β_ｏｌｓは制約無しの通常の最小二乗法（ＯＬＳ：Ordinary Least Squares）による推定値である．この41)式の左右両辺に，左側からＡを掛けると，

　

42)　 Ａβ_ｐ ＝ Ａβ_ｏｌｓ＋ Ａ(ＸＸ)^−１･Ａ'λ

　

になる．ここで，38)式より，Ａβ_ｐ ＝ ｂでなければならないから，これを42)式に代入して整理すると，

　

43)　ｂ ＝ Ａβ_ｏｌｓ＋ Ａ(Ｘ'Ｘ)^−１Ａ'λ

　　　Ａ(Ｘ'Ｘ)^−１･Ａ'λ ＝ ｂ − Ａβ_ｏｌｓ

　　　∴　　　λ ＝ ［Ａ(Ｘ'Ｘ)^−１･Ａ'］^−１（ｂ − Ａβ_ｏｌｓ）

　

になる．ここで，ｒ×ｒ行列Ａ(Ｘ'Ｘ)^−１Ａ'は位ｒであるから，非特異行列である．この43)式を41)式に代入すれば，

　

44)　 β_ｐ ＝ β_ｏｌｓ＋ (Ｘ'Ｘ)^−１Ａ'［Ａ(Ｘ'Ｘ)^−１･Ａ'］^−１（ｂ − Ａβ_ｏｌｓ）

　　　　　＝ β_ｏｌｓ− (Ｘ'Ｘ)^−１Ａ'［Ａ(Ｘ'Ｘ)^−１･Ａ'］^−１（Ａβ_ｏｌｓ − ｂ）

　

になる．この44)式が，ＣＬＳ推定値を求める式である．

　

　このガウス−ニュートン法の最大の長所は，計測対象となる式が複数存在しても，残差の均一分散を仮定すれば一度に測れてしまうことである．これは，例えば，二次式効用指標関数において，第２財が「その他一般の財」ではなく第１財同様に実態のある特定の財であり，第１財に関する需要関数

　

8)式再掲)　q_１ ＝｛(a_２２p_１ − a_１２p_２)Ｙ｝／｛a_２２p_１^２−2a_１２p_１p_２＋a_１１p_２^２｝

　　　　　　　　＋ ｛a_２p_１p_２ − a_１p_２^２｝／｛a_２２p_１^２−2a_１２p_１p_２＋a_１１p_２^２｝＋ u_ｔ

　

　

と，第２財に関する需要関数

　

45)　 q2　＝ ｛（−a_１２p_１ ＋ a_１１p_２）Ｙ｝／｛a_２２p_１^２−2a_１２p_１p_２＋a_１１p_２^２｝

　　　　　 ＋｛−a_２p_１^２ ＋ a_１p_１p_２｝／｛a_２２p_１^２−2a_１２p_１p_２＋a_１１p_２^２｝＋ u_ｔ

　

とを共に測らなくてはならず，かつ，8)式と45)式のように測りたい２本の式の双方に同一のパラメタが現れる場合に，一度に測れてしまうので便利な特徴であると考える．ただし，前提条件として，攪乱項（残差）u_ｔの分散が8)式でも45)式でも均一であることが要請される．8)と45)は一見して分かるごとく，二次式効用指標関数に基づく需要関数は，限界効用均等式と収支制約式という２本の構造方程式から導かれた誘導型方程式でありながら，「誘導型パラメタを測定してそれから構造パラメタを逆算する」という通常の間接最小二乗法を用いることが式の形状からできず，誘導型方程式上で構造パラメタを直接推定する必要があり，かつ，構造パラメタが複数本の需要関数に同時に現れるという面倒な特徴を持つ．このようなときに，ガウス−ニュートン法を使うと，複数の誘導型方程式にまたがる同一の構造パラメタを一度に計測できてしまうので，便利であろう．筆者の修士学位論文研究においては，筆者は２費目分割を行って第１財に関する需要関数8)式のみを計測したから，この便利な特徴を活用しなかったが，次回の研究テーマとして「自動車と自動車保険」のような補完関係にある財を計測する場合にはモデルを連立モデル（計量経済学の用語では同時モデル）にする場合，ガウス−ニュートン法によるこの特徴は，極めて大きな効力を発揮するものと考える．当然，この特徴は，２費目分割よりもよりも分割数の大きなｎ費目分割の二次式効用指標関数を想定して，需要関数の具体形が8)式，45)式と異なってくる場合でも，本質的な相違点はなんら発生するものではない．

　

　また，ガウス−ニュートン法の最も大きな利点は，偏微分を行う必要があるのが１階の微分までで済む，という点である．これは，ニュートン法において２階の微分まで求めなければならなかったことと比べると，計算上の操作容易性が向上している点である．筆者が開発したガウス−ニュートン法のプログラムでは，この１階の偏微分を数式の形でプログラムに与える構造を採っており，偏微分の数式の形を具体的に知っていなくても，微分の定義式を用いて近似計算してしまったニュートン法のプログラムの方法は採用しなかった．これにより，新たに計測したい式が増えた場合のプログラム・メンテナンス・ロードは増大してしまうが，計算の正確性は先述のニュートン法よりも良いようであり，本研究では，ニュートン法だと収束しなかったモデルでもガウス−ニュートン法では収束する場合がほとんどであった．これは，恐らくは丸めの誤差が影響したものと推測しているが，この丸めの誤差は演算の精度に起因するものであるので，プログラム実行の処理系や環境によっては状況が変わるであろうから，微分の値の与え方としてどちらが好ましいかという点については，一般論としては何も言えない．

　

　

５．計測に使用する，モデル選定基準統計量

　

　モデル選定基準統計量とは，複数のモデルの内，どちらの方が当てはまり良いかという判定基準を与えてくれる統計量である．本研究では，非線形回帰推定で求められた二次式効用指標関数による実測結果の中から最も当てはまりの良いモデルを選択したり，前回の研究（神山［1989］）でのベルヌイ−ラプラス型による需要関数と今回の二次式による需要関数とを比較したりする際の判定基準を与えることになる．

　

　モデルの当てはまり具合を検証する統計量の代表としては決定係数R^２があるが，R^２には説明変数を増やすとR^２の値が上がってしまうという欠点があるため，（線形回帰が通用す

る範疇の計測でも）複数の異なるモデル間の比較には使えなかった．この欠点を補うもの

として自由度修正済決定係数 R^２があり，線形回帰分析の世界ではこれがよく使われてきた．筆者の修士学位論文研究においては，計測対象モデルが非線形であることに鑑み，タイルのＵ，雨宮のPC基準（Prediction Criterion），赤池の情報量基準AIC（Akaike's Information Criterion）を用いることにした．外挿テストにおいては，だいたいの当てはまり具合の勘所をつかむために標本決定係数も用いたが，モデルの選定基準としてはあくまでもこの３種類の統計量を，�@雨宮のPC基準，�A赤池の情報量基準AIC，�BタイルのＵの優先順位で用いた．ただし，雨宮のPC基準と赤池の情報量基準とでは，Takeshi Amemiya, "Selection of Regressors", International Economic Review, Vol.21, No.2, June, 1980，p.353に'Note that PC 　and Aic are complementary....(continue)'と記されているとおり，両基準でモデルの指標が逆転することがなかったので，両者間で優先順位付けをする必要はなかった．モデルの選定基準として何を選んだかという点も実験計画の一環としては欠かせないと考えたので，この節では，これら３種類の統計量について説明する．ただし，これら統計量の証明や詳細な解説については省略する．

　

　まずは，タイルのＵについて説明する．タイルのＵは，従来から非線形回帰分析の際には決定係数の代わりとなる当てはまり判定基準としてよく用いられてきたものである．タイルのＵは，平均平方誤差ＲＭＳＥ（Root Mean Square Error）

　

　

には，変数の絶対値の大きさに依存して値が代わる欠点があることを修正するべく提唱された統計量であり，

　

　

となるものである．このタイルのＵは，その値が小さければ小さいほど，モデルの当てはまり具合が良いという指標である．

　

　次に，雨宮のPC基準について述べる．雨宮のPC基準とは，

　

3)　 PC＝σ_{ｅｓｔｉｍａｔｅｄ}^２･（１＋m／n）

　

というものである．これは損失関数(y_ｐｉ−y_ｉ)^２を最小化ならしめるべき指標となるもので，mは測定すべき回帰係数の個数，nは標本の大きさ（サンプル・サイズ）である．雨宮のPC基準は，Takeshi Amemiya, "Selection of Regressors", International Economic Review, Vol.21, No.2, June, 1980，pp.337-338,およびpp.352-353によれば非線形モデルにも使用可能なものであり，タイルのＵ同様に，その値が小さい方が望ましいという指標である．

　

　最後に，赤池の情報量基準AICについて述べる．赤池の情報量基準AICはTakeshi Amemiya, "Selection of Regressors", International Economic Review, Vol.21, No.2, June, 1980，pp.352-353によれば，尤度関数さえ定義できればどのモデルでも使えるため，やはり非線形モデルでも適用可能なものであり，雨宮のPC基準の代わりになるものである．赤池の情報量基準AICは

　

4)　 AIC ＝ −2logＬ(θ|y) ＋ 2k

　

というものである．ただし，Ｌ(θ|y)はモデルの最大尤度であり，kはモデルf(y|θ)のパラメタθの次元，すなわちパラメタの個数である．ここで，この4)式において，右辺第１項は最大対数尤度に−２を掛けたものであるから，モデルの適合度(goodness of fit)の高低を測る量である．第２項は，モデルの含む未知母数の個数の２倍であるので，パラメタの増加に対するペナルティーである．したがって，赤池の情報量基準AICも，タイルのＵや雨宮のPC基準と同様に値が小さい程望ましい指標である．ここで，回帰分析において攪乱項（誤差項）u_ｔは多変量正規分布Ｎ(0,σ^２I)にしたがうことを仮定している．この仮定の下での回帰モデルの誤差項の尤度関数は，

　

　

である．このままでは扱いにくので対数尤度をとると，

　

　

ここで，非線形回帰の場合でも線形回帰の場合でも同様に，

　

　

である．7)式から，

　

　

であるから，6)式においてσ^２_{ｅｓｔｉｍａｔｅｄ}≒σ^２とした上で，8)式を6)式に代入すれば，

　

9)　 logＬ ＝ − (n/2)･log(2π) − (n/2)･logσ^２_{ｅｓｔｉｍａｔｅｄ} − n/2

　

となる．これが回帰モデルの誤差項の最大対数尤度である．したがって，これを4)式に代入して，

　

10)　AIC ＝ n･log(2π) ＋ n･logσ^２_{ｅｓｔｉｍａｔｅｄ} ＋ n ＋ 2k

　

これが回帰モデルの当てはまり具合を判定するための赤池の情報量基準AICの具体形になる．

　

　

　

　

　

学術研究の部屋へ　　

　

このホームページのトップへ

　

　

このHP作者に関係のある別HPです。

本稿は，筆者の修士学位論文の中から，非線形最小二乗法（非線形回帰分析；ニュートン法，パターン法，ガウス−ニュートン法）の解法の解説を行っている部分を抜き出して，一般的な非線形最小二乗法の解説としたものである．；keywords = 非線形,最小二乗法,ガウス法,パターン法,ガウス−ニュートン法,回帰分析