はじめに

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

Sigmoid関数は、Wikipediaなどに記載されているように


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

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


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を使用します。
他にも必要なライブラリをインポートしておきます。

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

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

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関数のパラメータがそれぞれ

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

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

# 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

データ生成

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

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

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

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

フィット実行

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

# 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$は次のような形で表せます。

\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$の誤差を求めることができます。
（微分など間違っていたら教えてください）

\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}

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

# 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()

curve_fitのsigmaパラメータ

公式ドキュメントによると、sigmaパラメータはデータの$y$方向の誤差（あるいは不定性）を表すパラメータで、が1次元の場合、curve_fit の内部でフィット時の$\chi^{2}$を

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

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

# 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]

# 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()

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

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

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

(\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}

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

# 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()

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

まとめ

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