pythonでゆるふわカーブフィッティング


この記事は東京大学航空宇宙工学科／専攻 Advent Calendar 2018 - Adventarの14日目のものです。

adventar.org

カーブフィッティングとは、データに最もよく当てはまる曲線を求めることです。得られた実験データに対してモデル曲線をフィッティングし、そのフィッティングパラメータからデータの持つ特性を求めたり、あるいは実験データの無い領域の値を予測するといった用途で用いられることが多く、様々な分野でごく当たり前のように使われています。

航空宇宙アドベントカレンダーということで何を書こうか非常に迷ったのですが、自分や同期達の卒論にも頻繁に登場していた解析手法ということで最終的にフィッティングの話に落ち着きました。カーブフィッティングはその汎用性の高さにも関わらず、理論的な背景と実際用いられているアルゴリズムの実装とをわかりやすくまとめた記事はそれほど多くないと感じています。そのため本記事では、カーブフィッティングに用いられている代表的な手法を、簡単な理論的背景とpythonによる実装とを合わせてゆるゆる解説していこうと思います。

0. はじめに

今回の記事では、時間・スペースの都合上、最小二乗法によるカーブフィッティングのみを扱っています。
それ以外の手法を用いたものについては他の記事を参照してください。

1. 最小二乗法によるカーブフィッティングの一般論

まず、データが説明変数{\displaystyle \textbf{x}^{(1)},...,\textbf{x}^{(n)}}と目的変数{\displaystyle y^{(1)},...,y^{(n)}}の組{\displaystyle (\textbf{x}^{(1)}, y^{(1)}),...,(\textbf{x}^{(n)}, y^{(n)})}の形で得られるとします。ただし、説明変数{\displaystyle \textbf{x}}はk次元ベクトルの形
$$\textbf{x}=\left(
\begin{array}{c}
x_1 \\
x_2 \\
\vdots \\
x_k
\end{array}
\right)$$
で表されているとします。
このデータを何らかの曲線でフィッティングするとは、説明変数{\displaystyle \textbf{x}}, 目的変数{\displaystyle y}の関係{\displaystyle y=g(\textbf{x}) }を求めることであると言えます。そのために、m個のパラメータ{\displaystyle a_{1},...,a_{m}}と説明変数{\displaystyle x_{1},...,x_{k}}を変数に持つ{\displaystyle f(\textbf{x}, \textbf{a})}という関数形を与えてやり、{\displaystyle a_{1},...,a_{m}}を変化させて{\displaystyle f(\textbf{x}, \textbf{a})}が{\displaystyle \textbf{y}}に対して最も当てはまりが良くなるような{\displaystyle a_{1},...,a_{m}}を求めることを考えます。この問題は、最終的に
$$
S(\textbf{a}) = \sum_{i=1}^{n} ||y^{(i)}-f(\textbf{x}^{(i)}, \textbf{a})||^2
$$
が最小となるような{\displaystyle \textbf{a}}を求める問題に帰着され、この
$$
\textbf{a}=\left(
\begin{array}{c}
a_1 \\
a_2 \\
\vdots \\
a_m
\end{array}
\right)
$$
のことをフィッティングパラメータといいます。
以降、この二乗和S{\displaystyle (\textbf{a})}を最小化する問題を、線形の場合、非線形の場合それぞれについて見ていくこととします。

2. 関数形が線形の場合

