Pythonでフィッティング (線形/非線形/多変数を例示し、誤差評価まで)

2020 7/05
Pythonでフィッティング (線形/非線形/多変数を例示し、誤差評価まで)

今回はpythonを使ったフィッティング(fitting)のやり方をご紹介します。

実験や観測で、直線っぽい、二次曲線っぽい、2変数のガウシアンっぽいデータが得られた時、機械学習なんて使わずに最小二乗法を使ったfittingを行うと大体データを数式で表せます。

pythonでfittingを行うのは、scipyのcurve_fitを使えば、かなり簡単です。

得られた数式の誤差評価(今回紹介するのは、カイ二乗評価と決定係数評価)もpythonなら数行で計算してくれます。

下記では、最小二乗法を使用したfittingの原理や評価方法の概念の紹介を行い、線形fittingを例にpythonでの実装方法を示した後に、多変数の場合も2Dガウシアンを例に説明していきます。同時にfittingの良し悪しを判断できるカイ二乗評価と決定係数評価も実装していきます。

目次

「Fitting」とは

Fittingは、日本語にすると「曲線あてはめ」などとなります。

測定や観測したデータに対して、モデル式のパラメータを変えていき、最もデータに合う(fitする)パラメータ値を決定するのがfittingです。

Fittingの何が嬉しいのかというと、得られたモデル式を使って、実測データがない部分の値を推定(内挿、外挿)できるのです。例えば1日の気温を測定した時に、7時(20℃)と9時(22℃)のデータがあれば、8時の気温が21℃だと推定できます。

fittingの場合、機械学習とは異なり、モデル式は自分で考える必要があります。

モデル式は、理論的に「この式だ!」という式がある場合もあれば、よくわからないけど直線ぽいからとりあえず直線の式で当てはめてみるとか、全くどんな式か見当もつかないからべき乗関数式(テーラー展開)でやってみる、など様々です。

イメージとしては次の図のようになります。

fitting_001

Fitting曲線(モデル式)内のパラメータ\( a, b, c, d \)を変化させて、実測データとモデル式の差が最も少なくなるようなパラメータセットを決定するのが最小二乗法によるfittingです。

ここで「実測データとモデル式の差」とは、次の図中の赤い線の長さの合計(残差平方和)のことです。

この差が少なければ確かに、モデル式と実測データは合っていて(fitしていて)、モデル式を使えば実測データがない部分での\( y \)の予測もできそうなのがわかります。

fitting_002

Fittingの良し悪しの評価方法

Fitting後に、どれだけモデル式がデータに合ったかを定量的に示す方法を2つ紹介します。

カイ二乗

データ値\( y_o \)と、モデル式から計算される予測値\( y_e \)と、を使って次の「カイ二乗」を計算します。

$$ \chi^2 = \sum\frac{(y_o-y_e)^2}{y_e} $$

この数値は、「二乗」と言いつつ負も許される値ですが、これが0に近ければ近いほどに、モデル式がデータ値に合っている(fitしている)ことになります。

色々と難しい数式を追うと、この\( \chi^2 \)の値は、「自由度」の数\( d \)と同程度か小さければ、モデル式の精度がよく、「自由度」の数\( d \)よりも非常に大きければ精度が悪いです。

$$ \chi^2 \leq d \Rightarrow モデル式の精度高い \\ \chi^2 \gg d \Rightarrow モデル式の精度低い $$

「自由度」とは、「データ数 – モデル式のパラメータ数」です。前節の図の例だと、データの数は16で、パラメータの数は\( a, b, c, d \)の4つなので、自由度は12です。

カイ二乗評価の枠組み内には実はもう一つ「p値(p-value)」というものがあり、次のような式で表されます。

$$ p = \frac{2}{2^{d/2}\Gamma\left(\frac{d}{2}\right)}\int_{\chi_0}^{\infty}t^{d-1}e^{-\frac{t^2}{2}}{\rm d}t $$

ここで\( \chi_0 \)は「カイ二乗」の平方根で、\( \Gamma(z) \)はガンマ関数で、\( d \)は自由度です。

このp値は0(=0%)から1(=100%)の間の数値で、モデル式がデータにfitする確率のようなものです。

$$ p\simeq 1 \Rightarrow モデル式の精度高い \\ p\simeq 0 \Rightarrow モデル式の精度低い $$

ここまでガンマ関数が出てきたりと、かなり数学が濃くなってしまいましたが、pythonでこれを計算する時には、何も考えずに次のようにコードを書き込めばい良いだけです。pythonバンザイ。

