【データ分析】過渡現象のグラフをフィッティングしてみた！

#概要
学生実験などでレポートを書く時には、手軽に使えるExcelを使う機会も多いと思います。
Excelを使えば、グラフも簡単にかけ、簡単なグラフであれば、近視曲線まですぐに計算してくれて便利ですよね。
多項式で近似できる曲線であればよいのですが、指数関数とか対数関数とかだと
Excelで近似曲線を出そうと思ったら少しコツが必要だったり、Excelでは計算できないなんてこともありますよね。

今回は、電気回路などでよく出てくる過渡現象の近似曲線をPythonを使って計算し、近似曲線の式まで
計算してみようと思います。これでexpを含む式もすぐに近似できるようになりますね。

#動作環境
- Windows10(64bit)
- Python 3.7.2

#近似曲線を計算する関数
Pythonには、近似曲線を計算できる関数がいくつかあるようです。

結論から言うと、scipy.optimizeの**curve_fit**を使うといいと思います。

Pythonで近似曲線を計算してくれる関数は主に以下の3つがあるようです。(そのほかにもいくつかありそうです)
方法は全て最小二乗法っぽいです．

- numpy.polyfit()
- scipy.optimize.leastsq()
- scipy.optimize.curve_fit()

一つずつ順番に説明していきます。

- numpy.polyfit()

これは名前からもわかるように多項式(polynomial)で近似する関数です。なので、非常に簡単に使える反面、使える場面が限られてくると思います。

使い方は簡単で、第一引数に独立変数xの配列、第二引数に従属変数yの配列、第三引数に近似式の次数(degree)を取ります。
返り値は、近似多項式の係数(coefficient)で、高次のものから低次のものの順になります．

たとえば，以下のように使います．

```py
x = np.array([1,2,3])
y = np.array([2,7,16]) # y=2x^2-x+1

coeff = np.polyfit(x, y, 2) # coeff = [ 2. -1.  1.]
```

- scipy.optimize.leastsq()

これは名前からまさに最小二乗法っぽいですね．ただ，この関数使い方が少し独特だと私は感じました．

引数は2つ．第1引数は関数，第2引数は推定における初期値です．返り値は，係数の配列と次数のようです．

たとえば，以下のように使います．

```py
x = np.array([1,2,3])
y = np.array([2,7,16]) # y=2x^2-x+1

def func(param, x, y):
    residual = y -(param[0]*x**2 + param[1]*x + param[2])
    return residual

param[0,0,0]
coeff = optimize.leastsq(func, param, args = (x, y)) # coeff = (array([ 2., -1.,  1.]), 2)
```

- scipy.optimize.curve_fit()

最後に本命のcurve_fit()です．

使い方は，第1引数にパラメータをもつ関数，第2引数に独立変数x，第3引数に従属変数yを取ります．
返り値は，近似式のパラメータの配列と共分散です．

具体的には以下のように使います．

```py
x = np.array([1,2,3])
y = np.array([2,7,16]) # y=2x^2-x+1

def func(x, a, b, c):
    return a*x**2 + b*x + c

param, cov = curve_fit(func, x, y) 
# param = [ 2. -1.  1.]
# cov = 
[[inf inf inf]
 [inf inf inf]
 [inf inf inf]]
```

#コード
では，早速実装していきましょう．

```py:transient_phenomena.py
from scipy.optimize import curve_fit
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt

def fit(x,K,T):
    return K*(1-np.exp(-x/T)) # 近似したい関数の式

K = 10 # 比例定数
T = 1 # 時定数

list_x = range(11)
array_x = np.array(list_x)

list_y_data = []
for num in array_x:
    list_y_data.append(K*(1-np.exp(-num/T)) + np.random.rand())
array_y_data = np.array(list_y_data)

param, cov = curve_fit(fit, array_x, array_y_data)
print(param)
array_y_fit = param[0] *array_x + param[1]

list_y_fit = []
for num in array_x:
    list_y_fit.append(param[0]*(1-np.exp(-num/param[1])))
array_y_fit = np.array(list_y_fit)

sns.pointplot(x=array_x, y=array_y_data, join=False)
sns.pointplot(x=array_x, y=array_y_fit, markers="")

plt.show()
```

#解説

- フィッティングしたいデータの作成

```py
list_x = range(11)
array_x = np.array(list_x)

list_y_data = []
for num in array_x:
    list_y_data.append(K*(1-np.exp(-num/T)) + np.random.rand())
array_y_data = np.array(list_y_data)
```

まずは，フィッティングしたデータを用意します．実験を行ったのであれば，そのデータになりますが，
今回は，フィッティングの精度も確かめたいため乱数を使ってデータを作成します．

np.random.rand()で0.0以上1.0未満の一様分布の乱数を発生できます．

- フィッティングを行う

```py
param, cov = curve_fit(fit, array_x, array_y_data)
print(param)
array_y_fit = param[0] *array_x + param[1]

list_y_fit = []
for num in array_x:
    list_y_fit.append(param[0]*(1-np.exp(-num/param[1])))
array_y_fit = np.array(list_y_fit)
```
先に述べたcurve_fit()を用いてフィッティングを行います．

# 実行結果
実行結果は以下のようになりました．乱数を用いているので，実行時によって結果は多少変わると思います．

![transient_phenomena.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/433790/2d57540a-2347-0a18-d19d-13cf1f2bf6d7.png)

今の場合，
coeff = [10.5086373   1.00943089]
となっています．うまくいっていますね．

#まとめ
いかがだったでしょうか？今回は過渡現象のフィッティングを行ってみました．
次はもっと複雑な関数の近似を行ったり，物理現象に限らず回帰分析分野にも
挑戦していきたいと思います．

今回のコードのなどについて改善点や建設的なフィードバックなどありましたらコメントに書いていただけると
助かります！