当てはめたい関数形{\displaystyle f(\textbf{x}, \textbf{a})}が{\displaystyle \textbf{x}}の1次多項式で書ける時、1節で述べた手法は線形最小二乗法といいます。これを3節で述べる非線形最小二乗法の解法を用いて一般的に解くことも可能ですが、線形最小二乗法の場合は、ある条件においてフィッティングパラメータを解析的に求めることができます。
線形最小二乗法では、当てはめたい関数形は
$$
f(\textbf{x}, \textbf{a}) = a_{0} + a_{1}x_{1} + ...... + a_{k}x_{k}
$$
と書けます。
よって
$$
X = \left(
\begin{array}{cccc}
1 & x_1^{(1)} & ... & x_{k}^{(1)} \\
... & ...& ...&...\\
1 & x_1^{(n)} & ... & x_{k}^{(n)} \\
\end{array}
\right)
$$
とすれば、S{\displaystyle (\textbf{a})}の最小化問題は
$$
\|X\textbf{a} - \textbf{y}\|
$$
の最小化問題として書くことができます。
ただし
$$
\textbf{y}=\left(
\begin{array}{c}
y^{(1)} \\
y^{(2)} \\
\vdots \\
y^{(n)}
\end{array}
\right)
$$
と表すとします。
上の式を2乗して変形すると
\begin{align}
\|X\textbf{a} - \textbf{y}\|^2 &= (X\textbf{a} - \textbf{y})^T(X\textbf{a} - \textbf{y}) \\
&= (\textbf{a}^TX^T-\textbf{y}^T)(X\textbf{a}-\textbf{y}) \\
&= \textbf{a}^TX^TX\textbf{a} - 2\textbf{a}^TX^T\textbf{y} + \textbf{y}^T\textbf{y}
\end{align}
となります。
これを{\displaystyle \textbf{x}}で微分する、すなわち各要素で偏微分して勾配ベクトルを求めると*1、{\displaystyle 2X^TX\textbf{a}-2X^T\textbf{y}}となることから、{\displaystyle \|X\textbf{a} - \textbf{y}\|^2}が最小となる必要条件として、次の正規方程式
$$
X^TX\textbf{a} = X^T\textbf{y}
$$
が得られます。
したがって、{\displaystyle X^TX}が正則であれば、フィッティングパラメータ{\displaystyle \textbf{a}}を
$$
\textbf{a} = (X^TX)^{-1}(X^T\textbf{y})
$$
と解析的に求めることができます。

ではこの結果を用いて、実際にデータにフィットする直線を求めてみます。
ここでは有名なIrisのデータを用いて直線フィッティングします。

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd


def load_data():
    """irisデータの読み取り"""
    path = 'https://archive.ics.uci.edu/ml/machine-learning-databases/iris/'
    data = 'iris.data'

    df = pd.read_table(path + data, sep=',', header=None)
    df.drop([4], axis=1, inplace=True)
    return df


def preprocess(df):
    """バイアス項の追加などの前処理"""
    X = df[[2]].values
    y = df[[3]].values

    X = np.insert(X, 0, np.ones((1, X.shape[0]), dtype=int), axis=1)

    return X, y


def normal_equation(X, y):
    """正規方程式を解く"""
    A = np.dot(X.T, X)
    B = np.dot(X.T, y)
    a = np.linalg.solve(A, B)
    return a


def plot(X, y, params):
    """結果のプロット"""
    fig, ax = plt.subplots(figsize=(5, 5))
    fig.subplots_adjust(bottom=0.2)
    fig.subplots_adjust(left=0.2)
    x = np.arange(-0, 8, 0.01)
    best_fit = params[1]*x+params[0]
    ax.set_xlabel('petal_length')
    ax.set_ylabel('petal_width')
    ax.set_xlim([0.0, 8.0])
    ax.set_ylim([0.0, 3.0])
    ax.scatter(X[:, 1],y[:, 0], color='blue', s=5)
    ax.plot(x, best_fit, color='red')
    plt.savefig('result.jpg')
    plt.show()


if __name__ == '__main__':

    df = load_data()
    X, y = preprocess(df)
    params = normal_equation(X, y)
    plot(X, y, params)

f:id:saefro:20181212153047j:plain
irisデータの直線フィッティング結果

3. 関数形が非線形の場合

非線形の場合、また{\displaystyle X^TX}の正則性が成り立たない場合は先ほどのようにフィッティングパラメータを解析的に求めることができません。そのため、反復計算により数値的に求めてやる必要があります。
今回はNewton法に由来するアルゴリズムであり、非線形最小二乗法問題にのみ適用可能なGauss-Newton法、Levenberg-Marquardt法について解説したいと思います。

3-0. 準備 (Newton法について)

