通常用いられる最小２乗法は，計算を単純にするため「縦軸のデータだけに誤差がある場合で，横軸のほうの誤差は全くないかあるいは無視できるものと見なすことができる」という前提をおいて，目的関数を

　　２次元：ｓ＝Σｗi（ｙi−ｆ（ｘi））^2　　　（ｗi：重み）

　　３次元：ｓ＝Σｗi（ｚi−ｆ（ｘi，ｙi））^2 　　　（ｗi：重み）

などとし，誤差はあくまでもｙのみ（３次元の場合はｚのみ）に含まれるものとしています．それとは逆に，ｘのみに誤差がある場合は従属変数と独立変数を入れ替えるだけで，この解析法がそのまま適用できます．

　実用上はこれで十分な場合も多いのですが，横軸側のデータといえども必ずしも理想的には管理されず，誤差の混入が常ですから，どうしてもｘ，ｙ双方に誤差を考えないといけない場合があります．たとえば，ある物質と受容体の結合親和性や受容体量を求めるためのラングミュア・プロット（Langmuir plot）では，縦軸・横軸とも測定濃度をプロットするので，ｘ軸側，ｙ軸側ともに誤差があり，その値が大きいほどその誤差の絶対値は大きくなります．また，最高血圧（東京の最高気温）をｘ，最小血圧（仙台の最高気温）をｙとした散布図ではｘがｙの原因である（ｙがｘの原因である）という関係が想定しにくく，通常の回帰の観念を用いるのは適当ではありません．

　ｘ側，ｙ側どちらにも誤差がある場合は独立変数，従属変数の区別がなく対等と考えられますから，ｆ（ｘ，ｙ）＝０のように陰関数表現しておく方が楕円のような多価関数となるデータに対しても曲線のあてはめができるので何かと都合がよくなります．そのような方法に，青山学院大学理工学部・物理を定年退官された粟屋隆先生による２次元フィッティングの方法があります．

　読者が粟屋の方法に興味を持ち使って下さることを願って，ここで宣伝しておきますが，詳細は参考文献

　　粟屋　隆「データ解析」，学会出版センター

をご参照願います．

　ところが，先日，東海大学・海洋学部の大西修平先生から粟屋の方法に対するひとつの問題が提起されました．粟屋の計算式には母分散σxi^2，σyi^2が含まれるのですが，これらが既知の場合はいいとして未知の場合はどうしたらいいかというものです．実際問題として未知であるほうが当り前ですから，今回のコラムでは母分散σxi^2，σyi^2が未知のときの粟屋の２次元フィッティング法の取り扱いについて考えてみます．

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

【１】ラグランジュの未定係数法

　まず最初に粟屋の方法について説明します．粟屋の方法ではラグランジュの未定係数法が使われていますが，ラグランジュの未定係数法は，たとえば，「ｘ^3−３ｘｙ＋ｙ^3＝０の条件のもとでｘ^2＋ｙ^2の極値を求めよ」といった条件つき極値問題や制約条件付きの最小２乗法などの解を得るために導入された方法です．

　制約条件付きの最小２乗法は１９世紀の測量士たちのおこなった三角測量でもっともよく用いられましたが，以下にその簡単な例を示します．

　「三角形の三つの内角を等精度で測定し，α＝５４°０５′，β＝５０°０１′，γ＝７６°０６′を得た．内角の和は２直角になるべきであるが，測定誤差のため１８０°１２になった．」

　このような場合に内角の最確値ｘ，ｙ，ｚおよび確率誤差を求めるにはどうしたらよいかを考えてみましょう．ｘ＋ｙ＋ｚ＝πが要請されている条件です．また，測定精度は等しいので荷重をｗi ＝１とします．

　このような制約条件付きの問題では，変数λを新たに導入して，目的関数を

　　ｓ＝Σ（測定誤差）^2−λ（条件式）

　　　＝（ｘ−α）^2＋（ｙ−β）^2＋（ｚ−γ）^2−λ（ｘ＋ｙ＋ｚ−π）

としてｓを最小にするｘ，ｙ，ｚおよびλを求めます．極値の必要条件により，　　∂ｓ／∂ｘ＝２（ｘ−α）−λ＝０・・・・・（１）

