フィット結果のグラフをエラー付きで表示する（Sigmoid編）

# はじめに

データ点をSigmoid関数でフィットした結果をエラー付きで表示したいと思い、やってみました。
なお、エラーの計算は地道に手で微分して実装しています。微分や誤差の計算ツールを紹介するような記事ではないので、ご了承ください。

Sigmoid関数は、[Wikipedia](https://ja.wikipedia.org/wiki/%E3%82%B7%E3%82%B0%E3%83%A2%E3%82%A4%E3%83%89%E9%96%A2%E6%95%B0)などに記載されているように

```math

f(x) = \frac{1}{1+e^{-ax}}

```

という関数です。このままだと味気ないので、下のように

```math

f(x) = \frac{p_{0}}{1+e^{-p_{2
}(x - p_{1})}}

```

3つのパラメータ（$ p_{0}, p_{1}, p_{2}$）をもつ関数をフィットに使いたいと思います。

ここでやりたいことは、以下の3つになります。

- データ点に対して関数をフィットし、パラメータを決める
- フィット時のパラメータの共分散を誤差伝搬して、関数のエラーを求める
- グラフとして表示する



# 環境

実行した環境は以下の通り。

```
$sw_vers
ProductName:    Mac OS X
ProductVersion: 10.13.6
BuildVersion:   17G14042
```

Jupyter Notebookを使用しました。

```
The version of the notebook server is: 5.7.8
The server is running on this version of Python:
Python 3.7.3 (default, Mar 27 2019, 16:54:48) 
[Clang 4.0.1 (tags/RELEASE_401/final)]
```

# やったこと

## 準備

フィットには```scipy.optimize.curve_fit```を使用します。
他にも必要なライブラリをインポートしておきます。

```python
from scipy.optimize import curve_fit
import numpy as np
import matplotlib.pyplot as plt
```

フィット用のクラスを用意しました。
```fit```でSigmoid関数でデータ点をフィットしています。```func_sigmoid_error```でSigmoid関数のエラーを計算していますが、詳しいことは後述します。

```python
class Fitter:
    nparams = 3
    isFitted = False
    
    par = {}
    cov = np.zeros((nparams,nparams))
    err = np.zeros(nparams)
    
    def __init__(self):
        self.isFitted = False
        for i in range(self.nparams):
            self.par['p'+str(i)] = 0
    
    # フィットに使用する関数
    def sigmoid(self, x, a):
        return 1.0 / (1.0 + np.exp(-a * x))

    def func_sigmoid(self, x, a, b, c):
        return a * self.sigmoid(x-b, c)
    
    # フィット実行
    def fit(self, xdata, ydata, p0=None, sigma=None, bounds=None):
        par, self.cov = curve_fit(f=self.func_sigmoid, xdata=xdata, ydata=ydata, p0=p0, sigma=sigma, bounds=bounds)
        for i in range(self.nparams):
            self.par['p'+str(i)] = par[i]
        self.err = np.sqrt(np.diag(self.cov))
        self.isFitted = True
        
        print("params: ", self.par)
        print("cov. m: ", self.cov)
        print("error : ", self.err)
        
    
    # Sigmoidの微分
    def dfdp0(self, x):
        return (self.func_sigmoid_fitted(x) / self.par['p0'])
    
    def dfdp1(self, x):
        return (self.func_sigmoid_fitted(x) * (self.par['p0'] - self.func_sigmoid_fitted(x)) / 
                (self.par['p0'] * self.par['p1']) *
                np.log( (self.par['p0'] - self.func_sigmoid_fitted(x)) / self.func_sigmoid_fitted(x)))
    
    def dfdp2(self, x):
        return (self.func_sigmoid_fitted(x) * 
                (self.par['p0'] - self.func_sigmoid_fitted(x)) *
                self.par['p1'] / 
                self.par['p0'])

    # フィット後の関数
    def func_sigmoid_fitted(self, x):
        if not self.isFitted:
            return 0
        return self.func_sigmoid(x, *self.par.values())
    
    # フィット後の関数のエラー
    def func_sigmoid_error(self, x):
        if not self.isFitted:
            return 0
        
        # (df/dp0)^2 * dp0^2
        df2 =  np.power(self.dfdp0(x), 2) * self.cov[0][0]
        
        # (df/dp1)^2 * dp1^2
        df2 += np.power(self.dfdp1(x), 2) * self.cov[1][1]
        
        # (df/dp2)^2 * dp2^2
        df2 += np.power(self.dfdp2(x), 2) * self.cov[2][2]
        
        # 2 * df/dp0 * df/dp1 * (dp0dp1)^2
        df2 += 2.0 * self.dfdp0(x) * self.dfdp1(x) * self.cov[0][1]
        
        # 2 * df/dp1 * df/dp2 * (dp1dp2)^2
        df2 += 2.0 * self.dfdp1(x) * self.dfdp2(x) * self.cov[1][2]
        
        # 2 * df/dp0 * df/dp2 * (dp0dp2)^2
        df2 += 2.0 * self.dfdp0(x) * self.dfdp2(x) * self.cov[0][2]
        
        return np.sqrt(df2)
    
    def func_sigmoid_error_noCov(self, x):
        if not self.isFitted:
            return 0
        
        # (df/dp0)^2 * dp0^2
        df2 =  np.power(self.dfdp0(x), 2) * np.power(self.cov[0][0], 2)
        
        # (df/dp1)^2 * dp1^2
        df2 += np.power(self.dfdp1(x), 2) * np.power(self.cov[1][1], 2)
        
        # (df/dp2)^2 * dp2^2
        df2 += np.power(self.dfdp2(x), 2) * np.power(self.cov[2][2], 2)
                
        return np.sqrt(df2)
```

データをテキトーに生成するための関数を用意しておきます。上記のSigmoid関数のパラメータがそれぞれ

```math
\begin{align}
p_{0} &= 30 \\
p_{1} &= 0.4 \\
p_{2} &= 15
\end{align}
```

ということを想定しており、その上で乱数でデータ点をばらつかせています。

```python
# yのデータ生成用の関数
def ans_y(x):
    ey = np.random.normal(size=len(x))
    y = 30.0/ (1.0 + np.exp((x-15)*(-0.4)))
    return y + y * 0.1 *ey
```

## データ生成

簡単ですが、データを作成しておきます。

```python
# データ準備
x = np.linspace(0, 19, 20) 
y = ans_y(x)
```

データ点を散布図として図示すると下図のようになります。

```python
plt.plot(x, y, 'o')
plt.show()
```


![output_4_0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/329888/7ecffc45-50a3-f356-63f5-6d9435531d11.png)

## フィット実行

さっそくフィットを実行します。初期パラメータ```p0```とパラメータ境界```bounds```を定義して渡しておきます。パラメータ境界については、全パラメータが0以上の値を取るとしています。

```python
# sigmaパラメータなし
fitter = Fitter()

p0 =[1., 1., 1.] 
bounds = ((0.0, 0.0, 0.0), (np.inf, np.inf, np.inf))

fitter.fit(xdata=x, ydata=y, p0=p0, bounds=bounds)
```

    params:  {'p0': 24.313812849093022, 'p1': 13.4398153673723, 'p2': 0.5288927431932293}
    cov. m:  [[ 0.73058063  0.15976218 -0.02689416]
     [ 0.15976218  0.0458503  -0.00592103]
     [-0.02689416 -0.00592103  0.00168229]]
    error :  [0.8547401  0.21412682 0.04101571]

## 誤差伝搬

上記のフィットで3つのパラメータの値と、共分散が求められたので、関数のエラーを計算します。
3つのパラメータがあるので、関数$f$の誤差$\Delta f$は次のような形で表せます。

```math
\begin{align}
(\Delta f)^{2} &=
\Bigl(\frac{\partial f}{\partial p_{0}}\Bigr)^{2} (\Delta p_{0})^{2} +
\Bigl(\frac{\partial f}{\partial p_{1}}\Bigr)^{2} (\Delta p_{1})^{2} +
\Bigl(\frac{\partial f}{\partial p_{2}}\Bigr)^{2} (\Delta p_{2})^{2} + 
2 \frac{\partial f}{\partial p_{0}} \frac{\partial f}{\partial p_{1}} \Delta p_{0} \Delta p_{1} +
2 \frac{\partial f}{\partial p_{0}} \frac{\partial f}{\partial p_{2}} \Delta p_{0} \Delta p_{2} +
2 \frac{\partial f}{\partial p_{1}} \frac{\partial f}{\partial p_{2}} \Delta p_{1} \Delta p_{2} \\
&=
\Bigl(\frac{\partial f}{\partial p_{0}}\Bigr)^{2} \sigma_{p_{0}p_{0}}^{2} +
\Bigl(\frac{\partial f}{\partial p_{1}}\Bigr)^{2} \sigma_{p_{1}p_{1}}^{2} +
\Bigl(\frac{\partial f}{\partial p_{2}}\Bigr)^{2} \sigma_{p_{2}p_{2}}^{2} + 
2 \frac{\partial f}{\partial p_{0}} \frac{\partial f}{\partial p_{1}} \sigma_{p_{0}p_{1}}^{2} +
2 \frac{\partial f}{\partial p_{0}} \frac{\partial f}{\partial p_{2}} \sigma_{p_{0}p_{2}}^{2} +
2 \frac{\partial f}{\partial p_{1}} \frac{\partial f}{\partial p_{2}} \sigma_{p_{1}p_{2}}^{2}
\end{align}
```

ここで、$\sigma^{2}$は共分散、各パラメータについての偏微分については下のように計算しておき、上の式に代入すれば、関数$f$の誤差を求めることができます。
（微分など間違っていたら教えてください）

```math
\begin{align}
\frac{\partial f}{\partial p_{0}} &= \frac{1}{1 + e^{-p_{2}(x-p_{1})}} = \frac{f}{p_{0}} \\

\frac{\partial f}{\partial p_{1}} &= \frac{p_{0}}{1+e^{-p_{2}(x-p_{1})}} e^{-p_{2}(x-p_{1})} p_{2} = \frac{f(p_{0} - f)p_{2}}{p_{0}} \\

\frac{\partial f}{\partial p_{2}} &= \frac{p_{0}}{1+e^{-p_{2}(x-p_{1})}} e^{-p_{2}(x-p_{1})} (-(x - p_{1})) = \frac{f(p_{0} - f)}{p_{0} p_{2}} \log (\frac{a-f}{f})
\end{align}
```

フィット後の関数と、関数のエラーを図示すると下図のようになります。

```python
# sigmaパラメータなし & 共分散あり
pseudox = np.linspace(0, 49, 50) 

res = fitter.func_sigmoid_fitted(pseudox)
err = fitter.func_sigmoid_error(pseudox)

plt.plot(x, y, 'o')
plt.plot(pseudox, res, 'r-',label='fit: p0=%.2f, p1=%.2f, p2=%.2f' % tuple(fitter.par.values()))
plt.fill_between(pseudox, res - err, res + err, facecolor='r', alpha=0.2)
plt.legend()
plt.xlim(0, 50)
plt.grid()
plt.show()
```


![output_6_0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/329888/cc819a7b-f04a-81c9-3d29-9d264b955d08.png)

## curve_fitのsigmaパラメータ

[公式ドキュメント](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html)によると、sigmaパラメータはデータの$y$方向の誤差（あるいは不定性）を表すパラメータで、が1次元の場合、```curve_fit``` の内部でフィット時の$\chi^{2}$を

```math
\chi^{2} = \sum \Bigl( \frac{y_{obs} - f(x)}{sigma} \Bigr)^{2}
```

このような形で求めているとのこと。デフォルトだと$sigma = 1$として計算されるので、各データの$y$方向に誤差がある場合はsigmaパラメータに値を渡しておかないと、正しいフィット結果を得られないです。
ここでは仮に$y$方向の不定性が$\sqrt(y)$程度あるとしてsigmaパラメータに渡してフィットしてみます。

```python
# sigmaパラメータあり
fitter = Fitter()

p0 =[1., 1., 1.] 
bounds = ((0.0, 0.0, 0.0), (np.inf, np.inf, np.inf))
sigma = np.sqrt(y)

fitter.fit(xdata=x, ydata=y, p0=p0, sigma=sigma, bounds=bounds)
```

    params:  {'p0': 26.023718633299193, 'p1': 13.885172632785503, 'p2': 0.4506783056435821}
    cov. m:  [[ 1.34791009e+00  2.89911256e-01 -1.58552587e-02]
     [ 2.89911256e-01  7.53469527e-02 -4.30194913e-03]
     [-1.58552587e-02 -4.30194913e-03  3.70846221e-04]]
    error :  [1.1609953  0.27449399 0.01925737]



```python
# sigmaパラメータあり & 共分散あり
res = fitter.func_sigmoid_fitted(pseudox)
err = fitter.func_sigmoid_error(pseudox)

plt.plot(x, y, 'o')
plt.plot(pseudox, res, 'r-',label='fit: p0=%.2f, p1=%.2f, p2=%.2f' % tuple(fitter.par.values()))
plt.fill_between(pseudox, res - err, res + err, facecolor='r', alpha=0.2)
plt.legend()
plt.xlim(0, 50)
plt.grid()
plt.show()
```


![output_8_0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/329888/f1028e13-b5c5-504d-aa9b-e069a4af2ce9.png)

多少フィット結果とエラーが変わっている感じがします。

## 共分散の非対角成分を考慮しない場合

共分散の非対角成分が小さい値の場合は

```math
(\Delta f)^{2} 
\sim
\Bigl(\frac{\partial f}{\partial p_{0}}\Bigr)^{2} \sigma_{p_{0}p_{0}}^{2} +
\Bigl(\frac{\partial f}{\partial p_{1}}\Bigr)^{2} \sigma_{p_{1}p_{1}}^{2} +
\Bigl(\frac{\partial f}{\partial p_{2}}\Bigr)^{2} \sigma_{p_{2}p_{2}}^{2} 
```

と表せるので、試しにこの場合どうなるのか計算してみます。

```python
# sigmaパラメータあり & 共分散なし

res = fitter.func_sigmoid_fitted(pseudox)
err = fitter.func_sigmoid_error_noCov(pseudox)

plt.plot(x, y, 'o')
plt.plot(pseudox, res, 'r-',label='fit: p0=%.2f, p1=%.2f, p2=%.2f' % tuple(fitter.par.values()))
plt.fill_between(pseudox, res - err, res + err, facecolor='r', alpha=0.2)
plt.legend()
plt.xlim(0, 50)
plt.grid()
plt.show()
```


![output_9_0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/329888/163db67c-6f63-4d69-7233-60648212eb37.png)

この場合は誤差が若干大きく見積もられてしまっています。

# まとめ

データ点をSigmoid関数でフィットした結果をエラー付きで表示してみました。
ちょっと偏微分と誤差伝搬は自信ないですが、間違いがあれば教えてください。