Gauss-Newton法、Levenberg-Marquardt法の解説に移る前に、まずNewton法について簡単に述べておきます。
まず関数{\displaystyle f(\textbf{x})}が2次関数である時、以下のような形式で書き表す事が出来ます。
$$
f(\textbf{x}) = \frac{1}{2} \textbf{x}^T H \textbf{x} + c\textbf{x} + c_{0}
$$
この時
$$
\nabla f(\textbf{x}) = H\textbf{x} + c
$$
となることから、停留点は{\displaystyle \textbf{x}^{\ast} = -H^{-1}c}のみであり、Hが正定値であればこれが最小解となります。
したがって2次関数{\displaystyle f(\textbf{x})}に対して、ヘッセ行列Hが正定値であれば最小解を簡単に求めることが出来ます。

Newton法とは、元の関数のテイラー近似と上で述べた性質から最小解を求める方法です。
まず一般の関数{\displaystyle f(\textbf{x})}を下のように2次の項までテイラー展開します。
$$
\tilde{f}(\textbf{x}) = f(\bar{\textbf{x}}) + \nabla f(\bar{\textbf{x}})^T (\textbf{x}-\bar{\textbf{x}}) + \frac{1}{2}(\textbf{x}-\bar{\textbf{x}})^T H (\textbf{x}-\bar{\textbf{x}})
$$
ヘッセ行列Hが正定値のとき、{\displaystyle \tilde{f}(\textbf{x})}の最適値は
\begin{align}
\textbf{x} &= \bar{\textbf{x}} - H^{-1}\nabla f(\bar{\textbf{x}}) \\
&= \bar{\textbf{x}} - H^{-1} \textbf{g}
\end{align}
となります。
したがって、{\displaystyle \bar{\textbf{x}} - H^{-1} \textbf{g}}は関数{\displaystyle f}の最適値の良い近似解になることが期待され、Newton法ではこの考え方を用いて現在の点{\displaystyle \textbf{x}}を{\displaystyle \textbf{x} - H^{-1} \textbf{g}}に移動させることを繰り返すことにより、最適値を求めます。
今回のようにフィッティングパラメータ{\displaystyle \textbf{a}}を求める場合、以下のような漸化式によってパラメータを更新していきます。
$$
\textbf{a}^{(s+1)} = \textbf{a}^{(s)} - H^{-1} \textbf{g}
$$
Newton法は最急降下法などと比較して解周辺での収束が速いという特徴がありますが、ヘッセ行列とその逆行列の計算負荷が大きい、ヘッセ行列が正定値であることを仮定しているなどの欠点があり、一般の非線形最小二乗法問題ではあまり実用的な方法であるとは言えません。

3-1. Gauss-Newton法

Gauss-Newton法は先述のNewton法を改良した手法です。まずSの定義から、勾配{\displaystyle \textbf{g}}は以下の形で表されます。
$$
\textbf{g} = 2J_{r}^{T} \textbf{r}
$$
ヘッセ行列はこの勾配{\displaystyle \textbf{g}}を{\displaystyle \textbf{a}}で微分することによって定義されるので
$$
H_{jk} = 2 \sum_{i=1}^{m}(\frac{\partial r_{i}}{\partial a_{j}} \frac{\partial r_{i}}{\partial a_{k}}
+
r_{i} \frac{\partial^2 r_{i}}{\partial a_{j} \partial a_{k}})
$$
この右辺第2項を無視することにより
$$
H \approx 2J_{r}^{T}J_{r}
$$
が得られます。
よって、これらをNewton法の漸化式に代入することにより
$$
\textbf{a}^{(s+1)} = \textbf{a}^{(s)} - (J_{r}^{T}J_{r})^{-1} J_{r}^{T} \textbf{r}
$$
という漸化式が得られます。
Gauss-Newton法では、この漸化式を用いて、初期値{\displaystyle \textbf{a}^{(0)}}から始めてパラメータを更新していきます。

ではGauss-Newton法を用いて実際に非線形カーブフィッティングを行ってみます。
まずは簡単なwikipediaにある反応速度の例からやってみます。

