楕円近似について

　例えば，惑星は太陽をひとつの焦点とする楕円軌道をとるが，そのようなデータに対して楕円軌道をあてはめたいという要望は多いだろう．そこで，今回のコラムでは，楕円近似について，３通りの方法を紹介することにした．

　このサンプルデータでは，楕円の全周にわたって２５個の２元データ（ｘ，ｙ）が得られているが，場合によっては，弧の一部だけに基づいて楕円を確定することが要請されることもあるだろう．今から２００年前，ガウスはそのようなデータに最小２乗法を適用して小惑星セレスの発見を導いた．最小２乗法は２００年にわたる英知の結晶ともいえるのだが，まず，この逸話から紹介することにしたい．

　今日，小惑星は火星と木星のあいだの空間に幾百となく散在していることが知られています．天文学史上，最初に発見された小惑星はセレスで，セレスは第１号であり，かつ，最大の小惑星です．セレスは１８０１年１月１日の夜，イタリアの天文学者ピアッツィによって発見されました．ピアッツィはボーデの予測した位置にある小惑星の運行を追い続けましたが，２月初旬，この天体は太陽に接近しそのまぶしい光に消されてしまったために観測できなくなってしまいました．

　ピアッツィが観測した４１日間だけのデータ（９°の弧）だけを使って楕円軌道を確定することは，当時の天文学者たちの計算能力の限界を超えていました．なぜなら，それまでの軌道決定法は豊かな資料に基づくものであり，セレスの場合，少ないデータからケプラー運動を推論することが要請されたからです．そこで，２４才の若きガウスはたった３回の完全な観測からその軌道を計算し，太陽の近くで姿を消してしまったセレスがその年の終わり頃再び姿を現す位置を計算しました．１８０１年１２月３１日，セレスはガウスの予測した位置に再び姿を現しました．

　その位置は粗い円軌道近似で推定したものよりも７°以上も東にずれていましたから，結局，ガウスの予測は非常に正確であることがわかり，若いガウスに最初の大きな名声を与えることになりました．この成果は予知と観測とニュートン力学による軌道計算の劇的な合流点を表す天文学史上の事件であったと考えられるのですが，いまでもガウスの最も知られた業績の一つになっています．セレス以降，ドイツでは小惑星の発見ラッシュとなり，１８０２年に２番目の小惑星パラス，１８０４年に３番目の小惑星ジュノー，１８０７年に４番目の小惑星ベスタが発見されています．

　１７９４年，ガウスは１８才のときすでに最小２乗法を考案していたと記録されていますが，ガウスによって天体の運動・軌道を決定するための新しい方法として創始され，ある時刻の位置を予測して再発見の手がかりを与えた方法が最小２乗法なのです．当時の望遠鏡の解像度を考えると，位置の予測なしに再発見は難しかったと思われますが，小惑星セレスの再発見によって最小２乗法は有名になり，実用に供されるようになりました．そして，１８２１年と１８２３年にガウスは最小２乗法を発表し，１８２０年代までに今日広く使われている最小２乗法の基本的な大筋が完成しています．

　また，ガウスはこの過程で実験データの期待値からのバラツキが一定の法則に従うことに注目し，その分布を理論的に計算しました．この分布が正規分布（ガウス分布）で，ガウスの結果を厳密に証明したものが中心極限定理です．正規分布はあらゆる種類のデータ解析において中心的な役割を果たしています．測定誤差の基礎となる誤差論も１９世紀の始めガウスによって始められ，これが端緒となって数理統計学が進歩しました．換言すれば，数理統計学は正規分布を中心として展開され，ガウス以来２００年，誤差，変動，撹乱，ばらつき，偏りのあるデータを適切に処理し情報を抽出する方法を開発してきたのです．

【補】小惑星（アステロイド）の中には，セレス，パラス，ジュノーのように直径が数百キロメートルのものもありますが，直径が１キロにも満たないものもあります．天体観測の精度があがるにつれて，ますます小さい天体が発見されていますが，どこまでが小惑星で，また，どこからが宇宙塵になるか区別することは天文学者にとっては興味のあることでしょう．しかし，私が調べた限り，その違いの明確な取り決めはなされていないようです．

　前述のサンプルデータには測定値ｚがないので，ｚi＝０を加える，すなわち（ｘi，ｙi，０）のように，測定値ｚiをダミー変数のようにして用いることにする．

とすることになる．通常の最小２乗法プログラムを多変数問題を解くことができるように書き換える必要がでてくるが，計算原理は１元でも２元でも独立変数が多変数の場合でもほぼ同じである．

　　Ｆ（ｚ）＝ａ1Φ1（ｘ，ｙ）＋ａ2Φ2（ｘ，ｙ）＋・・・＋ａ5Φ5（ｘ，の形に書くことができる．これに対して，多変量線形最小２乗法によるあてはめを行った．

