曲線あてはめ

提供: miniwiki
移動先:案内検索

曲線あてはめ(きょくせんあてはめ)またはカーブフィッティング: curve fitting[1][2][3][4]は、実験的に得られたデータまたは制約条件に最もよく当てはまるような曲線を求めること。最良あてはめ、曲線回帰とも。一般に内挿回帰分析を用いる。場合によっては外挿も用いる。回帰分析で曲線を求める場合、その曲線はデータ点を必ず通るわけではなく、曲線とデータ点群の距離が最小になるようにする。曲線あてはめによって得られた曲線を、近似曲線という。特に回帰分析を用いた場合には回帰曲線という。現実の実験データは直線的ではないことが多いため散布図、近似曲線を求める必要性は高い。

一般論

最小二乗法による最適関数の推定

我々が考えるべき問題は、実験データを実験を説明する「説明変数」と「目的変数」に分類した上で、説明変数[math]\textbf{x}[/math] と、目的変数yの関係

[math]y=g(\textbf{x})[/math]

を求めることである。説明変数としては測定条件を考えることが多く、目的変数としては、測定値を考えることが多い。説明変数、目的変数共にベクトル量である可能性があるが、測定値のほうは、多変数関数の微分が、値域側の成分に関して独立であることからスカラー量としても一般性を失わない。一方、多変数関数の微分は、定義域側の成分については独立でないため、一般論を述べる上ではベクトル量としておかなければならない。以下、測定条件は、k次元ベクトルの形で与えられているとする。成分で表記すると

[math]\textbf{x}=\left( \begin{matrix} x_1\\ \vdots\\ x_k\\ \end{matrix} \right)[/math]

となる。

実験データは、説明変数に関するデータ[math]{\textbf{x}}_1 ,\cdots ,\textbf{x}_n[/math] と目的変数に関するデータ[math]y_1 ,\cdots ,y_n[/math]の組、 [math](y_1 ,\textbf{x}_1 ),\cdots ,(y_n ,\textbf{x}_n)[/math] の形で得られる。また、j番目の測定条件 [math]\textbf{x}_i[/math] の第i成分を [math]x_{ij}[/math] で表すものとする。

我々が考えるべき問題は、適当な [math]\alpha[/math] 個のパラメータ [math]a_1, \cdots ,a_{\alpha}[/math] と、k+ [math]\alpha[/math] 変数の関数 [math]f(\textbf{x}, \textbf{a} )[/math] を考え、 [math]\textbf{a}[/math] の値を調整し

[math]S(\textbf{a} )=\sum_{j=1}^n |y_i -f(\textbf{x}_j ,\textbf{a} )|^2[/math] ・・・・ (1-1)

を最小とするような、[math]\textbf{a}[/math] を求める問題に帰着される[1]。このSの平方根 [math]\sqrt{S}[/math] のことを、「関数当てはめ時の誤差」という。ここで

[math]{\textbf{a}}=\left( \begin{matrix} a_1\\ \vdots\\ a_{\alpha}\\ \end{matrix} \right)[/math]

のことを、フィッティングパラメータと言う。また、関数Sを考えるときには [math](y_1 ,\textbf{x}_1 ),\cdots ,(y_n ,\textbf{x}_n )[/math] は、もはや定数ベクトルでしかないことに注意されたい。飽くまで関数Sの変数は [math]a_1,\cdots a_{\alpha}[/math] である。

なお、Sを定義するにあたり、各データに対して、適当な定数(正または0)[math]w_1 ,\cdots , w_n[/math] によって重みを付け、