i 1 2 3 4 5 6 7
[S] 0.038 0.194 0.425 0.626 1.253 2.500 3.740
v 0.050 0.127 0.094 0.2122 0.2729 0.2665 0.3317

上の表のようなデータをミカエリス・メンテンの式より得られる曲線
$$
v = \frac{V_{\rm{max}}[S]}{K_{M}+[S]}
$$
でフィッティングしたいと思います。

import numpy as np
import matplotlib.pyplot as plt
import sympy as sym


def michaelis(x, a, b):
    """ミカエリス・メンテンの式"""
    f = a*x / (x + b)
    return f


def jacobi_res(p1, p2):
    """ヤコビ行列と残差ベクトルの計算"""
    jacobi = []
    res = []
    for i in range(arr_S.shape[0]):
        f = michaelis(arr_S[i], a, b)
        r = arr_v[i] - f
        j1 = sym.diff(r, a)
        j2 = sym.diff(r, b)
        j1 = j1.subs([(a, p1), (b, p2)])
        j2 = j2.subs([(a, p1), (b, p2)])
        j = [j1, j2]
        r = r.subs([(a, p1), (b, p2)])
        jacobi.append(j)
        res.append(r)
    jacobi = np.array(jacobi).astype('float')
    res = np.array(res).astype('float')
    return jacobi, res


def initialize(init_a, init_b):
    """パラメータベクトルの初期化"""
    params = np.array([init_a, init_b])
    return params


def gauss_newton(params):
    """最大反復回数20回として反復計算"""
    n_iteration = 0
    res_list = []
    while n_iteration <= 20:
        J, r = jacobi_res(params[0], params[1])
        inv = np.linalg.inv(np.dot(J.T, J))
        new_params = params - np.dot(np.dot(inv, J.T), r)
        if np.linalg.norm(new_params - params)/np.linalg.norm(params) < 1.0*10**(-15):
            break
        params = new_params
        res_list.append(np.linalg.norm(r)**2)
        n_iteration += 1
    return params, res_list


def plot(params, res_list):
    """グラフの描画"""
    fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(16, 6))
    fig.subplots_adjust(wspace=0.5, hspace=0.2)

    x1 = np.arange(0, 4, 0.01)
    func = params[0] * x1 / (x1 + params[1])
    ax1.set_xlabel('[s]', fontdict={'fontsize':18})
    ax1.set_ylabel('v', fontdict={'fontsize':18})
    ax1.plot(x1, func, color='black')
    ax1.scatter(arr_S, arr_v, color='red')

    x2 = np.arange(0, len(res_list), 1)
    ax2.set_yscale('log')
    ax2.set_xlabel('number of iteration', fontdict={'fontsize':18})
    ax2.set_ylabel('S(a)', fontdict={'fontsize':18})
    ax2.plot(x2, np.array(res_list), color='black')
    ax2.set_ylim([1.0*10**(-3), 1.0])
    ax2.set_xlim([0, 20])
    ax2.set_xticks(np.arange(20))

    plt.savefig('gauss-newton.jpg')
    plt.show()


if __name__ == '__main__':

    arr_S = np.array([0.038, 0.194, 0.425, 0.626, 1.253, 2.500, 3.740])
    arr_v = np.array([0.050, 0.127, 0.094, 0.2122, 0.2729, 0.2665, 0.3317])

    (x, a, b) = sym.symbols('x a b')

    init_params = initialize(0.9, 0.2)
    params, res_list = gauss_newton(init_params)
    plot(params, res_list)

f:id:saefro:20181213154522j:plain
フィッティング結果(左)とS(a)の変化(右)


もう一つ例をやってみます。
{\displaystyle f(x) = \frac{a}{\sqrt{2\pi}\sigma} \exp{\frac{(x-\mu)^2}{2\sigma^2}}}の式で表されるガウス関数(期待値{\displaystyle \mu=150}, 標準偏差{\displaystyle \sigma=10}, {\displaystyle a=30})に[-0.05, 0.05]の一様分布に従う誤差項を付与したデータを生成し、これをガウス関数でフィッティングする事を考えます (ちなみにこのデータを生成するコードはLevenberg-Marquardt法のところで載せるものと同じなのでそっちを参照してください)。
初期値を{\displaystyle \mu^{(0)}=150}, {\displaystyle \sigma^{(0)}=5}, {\displaystyle a^{(0)}=5}として与え、先ほどと同様に計算すると、以下のようにフィッティングできます。