　　∂ｓ／∂ｙ＝２（ｙ−β）−λ＝０・・・・・（２）

　　∂ｓ／∂ｚ＝２（ｚ−γ）−λ＝０・・・・・（３）

　　∂ｓ／∂λ＝−（ｘ＋ｙ＋ｚ−π）＝０・・・（４）

前の３つの式を加え合わせて，

　　２（ｘ＋ｙ＋ｚ−α−β−γ）−３λ＝０

が得られます．（４）式より，ｘ＋ｙ＋ｚ＝πを代入すると

　　λ＝２／３（π−α−β−γ）

が得られます．このλを（１），（２），（３）式に代入して

　　ｘ＝α＋１／３（π−α−β−γ）＝５４°０１′

　　ｙ＝β＋１／３（π−α−β−γ）＝４９°５７′

　　ｚ＝γ＋１／３（π−α−β−γ）＝７６°０２′

が求める解となります．すなわち，等精度で三角形の内角を測定したときの各補正量は測定値の和と１８０°との差を３等分したものです．

　このように，制約条件のある最小２乗問題を解くには未定係数λを導入して，未知のパラメータと未定係数λについて方程式を解けばよいことになり，制約条件のない場合の手法を用いて容易に解くことができることが理解されるでしょう．

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

【２】粟屋の方法

　粟屋の方法は，品質管理・計量管理で有名なデミング（Deming）のクラシカルな方法を改良したものですが，デミングの方法が反復計算にそぐわない１回限りの計算法であるのに対し，粟屋の方法は反復計算によって精度を逐次高めていくことが可能になっています．

　この２次元フィッティングは実に念の入った膨大なものであり，未知数が多くなって計算はかなり複雑になります．この解法については

　　粟屋隆著「データ解析」（学会出版センター）

に詳しいのですが，特殊な問題にはこのような考察が必要になり，拡張した取り扱いを然るべく行なえば解が得られることをアウトラインだけでも知っておいたほうがよいと思われます．

　詳細は省略し理論の要点のみを記しますが，粟屋の関数フィッティングの方法は，最尤法に基づいていて，データ（ｘ，ｙ）の真の値を（ｘ0，ｙ0）とすると，目的関数ｓは

　　ｓ＝Σ｛（ｘ−ｘ0）^2／σx^2＋（ｙ−ｙ0）^2／σy^2｝

と書き表すことができます．これはｚ軸方向の残差２乗和ではなくて，ｘ−ｙ平面におけるデータ点と真値のマハラノビス距離の２乗和です．もし，σx＝σyならば，

　　ｓ＝Σ｛（ｘ−ｘ0）^2＋（ｙ−ｙ0）^2｝

ですから，ｘ−ｙ平面上のユークリッド距離の２乗和を最小化することと等価になります．すなわち，粟屋の方法は陽関数におけるガウス・ニュートン法を陰関数まで取り扱いを拡張したもので，線形関数・非線形関数いずれの場合もユニバーサルに取り扱うことができるのです．

−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−

　２次元フィッティングでは，取り扱う陰関数をパラメータも含めｆ（ｘ，ｙ，ａ）＝０と書くことにします．ここで，データ（ｘi，ｙi）の真の値を（ｘi0，ｙi0），真の係数を示すベクトルａ0＝（ａ10，・・・，ａm0）とすると，

　　ｆ（ｘi0，ｙi0，ａ0）＝０

が恒等的に成り立ちます．したがって，制約条件式にはｆ（ｘi0，ｙi0，ａ0）＝０を用いることになります．

　目的関数ｓは，

　　ｓ＝Σ（ｘi−ｘi0）^2／２σxi^2＋Σ（ｙi−ｙi0）^2／２σyi^2

　　　＝ΣΣ（ｘi−ｘi0）^2／２σxi^2＋Σ（ｙi−ｙi0）^2／２σyi^2−Σλiｆ（ｘi0，ｙi0，ａ0）　　　　　　（λi ：ラグランジュの未定係数）

と表わすことができますから，

　　∂ｓ／∂ｘi0＝−（ｘi−ｘi0）／σxi^2−λi∂ｆ（ｘi0，ｙi0，ａ0）／∂ｘi0＝０

