SciPyのcurve_fitでデータを任意の関数に近似する方法

  • このエントリーをはてなブックマークに追加

技術計算の分野では、測定されたデータを任意の関数にカーブフィットする需要が頻繁にあります。Pythonのscipy.optimize.curve_fitを使えば点列データを1次元や2次元の関数で簡単にフィッティングできます。ここでは様々な関数を例にcurve_fitを使ってみた内容を紹介します。

こんにちは。wat(@watlablog)です。本日は超便利な関数近似メソッドであるscipy.optimize.curve_fitを使ってみた結果を紹介します

本記事の目標

1Dと2Dデータに対するcurve_fitの使い方を理解する

線形近似から学ぶ

関数近似入門としての代表例は最小二乗法による線形近似です。線形近似については以下の記事で数式展開もしながら解説していますが、直線の方程式でデータを近似するだけでも多くの分析が可能なデータ分析の王道手法です。

Advertisements

まずはこの内容をscipy.optimize.curve_fitでできるようにして、基礎を理解しましょう。
Pythonでカーブフィット!最小二乗法で直線近似する方法

様々な関数で応用力を付ける

線形近似の次は多項式近似に進み、その後は指数関数対数関数による近似をしてみます。この辺りでExcelのグラフ上で行う関数近似は網羅できるのですが、この記事ではさらに三角関数確率密度分布がそれぞれ近似できるような内容も載せました。

scipy.optimize.curve_fitは任意関数として設定する方程式を変えるだけでこれらの近似ができるので、ここまでくれば1Dのどんなデータが来ても対応できるでしょう。

曲面フィッティングもやってみる

説明変数が2つの重回帰分析をはじめとする曲面で近似したい場面も多くあります。ここでは2Dのデータに対して自分で設定した曲面を表現する関数にどうフィッティングさせるか、その結果の可視化方法も併せて紹介します。

直接的な可視化限界である2変数問題が扱えるようになれば、あとは超平面へのカーブフィットもコーディングできるようになるはずです(ここでは扱いませんが、可視化は結局写像して3Dプロットとかかな?)。

色々やりますが、結局は公式で紹介されていることをベースにしているだけです。公式ドキュメントはこちら↓。いつでも1次情報を確認するようにしましょう。
https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html

動作環境

本記事のコードは以下の環境で動作を確認しました。

Python Python 3.9.6
PyCharm (IDE) PyCharm CE 2020.1
Numpy 1.23.0
Scipy 1.4.1
matplotlib 3.7.1

1Dデータの関数近似例

多項式近似

線形近似(初回詳細解説):\(ax+b\)

直線の方程式(1)を使ってコードの基礎を確認します。

\[ y=ax+b \tag{1} \]

線形近似(線形回帰)を行うコードは段階を追ってコードを見ていきましょう。
まずはimport文を書きます。SciPyは膨大なパッケージですので、from scipy.optimize import curve_fitと必要なモジュールだけ取り込みます。

任意関数として直線の方程式を関数の形で記載します。第1引数xが変数で、その後に式に使う係数を続けます。return文に計算後の結果を渡すという至極一般的な関数の書き方です。

テストするためのデータセット(点列データ)を準備しましょう。以下のコードは先ほど作った関数funcに引数を渡してデータを作っているものです。ただ、それだけだと方程式通りの綺麗なデータになってしまい面白くないので、np.random.normalでガウス分布を持つノイズを混入させています。

ランダム成分があるので毎回完全一致はしませんが、下図の結果が得られるはずです(可視化の方法は後述します)。

y=ax+bのデータセット

そしていよいよカーブフィットを行います。カーブフィットはcurve_fit()に近似に使う関数funcそのものとデータセットx, yを与えるだけです。

結果がpoptpcovの2つあります。poptは最適化されたパラメータ(\(p_{opt}\))、pcovはパラメータの分散共分散行列(共分散=covariance)です。

pcovの対角成分を使えば以下の文でパラメータ毎に誤差を調べることができます。

この誤差errorの値が大きいほど、データによくフィットできていないということがわかります。

最後に最適化されたパラメータを使って近似された関数の形を確認します。poptは任意に設定した関数の引数の数に応じてサイズが変わりますが、*poptとアスタリスクを付けて書くことで一度に複数の引数を設定できます。つまり、関数によってこの部分の書き方を変えなくても良いというメリットがあります。

こちらがデータと近似結果を比較したプロットです。イメージ通りですね。

y=ax+bのカーブフィット結果

matplotlibによる可視化方法も含めた全コードを以下に示します。コピペで動作するはずなので是非ご自分の環境で動かしてみてください。

\(ax^2+bx+c\)

基本的な使い方がわかったので、式を変更してどんどんカーブフィットを試してみましょう。次は2次方程式(2)を近似します。

\[ y = ax^2+bx+c \tag{2} \]

コードは式とデータセット生成部分を変更するだけですが、どこから読みはじめてもコピペで動作できるように毎回全コードを載せます。

良い近似結果です!

2次方程式のカーブフィット結果

\(ax^3+bx^2+cx+d\)

