Python
scipy

様々な補間法と最小2乗法をPythonで理解する

実験などにより得られた観測値は、普通は飛び飛びの値になりますが、その間の値を求めたい時があります。その時に用いるのが、種々の補間法(補完ではありません)と、その他のカーブフィッティング(曲線近似)です。これらの解法も、プログラミング演習の題材としてよく使われますが、実用上は強力なライブラリが多数存在しますので、自作するよりも既存のライブラリを使ったほうが便利です。

ただし、それらの手法の中身を理解しておくこともとても大事ですし、その実装はプログラミングの勉強にとても良いです。

まずは、強力なライブラリ群を紹介してその使い方を習得しましょう。


計算機実験

ここでは、計算機実験として、以下の3つの数学関数 f(x), g(x), h(x) を用意します。

import numpy as np

f = lambda x: 1/(1 + np.exp(-x))
g = lambda x: 1.0/(1.0+x**2)
h = lambda x: np.sin(x)

上記3つの数学関数が実際に起きている現象であるものとします。


実験による観測値

ですが実際には、実験で得られる観測点は、次のように飛び飛びの値です。

x_observed = np.linspace(-10, 10, 11)

x_observed

array([-10.,  -8.,  -6.,  -4.,  -2.,   0.,   2.,   4.,   6.,   8.,  10.])

このとき、実験によってどのような値が観測されるか図示してみましょう。

%matplotlib inline

import matplotlib.pyplot as plt

for func in [f, g, h]:
y_observed = func(x_observed)
plt.scatter(x_observed, y_observed)
plt.grid()
plt.show()

output_5_0.png

output_5_1.png

output_5_2.png

以上のような観測値から、どのようなカーブ(曲線)が想像できるでしょうか。


観測されない真実の姿

実験では観測されない真実の姿を図示してみましょう。

x_latent = np.linspace(-10, 10, 101)

x_latent

array([-10. ,  -9.8,  -9.6,  -9.4,  -9.2,  -9. ,  -8.8,  -8.6,  -8.4,

-8.2, -8. , -7.8, -7.6, -7.4, -7.2, -7. , -6.8, -6.6,
-6.4, -6.2, -6. , -5.8, -5.6, -5.4, -5.2, -5. , -4.8,
-4.6, -4.4, -4.2, -4. , -3.8, -3.6, -3.4, -3.2, -3. ,
-2.8, -2.6, -2.4, -2.2, -2. , -1.8, -1.6, -1.4, -1.2,
-1. , -0.8, -0.6, -0.4, -0.2, 0. , 0.2, 0.4, 0.6,
0.8, 1. , 1.2, 1.4, 1.6, 1.8, 2. , 2.2, 2.4,
2.6, 2.8, 3. , 3.2, 3.4, 3.6, 3.8, 4. , 4.2,
4.4, 4.6, 4.8, 5. , 5.2, 5.4, 5.6, 5.8, 6. ,
6.2, 6.4, 6.6, 6.8, 7. , 7.2, 7.4, 7.6, 7.8,
8. , 8.2, 8.4, 8.6, 8.8, 9. , 9.2, 9.4, 9.6,
9.8, 10. ])

fig_idx = 0

plt.figure(figsize=(12,4))
for func in [f, g, h]:
fig_idx += 1
y_observed = func(x_observed)
y_latent = func(x_latent)
plt.subplot(1, 3, fig_idx)
plt.scatter(x_observed, y_observed)
plt.plot(x_latent, y_latent)
plt.grid()
plt.show()

output_8_0.png

さて、ここで求めたいことは、限られた観測値から、この「真実の姿」にどれだけ近いものを導き出せるかということです。


Scipy.interpolate を使った様々な補間法

scipy には interpolate という強力なライブラリが用意されていて、様々な補間法が使えます。

from scipy import interpolate

ip1 = ["最近傍点補間", lambda x, y: interpolate.interp1d(x, y, kind="nearest")]