　　∂ｓ／∂ｙi0＝−（ｙi−ｙi0）／σyi^2−λi∂ｆ（ｘi0，ｙi0，ａ0）／∂ｙi0＝０

　　∂ｓ／∂ａj0＝−Σλi∂ｆ（ｘi0，ｙi0，ａ0）／∂ａj0＝０

　　　（ｉ＝１〜ｎ，ｊ＝１〜ｍ）

を解けばよいことになります．

　ここで解くべき連立方程式の未知数は，ｘi0，ｙi0，λiがｎ個ずつ，ａj0がｍ個で合計３ｎ＋ｍ個となり，未知数が多すぎて計算は大変複雑となります．そこで，まずλiを消去することを考えます．以下，表記を簡単にするため，

　　∂ｆ（ｘi0，ｙi0，ａ0）／∂ｘi0≡〜　ｆxi

　　∂ｆ（ｘi0，ｙi0，ａ0）／∂ｙi0≡〜　ｆyi

　　∂ｆ（ｘi0，ｙi0，ａ0）／∂ａj0≡〜　ｆji

と表わすことにします．

　アプロ−チの途中を省略しますが，誤差伝播の法則より，ｚ＝ｆ（ｘ，ｙ）の母分散σz^2は，分布形にかかわらず，

　　σz^2＝（∂ｆ／∂ｘ）^2σx^2＋（∂ｆ／∂ｙ）^2σy^2

で近似することができることより，重み係数に相当するのが

　　ｗi＝｛（σxiｆxi）^2＋（σyiｆyi）^2｝^(-1)

また，２変数関数のテイラー展開より

　　ｆ（ｘ＋Δｘ，ｙ＋Δｙ）＝ｆ（ｘ，ｙ）＋Δｘｆx＋Δｙｆy＋１／２（Δｘ^2ｆxx＋２ΔｘΔｙｆxy＋Δｙ^2ｆyy）＋・・・・・

　ニュートン法の漸化式ではｆ（ｘ＋Δｘ，ｙ＋Δｙ）＝０とおいてΔｘ，Δｙを求めるわけですから，残差の近似値に相当するのがｆi0＋（ｘi−ｘi0）ｆxi＋（ｙi−ｙi0）ｆyiになり，結局，

　　Ｆj≡Σｗi｛ｆi0＋（ｘi−ｘi0）ｆxi＋（ｙi−ｙi0）ｆyi｝ｆji＝０　　　（ｊ＝１〜ｍ）

を解けばよいことになります．

　これをガウス・ニュートン法と同様の方法で整理することにします．真の値ｘi0，ｙi0，ａ10，・・・，ａm0の近くの点（ｘi，ｙi，ａ1，・・・，ａm）で関数ｆ（ｘi0，ｙi0，ａ0）＝０をでテイラー展開し，１次までの項をとると

　　ｆi0＝ｆi＋ｆxiΔｘi＋ｆyiΔｙi＋ΣｆkiΔａk＝０

となり，また，微分係数は点（ｘi，ｙi，ａ1，・・・，ａm）における値で置き換えると，

　　Σｗi｛ｆi＋（ｘi−ｘi0）ｆxi＋（ｙi−ｙi0）ｆyi｝ｆji＝ΣｗiｆaiΣｆkiΔａk

が得られ，Δａkが計算されることになるのです．

　多価関数に対する取り扱い（２次元フィッティング）において，パラメータを計算するには最初に未知のパラメータの初期値を与えなくてはいけないことはガウス・ニュートン法と同じですが，粟屋の方法ではさらに測定値の真の値も与えなければなりません．測定値の真の値の初期値には測定値自身を用いることにし，パラメータと測定値の真の値の第１近似値を与えてやって逐次近似解法により漸近連立方程式を解きます．

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

【３】大西先生の問題提起

　前節で示したように「粟屋の方法」には母分散σxi^2，σyi^2が含まれています．必ずしもそれらが母分散といった量である必要はありませんが，それらの値があらかじめ与えられていないと計算できません．したがって，これらが既知の場合はいいとして未知の場合はどうしたらいいかというのが大西先生からの問題提起でした．