ここでは自由度を入力していませんが、その場合には「データ数 – 1」が自由度として設定されます。

結果は変数「chi2」に格納されて、chi2[0]がカイ二乗値、chi2[1]がp値となります。

自由度設定などを詳しく知りたい方はここをご参照ください。

決定係数

ここではもう一つの「fittingの良し悪しを評価できる数値」として、決定係数( \( R^2 \) )を紹介します。

決定係数は次のような式で定義されています。

$$ R^2 = 1- \frac{\sum \left( y_o – y_e \right)^2}{\sum \left( y_o – \overline{y_o} \right)^2} $$

ここで\( y_o \) は測定値で、\( y_e \) はモデル式による予測値です。\( \overline{y_o} \) は観測値の平均を表しています。

この数値はイメージがしやすいですね。

$$ R^2 \simeq 1 \Rightarrow モデル式の精度高い \\ R^2 \simeq 0 \Rightarrow モデル式の精度低い $$

実は、この\( R^2 \) も二乗がついているものの、負の値を取り得ます。もちろん1から離れれば離れるほどに、モデル式の精度が悪いことを表しています。

Pythonでこれを計算する時には、単純に次のように計算させます。

Pythonでfittingする方法 <線形関数を例に>

これはかなり簡単で、scipyのcurve_fitというのを使用します。

メインの計算は次の1行でできます。

linear_funcは定義された関数(モデル式)で、x_observedとy_observedが観測データです。

Fittingで得られたパラメータはpoptに格納されますが、その推定誤差の情報はpcovに格納されます。

pcovは実際には共分散行列なので、上記コードのように対角成分の平方根をとることで、推定誤差を算出します。

実際に線形関数をモデル式にした場合の例を示します。

このコードの書き方では、関数をまず定義した後に、if name=="main":の部分から処理が開始されます。

観測データがnumpy配列として「data」に格納されて、fitting_linearに渡されます。

fitting_linear関数内では、モデル式として、linear_funcが定義されていて、curve_fitを使ってfittingが行われています。超簡単!

curve_fitでは、最初の入力変数がモデル式で、2番目がx(独立変数、説明変数)で、3番目がy(従属変数、被説明変数)となっています。

(このxは変数の数としては1つしか入力できないのですが、中身は配列でもよいので、多変数にできます。詳細は後ほど。)

評価値計算とグラフ表示も含めた例示

先述のカイ二乗、p値、決定係数の計算と、測定データおよびfitting結果のグラフ化まで含めたコードを示しておきます。

結果のグラフは次のようになります。

fitting_003

モデル式(linear_func部分)を適当に変えれば、非線形でもなんでもモデル式を作成できます。

多変数のfittingを行う場合

上記の例ではyはxのみの関数でした。つまり独立変数(説明変数)は1つでした。

ここでは独立変数が複数の時にはどのようにfittingをすれば良いかを、2変数を例に説明していきます。

ちなみに、モデル式は2次元のガウシアンです。

$$ z = A\frac{1}{\sqrt{2\pi \sigma_x^2}} e^{-\frac{(x-\mu_x)^2}{2\sigma_x^2}} \frac{1}{\sqrt{2\pi \sigma_y^2}} e^{-\frac{(y-\mu_y)^2}{2\sigma_y^2}} $$

1変数と特に変わったところは、curve_fitへの入力とモデル式を表す関数の入力の仕方です。

curve_fitへの入力は次のようにします。

独立変数(説明変数)は2つをカッコにくくり、1つの変数として入力します。

モデル式の定義では次のようにします。

このようにすることで、入力は1変数Xにして、関数中で2つ(=x, y)に分けてあげます。

多変数の場合のfittingコード例 <評価値計算とグラフ表示を含む>

コードは次のようになります。

Fitting結果のグラフは次のようになります。

fitting_005

これで多変数のfittingもできるようになりました。

最後に

pythonでならfittingは簡単です。数行でかけます。

ただ、今回はカイ二乗値や決定係数といった評価値に関して復習をさせられました。

特にp値を示した式というのは、案外ネット上になく探すのを苦労しました。。。

どうぞ、ご参考にしてください!

関連記事

応援よろしくお願いします☆

この記事を書いた人

天文の博士号をもつ理系パパ。
3歳の娘を子育て中。
最近はダイエットに挑戦中!

コメント

コメントする

目次
閉じる