あいさつ

突然ですが，数式を使えるようにしてみたので，テスト的な記事です．
自分は Hexo という静的サイトジェネレータを使っています．「Hexo 数式」とかでググると mathjax をインストールして…みたいな記事が多くヒットしますが，僕の環境ではなぜか動きませんでした．結局，ここに書いているscriptタグを記事冒頭に埋め込むといい感じに動いてそうだったので，それでテストしています．

本題

さて，スペクトルなどをガウス関数でフィッティングしたいということはよくあるわけですが，単純に一つのガウス関数だけなら簡単にフィッティングする方法がいくつかあります．そのなかで，非線形最小二乗法について書きます．これは Caruana’s Algorithm とよばれるもので，昔参考にしていた Qiita の記事があったのにそれが見つからなくなってしまったので，ここに書いておくことにします．
観測したデータを，次のガウス関数でフィッティングしたいとします．

$$f(x) = {A \exp \frac{-(x - \mu)^2}{2 \sigma^2}}$$

両辺の対数をとると，

$$\ln f(x) = \ln A + \frac{-(x - \mu)^2}{2 \sigma^2} = \ln A - \frac{\mu^2}{2 \sigma^2} + \frac{\mu x}{\sigma^2} - \frac{x^2}{2 \sigma^2} = a + bx + cx^2$$

ここで，$a = \ln A - \frac{\mu^2}{2 \sigma^2}$，$b = \frac{\mu}{\sigma^2}$，$c = - \frac{1}{2 \sigma^2}$とおきました．
観測したデータを$g(x)$とすると，フィッティング問題は，次の誤差関数$S$を最小にする組$(a, b, c)$を探す問題に帰着します：

$${S(a, b, c) = \sum_{x} \{ \ln g(x) - (a + bx + cx^2) \}^2}$$

次に，${S}$を$a, b, c$で偏微分して，極値となる場所を考えます．

\begin{eqnarray}\frac{\partial S}{\partial a} &=& 2 \sum_{x} \{ \ln g(x) - (a + bx + cx^2) \}(-1) = 0 \\ &\Leftrightarrow& a \sum_{x} 1 + b \sum_{x} x + c \sum_{x} x^2 = \sum_{x} \ln g(x)\end{eqnarray} \begin{eqnarray}\frac{\partial S}{\partial b} &=& 2 \sum_{x} \{ \ln g(x) - (a + bx + cx^2) \}(-x) = 0 \\ &\Leftrightarrow& a \sum_{x} x + b \sum_{x} x^2 + c \sum_{x} x^3 = \sum_{x} x \ln g(x)\end{eqnarray} \begin{eqnarray}\frac{\partial S}{\partial c} &=& 2 \sum_{x} \{ \ln g(x) - (a + bx + cx^2) \}(-x^2) = 0 \\ &\Leftrightarrow& a \sum_{x} x^2 + b \sum_{x} x^3 + c \sum_{x} x^4 = \sum_{x} x^2 \ln g(x)\end{eqnarray}

したがって，データの個数を$N = \sum_{x} 1$とすれば，

\begin{equation} \begin{pmatrix} N & \sum_{x} x & \sum_{x} x^2 \\ \sum_{x} x & \sum_{x} x^2 & \sum_{x} x^3 \\ \sum_{x} x^2 & \sum_{x} x^3 & \sum_{x} x^4 \end{pmatrix} \begin{pmatrix} a \\ b \\ c \end{pmatrix} = \begin{pmatrix} \sum_{x} \ln g(x) \\ \sum_{x} x \ln g(x) \\ \sum_{x} x^2 \ln g(x) \end{pmatrix} \end{equation}

と書けます．今，$N$や$x$，$g(x)$は既知ですから，

\begin{equation} \begin{pmatrix} a \\ b \\ c \end{pmatrix} = \begin{pmatrix} N & \sum_{x} x & \sum_{x} x^2 \\ \sum_{x} x & \sum_{x} x^2 & \sum_{x} x^3 \\ \sum_{x} x^2 & \sum_{x} x^3 & \sum_{x} x^4 \end{pmatrix}^{-1} \begin{pmatrix} \sum_{x} \ln g(x) \\ \sum_{x} x \ln g(x) \\ \sum_{x} x^2 \ln g(x) \end{pmatrix} \end{equation}

とすれば，組$(a, b, c)$が求まりますね．
もともとの未知数$A, \mu, \sigma$について解けば，

$$A = \exp\left(a - \frac{b^2}{4c}\right), \mu = -\frac{b}{2c}, \sigma = \sqrt{-\frac{1}{2c}}$$

となります．
今日はこれまで…．また追記するかもです．

2020年10月17日追記
このアルゴリズムの利点は，反復的な処理が必要ではなく，ただ連立方程式を解くだけであるという点です．一方で問題点は，ノイズに弱いことです．このことを示すために，観測値にノイズ$\eta$が含まれていると考えます．$\eta$は平均$0$，標準偏差$\sigma_{\eta}$であるとします．

$$\hat{g}(x) = g(x) + \eta$$

ここで，$g(x)$は何らかのモデルに従っている真の値です．
すると誤差は，

$$\delta = \ln \hat{g}(x) - (a + bx + cx^2) = \ln (g(x) + \eta) - (a + bx + cx^2)$$