f:id:saefro:20181213154216j:plain
Gauss-Newton法によるガウシアンフィッティング結果({\displaystyle \mu^{(0)}=150})

先程の例では期待値の初期値として正確な値を与えましたが、これを変化させて、初期値を{\displaystyle \mu^{(0)}=130}, {\displaystyle \sigma^{(0)}=5}, {\displaystyle a^{(0)}=5}として与えてみます。

f:id:saefro:20181213155322j:plain
Gauss-Newton法によるガウシアンフィッティング結果({\displaystyle \mu^{(0)}=130})
完全に発散していますね.......

このようにGauss-Newton法でガウシアンフィッティングを行う場合、期待値の初期値として正確な値を与える必要があります。

Gauss-Newton法では、下式のような2階微分項を無視するという近似を行っています。
$$
\|r_{i} \frac{\partial^2 r_{i}}{\partial a_{j} \partial a_{k}}\| << \|\frac{\partial r_{i}}{\partial a_{j}} \frac{\partial r_{i}}{\partial a_{k}}\|
$$
この条件が正しく成り立つのは次の2つのような場合です。

  • {\displaystyle r_{i}}が十分小さい。
  • {\displaystyle\|r_{i} \frac{\partial^2 r_{i}}{\partial a_{j} \partial a_{k}}\|}が小さい、即ち関数の非線形性が小さい。

Gauss-Newton法はこのような2つの条件が満たされない場合には安定性が悪く、また残差平方和Sが反復ごとに必ずしも減少するわけではないことから、実用的にはなんらかの安定化を行う必要があります。先ほどのガウシアンフィッティングの例だと、初期値として正確な値を与えないと初期の時点での{\displaystyle r_{i}}が大きすぎて発散してしまうということになります。次節では、このGauss-Newton法の改善バージョンの一つであるLevenberg-Marquardt法を紹介します。

3-2. Levenberg-Marquardt法

Levenberg-Marquardt法はGauss-Newton法の安定性を向上させた手法です。二乗誤差の最小化の問題では頻繁に用いられる実用的な手法であり、フィッティングモジュールのscipy.optimize.curve_fitでも、拘束条件の無い問題ではこのLevenberg-Marquardt法がデフォルトで適用されます。
Levenberg-Marquardt法では、以下の漸化式によってパラメータを更新します。
$$
\textbf{a}^{(s+1)} = \textbf{a}^{(s)} - (J_{r}^T J_{J} + \lambda D)^{-1} J_{r}^T \textbf{r}
$$
ここでDは正定値対角行列です。
この正定値対角行列Dには、最も簡単な場合だと単位行列Iを用いることもありますが、より収束性をよくするためには{\displaystyle J_{r}^TJ_{r}}の対角成分のみを並べ、それ以外の成分は0であるような行列を用います。また、発散を回避するために縮小因子{\displaystyle \alpha}*2を導入すると、最終的な漸化式は以下のように書けます。
$$
\textbf{a}^{(s+1)} = \textbf{a}^{(s)} - \alpha(J_{r}^T J_{J} + \lambda \rm{diag}(J_{r}^TJ_{r}))^{-1} J_{r}^T \textbf{r}
$$
適度に{\displaystyle \lambda}を変化させながらパラメータを更新することで、発散を回避しながら最適値に到達します。
直感的に言えば、解が近くにない時は{\displaystyle \lambda}を大きくして最急降下方向に更新し、解の近くでは{\displaystyle \lambda}を小さくし、Gauss-Newton法で一気に収束させるようなイメージです。そのため、Levenberg-Marquardt法は最急降下法とGauss-Newton法を組み合わせた手法であると言えます。