ip2 = ["線形補間", interpolate.interp1d]
ip3 = ["ラグランジュ補間", interpolate.lagrange]
ip4 = ["重心補間", interpolate.BarycentricInterpolator]
ip5 = ["Krogh補間", interpolate.KroghInterpolator]
ip6 = ["2次スプライン補間", lambda x, y: interpolate.interp1d(x, y, kind="quadratic")]
ip7 = ["3次スプライン補間", lambda x, y: interpolate.interp1d(x, y, kind="cubic")]
ip8 = ["秋間補間", interpolate.Akima1DInterpolator]
ip9 = ["区分的 3 次エルミート補間", interpolate.PchipInterpolator]

for method_name, method in [ip1, ip2, ip3, ip4, ip5, ip6, ip7, ip8, ip9]:

print(method_name)
fig_idx = 0
plt.figure(figsize=(12,4))
for func in [f, g, h]:
fig_idx += 1
y_observed = func(x_observed)
y_latent = func(x_latent)
fitted_curve = method(x_observed, y_observed)
plt.subplot(1, 3, fig_idx)
plt.scatter(x_observed, y_observed, label="observed")
plt.plot(x_latent, fitted_curve(x_latent), c="red", label="fitted")
plt.plot(x_latent, y_latent, label="latent")
plt.grid()
plt.legend()
plt.show()

最近傍点補間

output_13_1.png

線形補間

output_13_3.png

ラグランジュ補間

output_13_5.png

重心補間

output_13_7.png

Krogh補間

output_13_9.png

2次スプライン補間

output_13_11.png

3次スプライン補間

output_13_13.png

秋間補間

output_13_15.png

区分的 3 次エルミート補間

output_13_17.png


Numpy.polyfit を使ったカーブフィッティング

Numpy の .polfyfit を用いると、最小2乗法によるカーブフィッティングが簡単にできます。改めて、実験により観測された飛び飛びの観測値を計算し直しておくと

x_observed = np.linspace(-10, 10, 11)

x_observed

array([-10.,  -8.,  -6.,  -4.,  -2.,   0.,   2.,   4.,   6.,   8.,  10.])

y_observed = f(x_observed)

y_observed

array([4.53978687e-05, 3.35350130e-04, 2.47262316e-03, 1.79862100e-02,

1.19202922e-01, 5.00000000e-01, 8.80797078e-01, 9.82013790e-01,
9.97527377e-01, 9.99664650e-01, 9.99954602e-01])

上記の x, y に対して、最小2乗法による線形回帰を行うと

np.polyfit(x_observed, y_observed, 1)

array([0.06668944, 0.5       ])

上記の戻り値がそのまま y = ax + b の係数 a, b に相当します。つまり

y = 0.06668944 x + 0.5

というのが、得られた回帰式です。

最後の引数で何次式に回帰するかが指定できます。たとえば3次式にフィッティングしたければ

np.polyfit(x_observed, y_observed, 3)

array([-0.00074604,  0.        ,  0.11980739,  0.5       ])

y = -0.00074604 x3 + 0.11980739 x + 0.5

に回帰できたことを示しています。

得られたカーブを描くためには、複数のxとそれに対応するyが必要ですが、それは以下のように得られます。

np.poly1d(np.polyfit(x_observed, y_observed, 3))(x_observed)

array([ 0.04796474, -0.07648735, -0.05770001,  0.0685169 ,  0.26635352,

0.5 , 0.73364648, 0.9314831 , 1.05770001, 1.07648735,
0.95203526])

では、先ほどの3つの数学関数 f(x), g(x), h(x) を1〜9次式でフィッティングして図示してみましょう。

cf1 = ["最小2乗法(1次式)", lambda x, y: np.polyfit(x, y, 1)]

