はじめに

一般に、最小二乗法と最尤推定法は同じようなものです。最尤推定法の特別な場合（誤差が正規分布）が最小二乗法です。基本的な内容ですが、あまり途中計算過程も含めてこれらをまとめている文献も少ないので、自分用にまとめました。

最小二乗法

まず、線形の連立方程式は以下に分類されます。$m$を式の数、$n$を変数の数とします。

方程式の数が未知変数より多い：$m > n$　⇒解が一意に定まらない。
方程式の数＝未知変数の数：$m=n$　⇒解析的に解ける
方程式の数が未知変数より少ない：$m < n$ ⇒解が求まらない
最小二乗法は1)の場合にもっともらしい解を求めます。なお、3)の近似解法はラグランジュの未定乗数法と呼ばれています。

$x$を未知変数として、連立方程式を書き下すと、

\left\{
\begin{array}{c}
a_{11}x_{1}+a_{12}x_{2}+\cdots+a_{1n}x_{n} &= y_{1}\\
a_{21}x_{1}+a_{22}x_{2}+\cdots+a_{2n}x_{n} &= y_{2}\\
\cdots\\
a_{m1}x_{1}+a_{m2}x_{2}+\cdots+a_{mn}x_{n} &= y_{m}\\
\end{array}
\right.

これを行列で書くと、

Ax = y

ここで、下記の2乗誤差関数を最小化する$x$が求める解となります。また、$y$は観測値に対応します。

J = \frac{1}{2}  \sum_{i=1}^{m}(\sum_{j=1}^{n} a_{ij}x_{j}-y_i)^2 \rightarrow min

2次関数なので、最小となる解を求めるには微分して0のところを求めればよいことになります。つまり、

\frac{\partial}{\partial x}(Ax-y)^{T} (Ax-y)=0 \\

微分しやすいように変形して、

\frac{\partial}{\partial x} (x^{T}A^{T}-y^{T})(Ax-y)=0 \\
\frac{\partial}{\partial x} (x^{T}A^{T}Ax-x^{T}A^{T}y-y^{T}Ax+y^{T}y)=0 \\

スカラーのベクトル微分に注意して、

2A^{T}Ax-A^{T}y-A^{T}y = 0 \\
2(A^{T}Ax-A^{T}y) = 0

したがって、

A^{T}Ax-A^{T}y = 0 \\
x = (A^{T}A)^{-1}A^{T}y

これはいわゆる正規方程式$A^T Ax=A^T b$の解となります。また、$A^{\dagger}=(A^{T}A)^{-1}A^{T}$を一般逆行列、またはMoore-Penroseの疑似逆行列と呼びます。そもそも、$m=n$であれば$A$は正則行列となり、逆行列が定義でき、$x=A^{-1}y$と解くことができますが、最小二乗法の場合は一般に$m>n$であるため、逆行列が定義できないが、その場合は一般逆行列を用いて

x = A^{\dagger}y

のように解くことが可能です。

最尤推定法(Maximum Likelihood Estimation)

最尤推定法とは、ベイズ統計に基づくベイズの定理を用いて、誤差のある観測値から状態量の推定値を求める理論です。まず、ベイズの定理とは、以下のように表されます。

p(\theta|y)=p(\theta)\frac{p(y|\theta)}{y}

ここで、$p(A|B)$は条件付き確率を表し、$B$が起こった時の$A$の起こる確率です。$\theta$は未知のパラメータ、$y$は観測値を表します。つまり、この式は、観測値$y$を観測したときの、$\theta$の分布（事後分布）を、与えられた$\theta$の事前分布$p(\theta)$と尤度$p(y|\theta)$をもとに更新するという意味になります。この考え方が、現代制御理論のキモであるカルマンフィルタに繋がります。

さて、状態量$x$の事後確率分布は、正規分布を仮定すると、以下のようになります。

p(x|y)=\frac{1}{\sqrt{2\pi\sigma^2}}exp\biggl(-\frac{(x-\hat{x})^2}{2\sigma^2}\biggr)

ここで、平均値$\hat{x}$、分散$\sigma^2$の正規分布となります。また、$x=\hat{x}$のとき事後確率分布関数の値は最大となります。事後確率分布を最大にする$x=\hat{x}$を推定値とする方法を最尤推定法といいますが、今回は正規分布を仮定していたので平均値が最尤推定値となることは自明でした。

さて、次に最小二乗法と一致することを示します。ある状態量$x$における観測値$y_{n} (n=1,2,\cdots,N)$のデータ点が$N$個、正規分布の形をした乱数で用意されているとします。すると、確率密度関数$P$は以下のように書けます。

P=\prod_{n}^{N}\frac{1}{\sqrt{2\pi\sigma^2}}exp\biggl(-\frac{(y_{n}-f(x))^2}{2\sigma^2}\biggr)\\
=\biggl(\frac{1}{2\pi\sigma^2}\biggr)^{N/2}exp\biggl(-\frac{1}{2\sigma^2}\sum_{n}^{N}(y_{n}-f(x))^2\biggr)

ここで、右辺のexpの中に、2乗誤差が含まれています。最小二乗法の章と同様、これを$J$とおくと、

P=\biggl(\frac{1}{2\pi\sigma^2}\biggr)^{N/2}exp\biggl(-\frac{1}{\sigma^2}J\biggr)

計算を簡単にするために、$\beta=1/\sigma^2$とおきます。また、2乗誤差$J$は未知状態量$x$に依存しています。したがって、パラメータの依存性を明確にすると、

P(\beta,x)=\biggl(\frac{\beta}{2\pi}\biggr)^{N/2}e^{-\beta J(x)}

これを最大化するには対数を取り、対数尤度関数にして最大値を求めます。

\ln P(\beta,x)=\frac{N}{2}\ln\beta - \frac{N}{2}\ln 2\pi - \beta J(x)

$\partial\ln P/\partial x = 0$の解は、すなわち$\partial J(x)/\partial x = 0$を解くことにほかならず、これは最小二乗法そのものです。したがって、正規分布の時の最尤推定法は最小二乗法に一致します。

なお、$\partial\ln P/\partial \beta = 0$の解は、

\frac{1}{\beta}=\frac{2J}{N}

となり、標準偏差は$\sigma=\sqrt{1/\beta}$は、$\sigma=\sqrt{2J/N}=E_{RMS}$となって、平均2乗誤差の平方根(RMS)に一致します。