ではこのLevenberg-Marquardt法を用いて先程のガウシアンフィッティングをやってみます。
{\displaystyle S(\textbf{a})}の値が前回のステップでの値より小さい場合は、更新がうまくいったということで{\displaystyle \lambda}を1/10する事によってよりGauss-Newton方向に修正し、{\displaystyle S(\textbf{a})}の値が前回のステップでの値より大きい場合は、更新がうまくいかなかったという事で{\displaystyle \lambda}を10倍し、最急降下方向に修正します。また、始めは縮小因子{\displaystyle \alpha=1.0}としてスタートし、パラメータのいずれかの絶対値が与えられたthreshold(今回は1000)より大きくなった場合、{\displaystyle \alpha}を1/10して計算をやり直すこととしました。
ソースコードは以下の通りです。

import numpy as np
import matplotlib.pyplot as plt
import sympy as sym


def gaussian(arr_x, mean, sigma, a):
    """ガウス関数 (引数に値をとる)"""
    return a/(np.sqrt(2*np.pi)*sigma) * np.exp(-(arr_x - mean)**2 / (2*sigma**2))


def fitting_func(x, mean, sigma, amp):
    """ガウス関数 (引数にsympyのsumbolを取る)"""
    f = amp / (np.sqrt(2*np.pi)*sigma) * sym.exp(-(x - mean)**2 / (2*sigma**2))
    return f


def generate_data():
    """フィッティング用データの生成"""
    rand = np.random.rand(100)
    rand = (2 * rand - 1) * 0.05
    arr_x = np.arange(100, 200, 1)
    arr_y = gaussian(arr_x, 150, 10, 30) + rand
    return arr_x, arr_y


def jacobi_res(p1, p2, p3):
    """ヤコビ行列と残差ベクトルの計算"""
    jacobi = []
    res = []
    for i in range(arr_x.shape[0]):
        f = fitting_func(arr_x[i], a, b, c)
        r = arr_y[i] - f
        j1 = sym.diff(r, a)
        j2 = sym.diff(r, b)
        j3 = sym.diff(r, c)
        j1 = j1.subs([(a, p1), (b, p2), (c, p3)])
        j2 = j2.subs([(a, p1), (b, p2), (c, p3)])
        j3 = j3.subs([(a, p1), (b, p2), (c, p3)])
        j = [j1, j2, j3]
        r = r.subs([(a, p1), (b, p2), (c, p3)])
        jacobi.append(j)
        res.append(r)
    jacobi = np.array(jacobi).astype('float')
    res = np.array(res).astype('float')
    return jacobi, res


def initialize(init_a, init_b, init_c):
    """パラメータベクトルの初期化"""
    params = np.array([init_a, init_b, init_c])
    return params


def levenberg(params):
    """最大反復回数100回として反復計算"""
    n_iteration = 0
    res2_list = []
    lam = 0.01
    alpha = 1.0
    res_old = np.inf
    threshold = 1000
    init_params = params
    while n_iteration <= 100:
        J, r = jacobi_res(params[0], params[1], params[2])
        JJ = np.dot(J.T, J)
        inv = np.linalg.inv(JJ + lam*np.diag(np.diag(JJ)))
        new_params = params - alpha*np.dot(np.dot(inv, J.T), r)
        determination = np.linalg.norm(new_params - params)/np.linalg.norm(params) 
        if determination < 1.0*10**(-5):
            break
        if np.linalg.norm(r) < res_old:
            lam = lam/10
        else:
            lam = lam*10
        params = new_params
        res_old = np.linalg.norm(r)
        res2_list.append(np.linalg.norm(r)**2)
        n_iteration += 1
        if abs(params[0]) > threshold or abs(params[1]) > threshold or abs(params[2]) > threshold:
            alpha = alpha / 10
            n_iteration = 0
            params = init_params
            res2_list = []
    return params, res2_list