cf2 = ["最小2乗法(2次式)", lambda x, y: np.polyfit(x, y, 2)]
cf3 = ["最小2乗法(3次式)", lambda x, y: np.polyfit(x, y, 3)]
cf4 = ["最小2乗法(4次式)", lambda x, y: np.polyfit(x, y, 4)]
cf5 = ["最小2乗法(5次式)", lambda x, y: np.polyfit(x, y, 5)]
cf6 = ["最小2乗法(6次式)", lambda x, y: np.polyfit(x, y, 6)]
cf7 = ["最小2乗法(7次式)", lambda x, y: np.polyfit(x, y, 7)]
cf8 = ["最小2乗法(8次式)", lambda x, y: np.polyfit(x, y, 8)]
cf9 = ["最小2乗法(9次式)", lambda x, y: np.polyfit(x, y, 9)]

for method_name, method in [cf1, cf2, cf3, cf4, cf5, cf6, cf7, cf8, cf9]:

print(method_name)
fig_idx = 0
plt.figure(figsize=(12,4))
for func in [f, g, h]:
fig_idx += 1
y_observed = func(x_observed)
y_latent = func(x_latent)
fitted_curve = np.poly1d(method(x_observed, y_observed))(x_latent)
plt.subplot(1, 3, fig_idx)
plt.scatter(x_observed, y_observed, label="observed")
plt.plot(x_latent, fitted_curve, c="red", label="fitted")
plt.plot(x_latent, y_latent, label="latent")
plt.grid()
plt.legend()
plt.show()

最小2乗法(1次式)

output_25_1.png

最小2乗法(2次式)

output_25_3.png

最小2乗法(3次式)

output_25_5.png

最小2乗法(4次式)

output_25_7.png

最小2乗法(5次式)

output_25_9.png

最小2乗法(6次式)

output_25_11.png

最小2乗法(7次式)

output_25_13.png

最小2乗法(8次式)

output_25_15.png

最小2乗法(9次式)

output_25_17.png


Scipy.optimize を使う方法

カーブフィッティングを行う方法は他にも、Scipy.optimize を使う方法があります。


Scipy.optimize.curve_fit

Scipy.optimize.curve_fit がカーブフィッティング用の関数ですが、Numpy.polyfit のほうが正直使いやすいです。図示は省略します。

from scipy import optimize

eq1 = lambda x, a0, a1: a0 + a1 * x
optimize.curve_fit(eq1, x_observed, y_observed)

(array([-2.18114415e-12, -1.86569661e-03]), array([[ 0.61148173, -0.        ],

[ 0. , 0.00138973]]))


Scipy.optimize.leastsq

Scipy.optimize.leastsq を使ってもカーブフィッティングできますが、Numpy.polyfit のほうが正直使いやすいです。図示は省略します。

from scipy import optimize

eq1 = lambda a, x, y: y - a[0] - a[1] * x
param = [0, 0]
optimize.leastsq(eq1, param, args=(x_observed, y_observed))

(array([-1.33897935e-16, -1.86569660e-03]), 1)


課題


課題1

インターネットで検索して、以下の問いに答えてください(昨年の資料を参考にしてもらっても結構です)。


  • 補間法と、最小2乗法によるカーブフィッティングの違いは何か、説明してください。

  • ラグランジュ補間とスプライン補間の違いは何か、説明してください。


課題2

(Scipy.interpolateもNumpy.polyfitもScipy.optimizeも使わない問題)


  • ラグランジュ補間を行う Python プログラムをインターネットで検索し、それを使って上記の f(x), g(x), h(x) を補間してください(昨年のfortranによるプログラムを参考に自作してもらっても結構です)。そのプログラムが記載してあったURLを明記してください。ただし、Scipy.interpolateもNumpy.polyfitもScipy.optimizeも使わないこと。


  • スプライン補間を行う Python プログラムをインターネットで検索し、それを使って上記の f(x), g(x), h(x) を補間してください(昨年のfortranによるプログラムを参考に自作してもらっても結構です)。そのプログラムが記載してあったURLを明記してください。ただし、Scipy.interpolateもNumpy.polyfitもScipy.optimizeも使わないこと。


  • 最小2乗法を行う Python プログラムをインターネットで検索し、それを使って上記の f(x), g(x), h(x) を1〜9次式で近似してください。そのプログラムが記載してあったURLを明記してください。ただし、Scipy.interpolateもNumpy.polyfitもScipy.optimizeも使わないこと。