pythonで多項式近似

曲線フィッティング

なにか観測されたデータに対して、規則性を見つけたり、
予測をしたりするのに有用なのがカーブフィッティングというものです。
そして今回はあてはめる曲線を多項式による曲線でやってみます。

Contents

1. 回帰分析
2. 最小二乗法
3. プログラムの実装・解説
4. まとめ

1. 回帰分析

複数のデータ系列が与えられたとき、それらの対応関係を表すのに尤もらしい関数を求めることが回帰分析です。

例えば、「アイスクリームの売上」と「気温」のデータがあったとき、
なんとなく気温が高いと売上が上がる気がします。
この現象は右肩上がりの1次関数で近似できそうです。

何かの関数に近似ができれば、
今は手元に0℃から30℃まで5℃刻みの気温に対する売上情報しかなくても、
27.5℃など間の数字や38℃などの範囲外の売上が予測できます。

アイスクリームの例については「アイスクリーム屋さんで学ぶ楽しい統計学（
向後 先生）」が元になっています。
直感的にわかりやすい例や易しい解説があり、統計初心者にとっていい入り口になると思います。
ちなみに書籍(amazon)も出ています。

2. 最小二乗法

例えばFig. 1のようなデータが与えられたとします。

観測されるデータは、ある（本来は見えない）数理モデルにもとづいて出力されたものに、
ノイズが加わったものとして現れることが多いです。
このある数理モデルが多項式にもとづいているとして、
近似曲線を求めることを多項式近似といいます。

最小二乗法は近似曲線を求める方法のひとつです。
簡単のために直線（1次多項式）で近似する場合を考えます。

あてはめたい直線と観測されたデータの差が誤差で、
この誤差を二乗して、全部足したものを基準に近似をします(Fig. 2)。

誤差の二乗の和、すなわち観測されたデータ系列と予測している関数の差が
小さくなるように多項式のパラメータを決定します。
誤差は正負とも現れるので、+5000兆と-5000兆という誤差が発生した場合に
そのまま足し合わせると全体の誤差が0になってしまうので、
2乗することで打ち消し合うことを避けています。

こうして得られた二乗誤差の和は、

\[\sum_{i=1}^{n} \{ y_{i} – f(x_{i})\}^2 \]

のように表せます。これを最小になるようなf(x)を見つけます。

近似する直線の式は、

\[y=ax+b\]

の形で表せます。したがって調整するパラメータはaとbです。
調整すべきパラメータは次数によって決まります。
（n次多項式ならn+1のパラメータ）

最小二乗法はシンプルでメジャーな方法なので、
深く知りたい人は調べるとたくさん資料が出てきます。
知らなくてもなんとかなってしまいますが…

3. プログラム

現象の元となる関数の生成

a = 1.0
b = 2.0
c = 3.0

target_min = 0.0
target_max = 4.0

f = np.poly1d([a, b, c], True)
x = np.linspace(target_min, target_max, int(target_max-target_min)*20)

print(f)
plt.plot(x,f(x))

numpyのpoly1dで多項式を作れます。
Trueなしでやると、

\[a \times x^2 + b \times x + c \]

が得られ、上に書いてるようにTrueありだと、

\[(x-a)(x-b)(x-c)\]

が得られます。後者ならy=0と交わる点をもとに簡単に決められます。

観測される（であろう）データの生成

上で得られた関数をもとにノイズを加えます。

ave = 0.0
sd = 1.0
data_N = 7

y=[]
data_x = np.linspace(target_min, target_max, data_N)
for i in range(len(f(data_x))):
    y.append(f(data_x)[i] + random.normal(loc=0.0, scale=1.0))

plt.xlabel("x")
plt.ylabel("y")
plt.scatter(data_x, y);

numpy.normalは正規分布にもとづく乱数を生成しています。
観測できるデータは有限なのでdata_Nとして設定しました。

多項式近似

order = 1

estimated_parameters = np.polyfit(data_x, y, order)
estimated_curve = np.poly1d(estimated_parameters)

plt.xlabel("x")
plt.ylabel("y")
plt.plot(x, estimated_curve(x))
plt.scatter(data_x, y);

numpyのpolyfitで簡単に多項式の近似曲線を求めることができます。
（中で動いているのは最小二乗法です。）
orderという変数によってあてはめる多項式の次数を決定しています。

本来は3次曲線にのっとっているのに、無理に複雑なモデル
（今回でいうと5次）で近似をすると間違った解釈をしてしまいます。
これは観測されたデータに対してのみ過剰に適合している、と言えます。

誤差を見る（評価）

過剰に適合している場合があるので、必ずしもモデルの評価になるとも
限りませんが平均の二乗誤差を求めてみるコードを追加しました。

rate = int(data_N / (target_max-target_min)*20)

square_error = 0
for i in range(data_N):
    square_error += (estimated_curve[int(rate*i)] - y[i]) ** 2
mse = square_error / data_N
print(mse)

誤差をとって、2乗して、それらを繰り返し足して、平均しただけです。

4. まとめ

polyfitを使って多項式近似をしました。
データ点数（data_N）や多項式次数（order）を変えて、
グラフや二乗誤差がどう変わるかなどみると面白いかも、と思いました。

今回のプログラムもGithubにアップロードされていますので、
ぜひこちらもご利用ください。

---