def plot(params, res2_list):
    """グラフの描画"""
    fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(16, 6))
    fig.subplots_adjust(bottom=0.2)
    fig.subplots_adjust(left=0.2)

    x1 = np.arange(100, 200, 1)
    func = gaussian(x1, params[0], params[1], params[2])
    ax1.set_xlabel('x', fontdict={'fontsize':18})
    ax1.set_ylabel('y', fontdict={'fontsize':18})
    ax1.plot(x1, func, color='red')
    ax1.scatter(arr_x, arr_y, color='blue')

    x2 = np.arange(0, len(res2_list), 1)
    ax2.set_yscale('log')
    ax2.set_xlabel('number of iteration', fontdict={'fontsize':18})
    ax2.set_ylabel('S(a)', fontdict={'fontsize':18})
    ax2.plot(x2, np.array(res2_list), color='black')
    ax2.set_ylim([1.0*10**(-3), 100.0])
    ax2.set_xlim([0, 200])

    plt.savefig('levenberg.jpg')
    plt.show()


if __name__ == '__main__':

    arr_x, arr_y = generate_data()

    (x, a, b, c) = sym.symbols('x a b c')

    init_params = initialize(130, 5, 5)
    params, res2_list = levenberg(init_params)
    plot(params, res2_list)

f:id:saefro:20181213164503j:plain
Levenberg-Marquardt法によるガウシアンフィッティング結果({\displaystyle \mu^{(0)}=130})
Gauss-Newton法で発散した条件でも、このように改良を加える事によってフィッティングできていることが分かります。

4. ライブラリの紹介

最後に、pythonで使えるフィッティングモジュールを3つ紹介したいと思います。

4-1. numpy.polyfit

多項式でカーブフィッティングしたい時には、numpyのpolyfitが簡単で便利です。
データ(x, y)を2次関数でフィッティングしたい場合、フィッティングパラメータを1行のコードで求めることが出来ます。

params = np.polyfit(x, y, 2)

params[0]には最高次の項の係数が入っています。

f:id:saefro:20181214113936j:plain
numpy.polyfitによるフィッティング結果

4-2. scipy.optimize.curve_fit

scipy.optimize.curve_fitでは、自分で定義した関数を用いて簡単にフィッティングが行えます。
自分で定義した関数をscipy.optimize.curve_fitの第一引数に渡すと、最適なフィッティングパラメータとその共分散行列*3がタプルで返ってきます。
先ほど扱ったガウシアンフィッティングの例でやってみましょう。

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


def gaussian(arr_x, mean, sigma, a):
    """ガウス関数 (引数に値をとる)"""
    return a/(np.sqrt(2*np.pi)*sigma) * np.exp(-(arr_x - mean)**2 / (2*sigma**2))


def generate_data():
    """フィッティング用データの生成"""
    rand = np.random.rand(100)
    rand = (2 * rand - 1) * 0.05
    arr_x = np.arange(100, 200, 1)
    arr_y = gaussian(arr_x, 150, 10, 30) + rand
    return arr_x, arr_y


def initialize(init_a, init_b, init_c):
    """パラメータベクトルの初期化"""
    params = np.array([init_a, init_b, init_c])
    return params


def plot(params):
    """グラフの描画"""
    fig, ax = plt.subplots(figsize=(6, 4))
    fig.subplots_adjust(bottom=0.2)
    fig.subplots_adjust(left=0.2)

    x1 = np.arange(100, 200, 1)
    func = gaussian(x1, params[0], params[1], params[2])
    ax.set_xlabel('x', fontdict={'fontsize':18})
    ax.set_ylabel('y', fontdict={'fontsize':18})
    ax.plot(x1, func, color='red')
    ax.scatter(arr_x, arr_y, color='blue')

    plt.savefig('levenberg.jpg')
    plt.show()


if __name__ == '__main__':

    arr_x, arr_y = generate_data()
    init_params = initialize(150, 5, 5)
    opt_params, cov = scipy.optimize.curve_fit(gaussian, arr_x, arr_y, p0=init_params)
    plot(opt_params)

f:id:saefro:20181214013130j:plain
scipy.optimize.curve_fitによるフィッティング結果

4-3. lmfit

