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

<a href="https://qiita.com/maskot1977/items/2a55c087974d0ac5e843">様々な補間法と最小2乗法をPythonで理解する</a> のうち、「Numpy.polyfit を使ったカーブフィッティング」を、実データっぽい模擬データを解析するように書き直したサンプルプログラムです。

# 実データっぽい模擬データ

こういうデータを多項式近似したいとしましょう。

```python
x_observed = [9, 28, 38, 58, 88, 98, 108, 118, 128, 138, 148, 158, 168, 178, 188, 198, 208, 218, 228, 238, 278, 288, 298]
y_observed = [51, 80, 112, 294, 286, 110, 59, 70, 56, 70, 104, 59, 59, 72, 87, 99, 64, 60, 74, 151, 157, 57, 83]
```

# まずはデータの図示

こういうデータですね。

```python
%matplotlib inline
import matplotlib.pyplot as plt

plt.scatter(x_observed, y_observed)
plt.grid()
```

![output_1_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/823fefba-ea01-ae52-3ef0-ea98d23c9535.png)


# フィッティンブカーブを描く x 座標を作る

実データの最小値から最大値までの間に等間隔に１００ポイント作りました。


```python
import numpy as np
x_latent = np.linspace(min(x_observed), max(x_observed), 100)
```

# Numpy.polyfit を用いた最小二乗法

たとえば１次式から９次式までフィッティングしてみましょう。


```python
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)]
```

# Sympy を用いた数式の表示の準備

Sympy の使い方は <a href="https://qiita.com/maskot1977/items/b4395da5f33f70cd4a09">微分や微分方程式をPythonで理解する</a> などを参考にしていただければと。

```python

import sympy as sym
from sympy.plotting import plot
sym.init_printing(use_unicode=True)

x, y = sym.symbols("x y")
```

# 最小二乗法によるカーブフィッティング（１次式から９次式まで）

数式の表示とグラフの表示を一気に行います。

```python
for method_name, method in [cf1, cf2, cf3, cf4, cf5, cf6, cf7, cf8, cf9]:
    print(method_name)
    # 係数の計算
    coefficients = method(x_observed, y_observed)
    
    # Sympy を用いた数式の表示
    expr = 0
    for index, coefficient in enumerate(coefficients):
        expr += coefficient * x ** (len(coefficients) - index - 1)
    display(sym.Eq(y, expr))
    
    # プロットと曲線の表示
    fitted_curve = np.poly1d(method(x_observed, y_observed))(x_latent)
    plt.scatter(x_observed, y_observed, label="observed")
    plt.plot(x_latent, fitted_curve, c="red", label="fitted")
    plt.grid()
    plt.legend()
    plt.show()
```

    最小2乗法（1次式）


![output_5_2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/724fc183-16fb-2416-05a8-375b3724b91a.png)


![output_5_3.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/c6edcd98-8326-6839-1bda-6a4103d52de2.png)



    最小2乗法（2次式）


![output_5_5.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/de9d20da-58b6-9b6f-f800-a5b9e74fa9b7.png)


![output_5_6.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/b9b86edd-8b68-925d-89ed-06260632c86d.png)




    最小2乗法（3次式）



![output_5_8.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/231964f1-5480-268c-5799-fe343b4b99ba.png)


![output_5_9.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/6f55f222-5f49-7e47-ee9f-da230acfa624.png)


    最小2乗法（4次式）



![output_5_11.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/adf50343-ebf4-7537-0235-9cf931c92706.png)


![output_5_12.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/c313c3b9-b7b6-0e55-5c4a-fe6e5fad662d.png)



    最小2乗法（5次式）


![output_5_14.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/1e2677a8-0d7d-5006-bfc8-071ce70a943c.png)


![output_5_15.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/313bad92-8e03-56cb-927a-fce395a64502.png)



    最小2乗法（6次式）



![output_5_17.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/46f590f2-f15b-1cf7-491f-1c7385cec2dd.png)


![output_5_18.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/5bbaad3d-a9fe-d937-6d5d-440ed2afa70b.png)



    最小2乗法（7次式）


![output_5_20.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/ca56bf73-de09-ca5d-0ff1-7a312dba90d0.png)

![output_5_21.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/7d441a45-90fa-e153-51e3-af78b4cc7a33.png)



    最小2乗法（8次式）


![output_5_23.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/5a0d2608-622d-0d55-2cf3-40be048e774a.png)


![output_5_24.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/86d2ac3a-1e93-b774-e0c2-5bf2829ab555.png)


    最小2乗法（9次式）



![output_5_26.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/39014058-7409-63f0-c1df-8b859a2f38b1.png)


![output_5_27.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/af8f0137-cb39-099c-1a98-c2afd170ea30.png)

# 参考

<a href="https://qiita.com/maskot1977/items/913ef108ff1e2ba5b63f">Scipy.interpolate を使った様々な補間法</a>もご参考にどうぞ。