　大西先生のご指摘の通り，確かに一般的なデータ解析では分散も未知として扱う必要があります．粟屋の方法は最尤法に基づいていますから，分散が未知の場合についても対数尤度を書き下し最大化を実行する必要があるのですが，データ（ｘi，ｙi）が独立な誤差をもつときの確率分布は，２次元正規分布：

　　ｐ（ｘi，ｙi）＝１／２πσxi^2σyi^2ｅｘｐ｛−（ｘi−ｘi0）^2／２σxi^2−Σ（ｙi−ｙi0）^2／２σyi^2｝

ですから，対数尤度関数は，

　　Ｌ＝Σ｛−（ｘi−ｘi0）^2／２σxi^2−ｌｏｇσxi−Σ（ｙi−ｙi0）^2／２σyi^2−ｌｏｇσyi｝−ｎｌｏｇ２π

と書けます．

　これをσxi，σyiで偏微分し０とおくことによって

　　∂Ｌ／∂σxi＝０　→　σxi^2＝（ｘi−ｘi0）^2

　　∂Ｌ／∂σyi＝０　→　σyi^2＝（ｙi−ｙi0）^2

が得られます．これを

　　ｗi＝｛（σxiｆxi）^2＋（σyiｆyi）^2｝^(-1)

に代入すると，σxi^2，σyi^2が陽には現れない形になって大西先生の提起された問題点が解消されます．

　実際に計算してみたところ，誤差の小さいデータに大きな重みのかかった形の解が得られました．また，計算所要時間は余分にかかりましたが，それは重みｗiのなかに未知のパラメータｘi0，ｙi0が潜在変数として含まれるために応答曲面が滑らかにならないためと考えられました．なお，繰り返し計算の初回にはｗi＝∞となるのを防ぐ工夫が必要でした．

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

【４】f(x,y,θ)=0型回帰の信頼区間

　以下，∂f/∂θk，∂f/∂x，∂f/∂yをそれぞれfk,fx,fyと略記することにしますが，この場合も，未知母数の信頼区間を与えてくれる分散共分散行列{cov(θi,θj)}は正規方程式の係数行列の逆行列として求めることができますし，関数ｆの誤差は，パラメータの回りでテイラー展開して，

　　(Δf)2=ΣΣfifjcov(θi,θj)

で与えられることも，y=f(x,θ)型回帰の場合と同様です．

　横軸側にも誤差がある場合の信頼区間表示で異なっている点は，この誤差を縦軸方向と横軸方向に分解する必要があることです．再び，誤差伝播の法則

　　(Δf)^2=(fxΔx+fyΔy)^2

を用いますが，この式は，いわば誤差の分散公式

　　σ2(x+y)=σ2(x)+σ2(y)+2rσ(x)σ(y)

であって，z=x+yとするとｚ軸方向の合成分散をｘ軸，ｙ軸の２つの方向に分配すると考えることができます．

　その際，ベクトル（Δｘ，Δｙ）は法線方向を向きますが，陰関数ｆ（ｘ，ｙ）＝０上の点（ｘ，ｙ）で接線の方程式を求めるには，２変数関数の微分の知識が必要で，ｆ（ｘ，ｙ）＝０のｙをｘの関数（ｆ（ｘ，ｙ（ｘ））＝０）とみなして，両辺をｘで偏微分すれば２変数関数の合成微分の公式によって

　　∂ｆ／∂ｘｄｘ／ｄｘ＋∂ｆ／∂ｙｄｙ／ｄｘ＝０

すなわち，ｆx＋ｆyｄｙ／ｄｘ＝０より，

　　ｙ’＝ｄｙ／ｄｘ＝−ｆx／ｆy

が得られます．

　これが陰関数定理ですが，これにより，

　　Δｙ／Δｘ＝ｆy／ｆx

となり，独立な誤差をもつときは

　　(Δf)^2=(fxΔx)^2+(fyΔy)^2

ですから

　　（Δｘ）^2＝（Δｆ）^2ｆx^2／（ｆx^4＋ｆy^4）

　　（Δｙ）^2＝（Δｆ）^2ｆy^2／（ｆx^4＋ｆy^4）

が得られます．

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