,最後にlmfitを紹介します。
lmfitはscipy.optimize.curve_fitのように自分で定義した関数を用いてフィッティングすることもできますが、主要な関数形のクラスがあらかじめ定義されているので、それらを呼び出すだけで良いので関数形を自分で定義する必要がない場合には非常に便利です。また、フィッティングパラメータに制限をかけるなどの操作が簡単に行えるというメリットもあります*4。
ではこのlmfitを用いてNISTのスペクトルデータ*5のフィッティングを行ってみます。
このデータは exp関数の上にガウシアンが2つ乗っているような比較的複雑な形をしていますが、 lmfitを使えば簡単にフィッティングが行えます。

import pandas as pd
import matplotlib.pyplot as plt

"""データの読み取り"""
data = pd.read_csv('gauss_data.csv', header=None)

x = data[1]
y = data[0]

"""フィッティングしたいモデルのインポート"""
from lmfit.models import ExponentialModel, GaussianModel

"""expモデルの定義"""
# パラメータオブジェクトparsの生成
exp_mod = ExponentialModel(prefix='exp_')
pars = exp_mod.guess(y, x=x)

# 1つ目のピーク
gauss1 = GaussianModel(prefix='g1_')
pars.update(gauss1.make_params())
pars['g1_center'].set(55, min=40, max=100)
pars['g1_sigma'].set(10, min=3)
pars['g1_amplitude'].set(2000, min=10)

# 2つ目のピーク
gauss2 = GaussianModel(prefix='g2_')
pars.update(gauss2.make_params())
pars['g2_center'].set(170, min=150, max=200)
pars['g2_sigma'].set(10, min=3)
pars['g2_amplitude'].set(2000, min=10)

# モデルの合成
mod = gauss1 + gauss2 + exp_mod

# 初期値
init = mod.eval(pars, x=x)

# 最適値
out = mod.fit(y, pars, x=x)

fig, ax = plt.subplots(figsize=(4, 3))
fig.subplots_adjust(bottom=0.2)
fig.subplots_adjust(left=0.2)

ax.scatter(x, y, s=5, color='blue')
ax.plot(x, out.best_fit, 'r-')
ax.set_xlabel('wavelength (nm)')
ax.set_ylabel('intensity (a.u.)')

plt.savefig('lm.jpg')
plt.show()

f:id:saefro:20181214020905j:plain
lmfitによるフィッティング結果

5. 終わりに

これまで、カーブフィッティングの理論を簡単な実装つきで解説し、実用的なライブラリの紹介も行ってきました。結果的にざっくりとした雰囲気はつかめるような記事にはなったかなあと思いますが、分量の関係で適合度や誤差の話ができず、フィッティング結果を定量的に評価できているとは言えないので、この記事の続編を書くとしたらそのあたりの話をしようかなと思います(記事のタイトルに「ゆるふわ」とつけたのは主にそういう理由です)。また、至らないところばかり*6だと思うので、間違いなどがあれば気軽に教えていただけるとありがたいです。それでは、駄文、長文に付き合っていただきありがとうございました。

出典

  1. 東京大学工学教程 最適化と変分法
  2. 非線形最小二乗法 - Wikipedia
  3. ガウス・ニュートン法 - Wikipedia
  4. Levenberg–Marquardt algorithm - Wikipedia
  5. 制約なし非線形最適化アルゴリズム - MATLAB & Simulink - MathWorks 日本
  6. numpy.polyfit — NumPy v1.16 Manual
  7. scipy.optimize.curve_fit — SciPy v1.3.0 Reference Guide
  8. Non-Linear Least-Squares Minimization and Curve-Fitting for Python — Non-Linear Least-Squares Minimization and Curve-Fitting for Python

*1:「ベクトルで微分・行列で微分」公式まとめ - Qiitaが参考になります。

*2:学習率とも言いますかね。

*3:次回は共分散行列を求めて誤差評価する話がしたいですね...

*4:これについては時間・スペース的な都合であんま解説できません...申し訳ないです

*5:StRD Nonlinear Least Squares Regression Datasets

*6:線型方程式を解くときに逆行列をそのまま計算していたり、収束判定が適当(&恣意的)だったりと計算の部分がガバガバです...続編を書くなら数値計算的にもっとマシな実装にしたいですね