[math]S(\textbf{a} )= \sum_{j=1}^n w_j | y_i -f( \textbf{x}_j , \textbf{a} )|^2[/math] ・・・・ (1-1')

のようにすることもある。この方法によって、y方向に誤差(Yエラーバー)がある場合や、「測定回数の異なるデータの平均」の比較が可能であるが、x方向にも誤差(Xエラーバー)がある場合には、対応できない。x方向にも誤差がある場合には、デミングの方法[1]を用いる。なお、(1-1)は「(1-1')において、全てのデータの重みづけが等しい状況」を意味することに注意されたい。

我々が考えるべき問題は、、(1-1)あるいは(1-1')の関数Sの極値問題[5]に他ならない。一般に、極値問題は解を持たない可能性があり、また、解が存在したとして、重解の可能性もあるが、一般論として、以下の定理が知られている。

「もしも、[math]\textbf{a}_0[/math] で、Sが極値をとるとすると、[math]{\rm grad} S( \textbf{a}_0 )=\textbf{0}[/math] である。」

この定理は、最適なフィッティングパラメータに対する必要条件を与える。極小値を与えるような [math]\textbf{a}[/math] の十分条件としては、

「Sの [math]\textbf{a}[/math] におけるヘッセ行列正定値(正値, 正定符号)となること」[5]

がある。極小値が仮に存在したとして、それらが必ずしも最小であるとは限らない。例えば、最適な [math]{\textbf{a}}[/math] が無限遠に存在する可能性もある。

(1-1)のSを最小とするようなフッティングパラメータ [math]\textbf{a}_0[/math] が得られた場合には、以下の [math]g(\textbf{x})[/math] を、最適関数(xがスカラーの場合には、最適曲線)という。

[math]g( \textbf{x} )=f( \textbf{x} , \textbf{a}_0 )[/math] (1-2)

このgは、説明変数 [math]\textbf{x}[/math] と目的変数yの間に1つの関数関係を与えている。つまり、このgは、 [math]\textbf{x}[/math] とyの関数であり、フィッティングパラメータは定数ベクトルと考える。

一般には、「必ず [math](y,\textbf{x} )=(0,\textbf{0} )[/math] を通る」といった付帯条件が付いている場合がある。このような場合には、ラグランジュの未定乗数法[5]が最適なフィッティングパラメータを探る上で手掛かりを与える。

なお、付帯条件のある場合、ない場合共に、実際の数値計算では、Levenberg-Marquardt algorithm(レーベンバーグ・マルカート法)が用いられることが多い。

様々な曲線あてはめ

直線または多項式曲線

ファイル:Curve fitting.png
正弦関数のデータ点群(黒)に対して、1次多項式(赤)、2次多項式(緑)、3次多項式(黄)、4次多項式(青)をあてはめた図

まず、次のような1次多項式を考える。

[math]y = ax + b\;[/math]

これは、傾斜 a の直線である。一般に、このような直線はx座標の異なる2点によって一意に定まる。したがって1次多項式は、x座標の異なるデータ点がちょうど2個ある場合に、正確にそれらを通る直線となる。最小二乗法を使う場合には、データ点が何個あっても、最適な直線が一意に定まる。ただし、最適な直線とは、残差平方和が最小というだけで、そのデータの素性を最もよく表しているとは限らない。

次数を上げて2次多項式にすると、次のような式になる。

[math]y = ax^2 + bx + c\;[/math]

この場合、x座標の異なる任意の3点に当てはめることができる。

さらに次数を上げて3次多項式にすると、次のような式になる。

[math]y = ax^3 + bx^2 + cx + d\;[/math]

この場合、x座標の異なる任意の4点に当てはめることができる。

より一般化すれば、4つの「制約」を正確に満足する、と言える。制約は点だけでなく、角度曲率(接する円の半径の逆数)などもある。角度や曲率の制約は曲線の端に設定することが多く、それを端末条件 (end condition) と呼ぶ。多項式曲線を連結したスプライン曲線が滑らかな曲線となるには、連結する両方の多項式曲線で端末条件を同一にする必要がある。より高次の制約として、例えば「曲率の変化率」といった制約を与えることもある。これは例えば、クローバー型インターチェンジで通行する自動車にかかる力を決め、制限速度を決定するのに役立つ。

また、1次多項式は1つの点と角度の制約に当てはめることができ、3次多項式は2つの点と1つの角度、1つの曲率の制約に当てはめることができる。他にも様々な制約の組み合わせで各種次数の多項式を当てはめることができる。

n + 1 個より多い制約があるときでも(n は多項式の次数)、それらを満足する多項式曲線を描くことができる場合がある。例えば3つの点が同一直線上に並ぶような配置であれば、1次多項式を正確にあてはめることができる。しかし、このような配置は例外的であって稀である。通常は全制約を正確に満足することは期待できない。よって、一般には近似度を評価する方法を必要とすることになる。最小二乗法は最も一般的な方法である。

ここで、なぜ近似ではなく多項式の次数を高くして正確に当てはめようとしないのかという疑問が生じる。それには、以下のような理由がある。

  • 正確な一致が存在するとしても、それを計算できるとは限らない。使用しているアルゴリズムによっては計算が発散してしまって解を求められなかったり、非常に時間がかかることがある。
  • データそのものに誤差がある場合、各点を正確に通る曲線よりも近似的な曲線の方が好ましい場合がある。
  • ルンゲ現象が起きやすい。n次多項式曲線の変曲点の数は最大 n-2 である。したがって、一般に次数が低いほうが曲線はより滑らかになる。次数が高くても滑らかな曲線にすることは可能だが、次数が低い方が簡単である。

ここまで、多項式の次数が制約数より少ない場合を述べてきたが、逆に多項式の次数が制約数より大きい場合はどうなるだろうか。上述の高次多項式の問題が全て生じることになるだけでなく、解が一意に定まらないという問題も生じる。そこから1つの曲線を選択するのはソフトウェアや人間の役割となる。このため、近似でよい場合は次数をなるべく低く設定するのが一般に最善とされている。


その他の曲線

場合によっては、円錐曲線(円、楕円、放物線、双曲線など)や三角関数(サイン、コサインなど)の曲線を使うこともある。例えば、空気抵抗を無視すると、重力の影響下にある物体の軌跡は放物線を描く。したがって、そのような物体の観測結果に放物線を当てはめることは意味がある。潮の干満は正弦波のパターンを描くので、干満に関するデータには正弦曲線を当てはめることができる。あるいは、より正確には月と太陽の影響を考慮して、2種類の正弦曲線の合成を当てはめることもある。

代数的曲線あてはめと幾何学的曲線あてはめ

代数的なデータ解析においては、曲線とデータ点のY軸方向の距離を最小化するような曲線を求める。しかし、グラフィックスや画像の分野では、幾何学的な曲線あてはめ、すなわち曲線とデータ点の直交する距離が最小になるような曲線を求める(あるいは、X軸とY軸の両方で曲線とデータ点の距離を最小化する)方が見た目がよい。しかし、幾何学的曲線あてはめは非線形かつ反復的な計算を必要とするため、あまり使われない。

円の幾何学的あてはめ

Coope[6] は、2次元のデータ点群に視覚的に最もうまくあてはまる円を求める問題を考えた。その手法は非線形な問題を数値解析的な反復なしに線形な問題にみごとに変換でき、従来の技法よりも劇的に性能を向上させることができる。

楕円の幾何学的あてはめ

上述の技法は非線形なステップを1つ追加することで楕円に一般化でき[7]、高速に様々な向きと離心率の楕円を見付けることができる。

曲面への応用

以上、2次元の曲線を論じてきたが、同じことは3次元空間における曲面にも適用できる。通常、uv という2つの方向を表すパラメータに沿った曲線の網で定義されるパッチ(あて布)を並べたものが曲面となる。

ソフトウェア

R言語のような統計パッケージ、GNU Scientific LibrarySciPy のような数値解析ソフトウェアには、曲線あてはめ用コマンドがある。また、曲線あてはめ専用のプログラムとして、TableCurve、Fityk などがある。詳しくは外部リンクを参照。

脚注

  1. 1.0 1.1 1.2 本間 仁,春日屋 伸昌「次元解析・最小二乗法と実験式」コロナ社(1989)
  2. 加川 幸雄,霜山 竜一「入門数値解析」朝倉書店(2000)
  3. John R. Taylor、林 茂雄、 馬場 凉「計測における誤差解析入門 」東京化学同人(2000)
  4. 吉沢 康和「新しい誤差論―実験データ解析法 」共立出版 (1989/10)
  5. 5.0 5.1 5.2 島 和久「多変数の微分積分学」近代科学社 (1991/09)
  6. Coope, I.D., Circle fitting by linear and nonlinear least squares, Journal of Optimization Theory and Applications Volume 76, Issue 2, New York: Plenum Press, February 1993
  7. Paul Sheer, A software assistant for manual stereo photometrology, M.Sc. thesis, 1997

外部リンク

実装

オンラインの解説

オンライン計算器、アプリケーション、デモ