続いて3次方程式(3)です。

\[ y = ax^3+bx^2+cx+d \tag{3} \]

良い近似…(以下略)

3次方程式のカーブフィット結果

指数関数近似/対数関数近似

\(e^{x}\)

続いて指数関数(4)の近似をしてみます。式全体に係数をつけていますが、xにも何かかけても良いかも知れません。

\[ y = a e^{x} \tag{4} \]
指数関数のカーブフィット結果

\(\log{x}\)

続いて対数関数(5)を近似します。

\[ y = a \log{x} \tag{5} \]

対数関数はxの範囲に0をとってしまうとエラーになるので、その部分もちょっと変更しています。

対数関数のカーブフィット結果

三角関数近似

\(\sin x\)

簡単過ぎてもはや単なる作業になってきましたが、近似したい関数の中でも人気(?)なサインカーブ(6)も近似してみましょう。波の特性がわかりやすくなるように、xを時間と仮定して周波数fの概念を入れてみました。

\[ y = a \sin(2 \pi f x) \tag{6} \]
サインカーブのカーブフィット結果

\(\cos x\)

サインをやっておいてコサインが無いのも気持ちが悪いので、コサインカーブ(7)も載せておきます。ただこちらは上記func()関数内sin部分をcosにするだけなのでコードは割愛し、結果だけお見せします。

\[ y = a \cos(2 \pi f x) \tag{7} \]
コサインカーブのカーブフィット結果

\(\tan x\)

三角関数シリーズでタンジェントだけ仲間はずれにするのは可哀想ということで、タンジェントカーブ(8)も近似しましょう。

\[ y = a \tan(2 \pi f x) \tag{8} \]

\(\tan\)はy軸が大きくなりすぎるので、matplotlibにスケール設定を追加してみました。

タンジェントカーブのカーブフィット結果

他にも三角関数や指数関数の親戚である双曲線関数(\(\sinh\)とか)も関数func部分を変更するだけで試すことができます。ここでは割愛しますが、是非お好みの関数で近似してみてください。

1次元確率密度関数近似

実験データを多数取得した時、データを任意の値域(ビン)でカウントしたヒストグラムを描くことが多いです。ヒストグラムに対して平均や標準偏差を調べることは統計的に実験データを解釈するのに必要なので、確率密度関数でデータを近似するのは非常に重要なことです。

式(9)はガウス分布を持つ確率密度関数です。代表的な例として、この式で近似するコードを紹介します。

\[ y = a e^{(-\frac{(x - \mu)^{2}}{(2 \sigma^{2})})} \tag{9} \]

このコードは先ほどまでと異なり、ヒストグラムを作成する関数から書いています。その後、その作成されたデータを用いてカーブフィットをするという流れです。

こちらが近似結果です。綺麗に近似できていますね!poptを見れば平均や標準偏差の値を調べることができます(パラメータとして設定してあるので)。

確率密度関数のカーブフィット結果

2Dデータの関数近似例

次は1次元ではなく、2次元のデータに対して近似を行う方法を試してみましょう。

2次元確率密度関数近似

先ほどの続きで確率密度関数の2次元版をやってみます。近似する式は式(10)です。

\[ z = a e^{(-\frac{(x - \mu_{x})^{2}}{(2 \sigma^{2}_{x})})} e^{(-\frac{(y - \mu_{y})^{2}}{(2 \sigma^{2}_{y})})} \tag{10} \]

こちらがコードです。2次元の場合は可視化するためにnp.meshgridを使っていたり、ややデータ整形にややこしいことがあります。地味にnp.histogram2d等も使っていますが、これは理想のデータを生成するためのコードなので、実際は実験データ等をそのまま使うことが可能です。

こちらが結果です。うーん、いいですね!

2D確率密度関数のカーブフィット結果

2D平面近似:\(z = w_{0} + w_{1} x + w_{2}y\)

次は重回帰分析でよくお世話になる2次元平面の式(11)です。

\[ z = w_{0} + w_{1} x + w_{2}y \tag{11} \]

こちらも綺麗にフィッティングできています!

2D平面のカーブフィット結果

\(z = \sin (x) \cos(y)\)

2Dデータに対するフィッティングがわかってきたので、最後はお気に入りの式(12)を近似してみます。

\[ z = a \sin(x) + b \cos(y) \tag{12} \]

コードはこちら。一度できたら小変更で済むのがプログラミングのすごいところですね。

2D波のカーブフィット結果

まとめ

本日はここまで!
この記事ではscipy.optimize.curve_fitによる関数近似の事例をまとめました。1次元だけではなく、2次元の点列データに対するカーブフィッティングも簡単に行うことができました。
これは久しぶりに汎用性の高いメソッドを見つけたと思うので、実際のデータ分析業務で活用ができそうです。

任意関数のカーブフィットはもはやこれで十分なのでは??
X(旧Twitter)でも関連情報をつぶやいているので、wat(@watlablog)のフォローお待ちしています!

  • このエントリーをはてなブックマークに追加

SNSでもご購読できます。

コメントを残す

*