テイラーの定理から，第一項について近似をとると（おそらくここの五つ目の近似，さらに$\eta$は$g(x)$に比べて十分小さいことを想定していると思う），

$$\delta \approx \ln g(x) - (a + bx + cx^2) + \frac{\eta}{g(x)}$$

先程は二乗誤差の和について考えましたが，今度は$\delta^2$の平均二乗誤差を考えます．先程の例では，偏微分の計算過程でデータ数$N$を払ってしまうので和でも良かったわけですが，普通は平均二乗誤差を考えることが多いですね．

$$E[\delta^2] = \{\ln g(x) - (a + bx + cx^2)\}^2 + \frac{\sigma_{\eta}^2}{g^2(x)}$$

ぶっちゃけここが一番わかりにくいところかと思いますが，ここの9合目セクションにおいて，$f(x) = \ln g(x) - (a + bx + cx^2)$，$\epsilon = \frac{\eta}{g(x)}$，$\hat{f}(x) = 0$を代入して計算すればいいんだと思います，たぶん．
さて，上の式で大事なのは第二項で，これはもし$g(x)$が小さい，すなわちガウス関数の場合平均から離れる（=0に近い）データ点を選んでしまった場合，大きな誤差を発生させてしまうということを表しています．実際，Caruana’s algorithm ではこの問題が発生してしまいます．
そこで，Guo’s algorithm というアルゴリズムを紹介します．発想はひじょうにシンプルで，上の誤差$\delta$に$g(x)$をかけた，すなわち$g(x)$で重み付けをしたものを新たに誤差として採用すれば，小さなデータ点に対する影響も小さくなるだろうというものです．つまり，

$$\delta \approx g(x)\{\ln g(x) - (a + bx + cx^2)\} + \eta$$ $$E[\delta^2] = [g(x)\{\ln g(x) - (a + bx + cx^2)\}]^2 + \sigma_{\eta}^2$$

とすればよかろうということになります．この関数を先ほどと同様に$a, b, c$で偏微分して，不明な$g(x)$の代わりに観測値$\hat{g}(x)$を用いれば，

\begin{eqnarray}\frac{\partial S}{\partial a} &=& 2 \sum_{x} \hat{g}(x)\{ \ln \hat{g}(x) - (a + bx + cx^2) \}(-\hat{g}(x)) = 0 \\ &\Leftrightarrow& a \sum_{x} \hat{g}^2(x) + b \sum_{x} x\hat{g}^2(x) + c \sum_{x} x^2\hat{g}^2(x) = \sum_{x} \hat{g}^2(x)\ln \hat{g}(x)\end{eqnarray} \begin{eqnarray}\frac{\partial S}{\partial b} &=& 2 \sum_{x} \hat{g}(x)\{ \ln \hat{g}(x) - (a + bx + cx^2) \}(-x\hat{g}(x)) = 0 \\ &\Leftrightarrow& a \sum_{x} x\hat{g}^2(x) + b \sum_{x} x^2\hat{g}^2(x) + c \sum_{x} x^3\hat{g}^2(x) = \sum_{x} x\hat{g}^2(x) \ln \hat{g}(x)\end{eqnarray} \begin{eqnarray}\frac{\partial S}{\partial c} &=& 2 \sum_{x} \hat{g}(x)\{ \ln \hat{g}(x) - (a + bx + cx^2) \}(-x^2\hat{g}(x)) = 0 \\ &\Leftrightarrow& a \sum_{x} x^2\hat{g}^2(x) + b \sum_{x} x^3\hat{g}^2(x) + c \sum_{x} x^4\hat{g}^2(x) = \sum_{x} x^2\hat{g}^2(x) \ln \hat{g}(x)\end{eqnarray}

禍々しくなってきましたが，さっきとやってることは一緒です．よって，

\begin{equation} \begin{pmatrix} \sum_{x} \hat{g}^2(x) & \sum_{x} x\hat{g}^2(x) & \sum_{x} x^2\hat{g}^2(x) \\ \sum_{x} x\hat{g}^2(x) & \sum_{x} x^2\hat{g}^2(x) & \sum_{x} x^3\hat{g}^2(x) \\ \sum_{x} x^2\hat{g}^2(x) & \sum_{x} x^3\hat{g}^2(x) & \sum_{x} x^4\hat{g}^2(x) \end{pmatrix} \begin{pmatrix} a \\ b \\ c \end{pmatrix} = \begin{pmatrix} \sum_{x} \ln \hat{g}^2(x)\hat{g}(x) \\ \sum_{x} x \hat{g}^2(x)\ln \hat{g}(x) \\ \sum_{x} x^2 \hat{g}^2(x)\ln \hat{g}(x) \end{pmatrix} \end{equation}
となります．めでたしめでたし．また追記すると思います．

参考

R. A. Caruana, R. B. Searle, T. Heller, and S. I. Shupack, “Fast algorithm for the resolution of spectra,” Analytical Chemistry, Vol. 58, No. 6, pp. 1162―1167, 1986.

H. GUO, “A Simple Algorithm for Fitting a Gaussian Function,” IEEE Signal Process, Mag.28, No. 5 (2011), 134–137.