　線形最小２乗法では計算は１秒以内で完了したが，後述する非線形最小２乗法に比べて，あてはめ精度がよくないのが欠点である．

は結構面倒であった．自分で式を立ててみれば，この問題の面倒さがわかるので，係数変換式については読者の演習問題としよう．

　　ｆ（ｘ，ｙ）＝｛（ｘ－ｘ0）ｃｏｓθ＋（ｙ－ｙ0）ｓｉｎθ｝^2／ａ^2＋｛（ｘ－ｘ0）ｓｉｎθ－（ｙ－ｙ0）ｃｏｓθ｝^2／ｂ^2－１

　この式では，未定係数が三角関数の中に含まれているため，どのように変形してもパラメータに関して線形とはならない．すなわち，本質的に非線形なので，線形最小２乗法は使えないのである．

　また，非線形最小２乗法では，５個のパラメータａ，ｂ，θ，ｘ0，ｙ0の初期値が必要になる．ここで用いるサンプルデータのように楕円の全周にわたって２元データ（ｘ，ｙ）が得られている場合は，かなり正確な初期値をプログラム上で設定することができるが，弧の一部のデータしかない場合もあろう．

　　ｆ（ｘ，ｙ）＝｛（ｘ－ｘ0）ｃｏｓθ＋（ｙ－ｙ0）ｓｉｎθ｝^2／ａ^2＋｛（ｘ－ｘ0）ｓｉｎθ－（ｙ－ｙ0）ｃｏｓθ｝^2／ｂ^2－１

　多変量の非線形最小２乗法を用いる方法では，計算が収束するまで数秒を要したが，あてはまりの良さに関しては満足のゆくものであった．

　また，この非線形モデルのあてはめに対しては，ガウス・ニュートン法を用いたが，初期値と計算アルゴリズムの選び方によっては収束しないこともあり得ることを注意しておく．

　【２】【３】項で述べた多変量最小２乗法は，最も一般的に用いられている楕円近似法であると思われるが，「ｚ軸のデータだけに誤差がある場合で，ｘ，ｙ軸のほうの誤差は全くないかあるいは無視できる．」という前提をおいて，目的関数をｚ軸方向の残差２乗和

　実用上はこれで十分な場合も多いが，ｘ，ｙ軸側のデータといえども必ずしも理想的には管理されず，誤差の混入が常であるから，どうしてもｘ，ｙ双方に誤差を考えないといけない場合がある．

　そこで，この節では，ｘ，ｙ軸側のデータにも誤差がある場合のモデル式のあてはめについて，青山学院大学・理工学部・物理を定年退官された粟屋隆先生の方法を用いてあてはめを行う．

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

　粟屋の関数フィッティングの方法は，最尤法に基づいていて，データ（ｘ，ｙ）の真の値を（ｘ0，ｙ0）とすると，目的関数ｓは

と書き表すことができる．これはｚ軸方向の残差２乗和ではなくて，ｘ－ｙ平面におけるデータ点と真値のマハラノビス距離の２乗和である．

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

　　ｆ（ｘ，ｙ）＝｛（ｘ－ｘ0）ｃｏｓθ＋（ｙ－ｙ0）ｓｉｎθ｝^2／ａ^2＋｛（ｘ－ｘ0）ｓｉｎθ－（ｙ－ｙ0）ｃｏｓθ｝^2／ｂ^2－１

をあてはめたが，陰関数ｆ（ｘ，ｙ）＝０のあてはめでは，独立変数，従属変数の区別がなく対等と考えられるから，ｘ側，ｙ側どちらにも誤差がある場合のあてはめ法として，粟屋の２次元フィッティング法を適用することができる．しかし，この計算例の場合，数分の計算所要時間を要した．

　あてはまりの良さをとるか，あてはめ精度を多少犠牲にしても計算速度を重視するかによって，アルゴリズムの優劣のポイントは変わってくる．まとめとして，楕円近似をする際の各方法の性能比較したものを掲げておこう．

　計算原理からいえば，粟屋の方法が優れていることは明らかであるが，計算に時間がかかるという欠点がある．また，適当な初期値を出発点として反復計算する非線形法と粟屋の方法では，ときに計算が収束せず，発散してしまう場合があるので留意されたい．

　その点，多変量の線形最小２乗法では初期値は不要で，収束安定性を考慮する必要もない．プログラムも簡単なため，計算時間が短くて済む．また，楕円の場合は

の係数変換式が利用できるので描画も容易である．ただし，線形法であるから，楕円の基本情報（半径や中心など）の誤差値を求めることは簡単ではない．

　楕円近似をコンピュータ上にインプリメントしようとする際は，要求されているあてはめ精度やどれくらいの計算所要時間が許容されるのか，楕円の基本情報の誤差を求められているかどうかにもよるが，総合的にみて多変量の線形最小２乗法が最もすぐれているものと思われた．