peak deconvolution の覚え書き①～5 点を通る三次式～

5 点に三次関数を当てはめる．
すなわち，5 点 $${(x_1, y_1), (x_2, y_2), (x_3, y_3), (x_4, y_4), (x_5, y_5)}$$ があるとき，もっともよく当てはまるような $${y=ax^3+bx^2+cx+d}$$ を与える定数 $${a, b, c, d}$$ の組を考える．

理屈

単回帰分析と同様に残差平方和を最小とすることを考えると，$${y_k'=ax_k^3+bx_k^2+cx_k+d}$$ として$${S=\displaystyle\sum_{k=1}^n(y_k-y_k')^2}$$ を最小化する (今回は $${n=5}$$)．

$${S}$$ について展開すると，
$${S=\displaystyle\sum_{k=1}^n(y_k-y_k')^2=\displaystyle\sum_{k=1}^n(y_k-ax_k^3-bx_k^2-cx_k-d)^2}$$

すなわち $${S}$$ は $${a, b, c, d}$$ について二次式なので最小値において
$${\cfrac{\partial S}{\partial a}=0, \cfrac{\partial S}{\partial b}=0, \cfrac{\partial S}{\partial c}=0, \cfrac{\partial S}{\partial d}=0}$$

$${\cfrac{\partial S}{\partial a}=\displaystyle\sum_{k=1}^n(-2x_k^3(y_k-ax_k^3-bx_k^2-cx_k-d))}$$
$${=2a\displaystyle\sum_{k=1}^nx_k^6+2b\displaystyle\sum_{k=1}^nx_k^5+2c\displaystyle\sum_{k=1}^nx_k^4+2d\displaystyle\sum_{k=1}nx_k^3-2\displaystyle\sum_{k=1}^nx_k^3y_k}$$

$${\cfrac{\partial S}{\partial b}=\displaystyle\sum_{k=1}^n(-2x_k^2(y_k-ax_k^3-bx_k^2-cx_k-d))}$$
$${=2a\displaystyle\sum_{k=1}^nx_k^5+2b\displaystyle\sum_{k=1}^nx_k^4+2c\displaystyle\sum_{k=1}^nx_k^3+2d\displaystyle\sum_{k=1}^nx_k^2-2\displaystyle\sum_{k=1}^nx_k^2y_k}$$

$${\cfrac{\partial S}{\partial c}=\displaystyle\sum_{k=1}^n(-2x_k(y_k-ax_k^3-bx_k^2-cx_k-d))}$$
$${=2a\displaystyle\sum_{k=1}^nx_k^4+2b\displaystyle\sum_{k=1}^nx_k^3+2c\displaystyle\sum_{k=1}^nx_k^2+2d\displaystyle\sum_{k=1}^nx_k-2\displaystyle\sum_{k=1}^nx_ky_k}$$

$${\cfrac{\partial S}{\partial d}=\displaystyle\sum_{k=1}^n(-2(y_k-ax_k^3-bx_k^2-cx_k-d))}$$
$${=2a\displaystyle\sum_{k=1}^nx_k^3+2b\displaystyle\sum_{k=1}^nx_k^2+2c\displaystyle\sum_{k=1}^nx_k+2d\displaystyle\sum_{k=1}^n1-2\displaystyle\sum_{k=1}^ny_k}$$

よって
$${\cfrac{\partial S}{\partial a}=0, \cfrac{\partial S}{\partial b}=0, \cfrac{\partial S}{\partial c}=0, \cfrac{\partial S}{\partial d}=0}$$
$${\iff \begin{pmatrix} \displaystyle\sum_{k=1}^nx_k^6&\displaystyle\sum_{k=1}^nx_k^5&\displaystyle\sum_{k=1}^nx_k^4&\displaystyle\sum_{k=1}^nx_k^3 \\\displaystyle\sum_{k=1}^nx_k^5&\displaystyle\sum_{k=1}^nx_k^4&\displaystyle\sum_{k=1}^nx_k^3&\displaystyle\sum_{k=1}^nx_k^2 \\\displaystyle\sum_{k=1}^nx_k^4&\displaystyle\sum_{k=1}^nx_k^3&\displaystyle\sum_{k=1}^nx_k^2&\displaystyle\sum_{k=1}^nx_k \\\displaystyle\sum_{k=1}^nx_k^3&\displaystyle\sum_{k=1}^nx_k^2&\displaystyle\sum_{k=1}^nx_k^1&\displaystyle\sum_{k=1}^n1 \end{pmatrix}\begin{pmatrix} a \\b \\c \\d \end{pmatrix}=\begin{pmatrix} \displaystyle\sum_{k=1}^nx_k^3y_k \\\displaystyle\sum_{k=1}^nx_k^2y_k \\\displaystyle\sum_{k=1}^nx_ky_k \\\displaystyle\sum_{k=1}^ny_k \end{pmatrix}}$$

実践

では実践してみよう．
$${x=-1.0, -0.5, 0.0, 0.5, 1.0}$$ の 4 点で 3 回ずつデータを測定したら Fig. 1 のような結果が得られたとする．

Fig. 1 等間隔の $${x}$$ 5 点で各 3 回データを取得した結果

なおこれは $${y=x(x-1)(x+1)}$$ の各結果に平均 0，分散 0.05 のガウス白色ノイズを付加してつくったダミーデータである．つまり「本当の」結果は Fig. 2 青点線のようになっているはずである．

Fig. 2 $${y=x^3-x}$$ を Fig. 1 に重ねたもの．

この時の $${(x, y)}$$ の組合せは Table 1のとおりである．

Table 1

これを前項目で示した結果に代入すると最も当てはまりのいい三次関数を与える $${(a, b, c, d)}$$ の組は Fig. 3 の行列を解いて求まる．

Fig. 3

4 行 4 列の行列の逆行列なので手計算だといささか大変だが，Excel の場合は逆行列を求める MINVERSE 関数と行列の積を求める MMULT 関数があるので，これを組み合わせれば簡単に求まる．

Fig. 3 より $${y\approx1.1x^3+0.0x^2-1.1x-0.0}$$ であり，青点線の $${y=x(x-1)(x+1)=x^3-x}$$ とほぼ同じ結果になっていることがわかる．

Fig. 4 実測値 (赤丸) に基づく予測値 (赤点線) と理論値 (青点線)

ちなみに

「Excel の近似式って直線以外もできなかったっけ？」と思ったあなたは非常に正しいです．
近似のグラフも近似式も $${R^2}$$ 値もワンクリックで出せます．
行列を二つも書いて計算する意味はほぼないです．
ほんとは Savitzky-Golay 法に応用したかったんだけどそこまでたどり着けなかったのでまたいつか．
(2024 年 1 月 1 日：誤字等修正しました)