要旨
日本語で、回帰式が二次以上の最小二乗法の導出について、丁寧に記述したものがなかったので書きます。
また、それぞれのデータ点に重みを付けることで、データ間にバイアスをつけ回帰曲線を求めることができる最小二乗法を拡張したアルゴリズムの紹介も行います。
導出
以下のような$n$個のデータセットがあるとします
\{(x_{1}, y_{1}), (x_{2}, y_{2}), ... , (x_{n}, y_{n})\}\label{a}\tag{1}
これを以下で示す$m$次の多項式$y$で近似することを考えます。
y(x)=\omega_{0}x^0+\omega_{1}x^1+\omega_{2}x^2+...+\omega_{m}x^m\label{b}\tag{2}
ここで、$\omega_{i}$は多項式の係数です。
最小二乗法は、$(\ref{a})$で示す$n$個のデータセットと、$(\ref{b})$で示す多項式$y$との二乗誤差和が最小となる$\omega$の組み合わせを導出するアルゴリズムです。
まずは簡単に、$(x_{1}, y_{1})$と$y$との誤差を考えてみます。誤差(error)は、$x_{1}$を$y$に入れたときに出てくる値と、実際の値$y_{1}$との差の絶対値なので、
\begin{align}
error&=|y_{1}-y(x_1)|\\
&=|y_{1}-(\omega_{0}x^0_1+\omega_{1}x^1_1+\omega_{2}x^2_1+...+\omega_{m}x^m_1)|
\end{align}
したがって、二乗誤差(squared error)は、
\begin{align}
squared\,error&=|y_{1}-(\omega_{0}x^0_1+\omega_{1}x^1_1+\omega_{2}x^2_1+...+\omega_{m}x^m_1)|^2
\end{align}
となります。これと同じことを$(x_{2}, y_{2}), ... , (x_{n}, y_{n})$でも行うと、
\begin{align}
|y_{1}-(\omega_{0}x^0_1+\omega_{1}x^1_1+&\omega_{2}x^2_1+...+\omega_{m}x^m_1)|^2\\
|y_{2}-(\omega_{0}x^0_2+\omega_{1}x^1_2+&\omega_{2}x^2_2+...+\omega_{m}x^m_2)|^2\\
&\vdots\\
|y_{1}-(\omega_{0}x^0_n+\omega_{1}x^1_n+&\omega_{2}x^2_n+...+\omega_{m}x^m_n)|^2\label{f}\tag{3}\\
\end{align}
二乗誤差和を$E$とすると、$E$は$(\ref{f})$の和となります。$E$は行列で以下のように書くことができます。
E=
\left|
\left[
\begin{array}{c}
y_{1} \\
y_{2} \\
\vdots \\
y_{n} \\
\end{array}
\right]
-
\left[
\begin{array}{cccc}
{x^0_{1}} & {x^1_{1}} & \cdots & {x^m_{1}} \\
{x^0_{2}} & {x^1_{2}} & \cdots & {x^m_{2}} \\
\vdots & \vdots & \ddots & \vdots \\
{x^0_{n}} & {x^1_{n}} & \cdots & {x^m_{n}} \\
\end{array}
\right]
\left[
\begin{array}{c}
\omega_{0} \\
\omega_{1} \\
\vdots \\
\omega_{m} \\
\end{array}
\right]
\right|^2\label{g}\tag{4}
$(\ref{g})$を簡単にするために、以下のように行列を表すことにします。
\begin{align}
\boldsymbol{Y}&=
\left[
\begin{array}{c}
y_{1} \\
y_{2} \\
\vdots \\
y_{n} \\
\end{array}
\right]
\in\mathbb{R}^{n},\\
\boldsymbol{X}&=
\left[
\begin{array}{cccc}
{x^0_{1}} & {x^1_{1}} & \cdots & {x^m_{1}} \\
{x^0_{2}} & {x^1_{2}} & \cdots & {x^m_{2}} \\
\vdots & \vdots & \ddots & \vdots \\
{x^0_{n}} & {x^1_{n}} & \cdots & {x^m_{n}} \\
\end{array}
\right]\in\mathbb{R}^{n\times (m+1)},\\
\boldsymbol{W}&=
\left[
\begin{array}{c}
\omega_{0} \\
\omega_{1} \\
\vdots \\
\omega_{m} \\
\end{array}
\right]
\in\mathbb{R}^{m+1}.
\end{align}
すると、$(\ref{g})$は、
E=|\boldsymbol{Y-XW}|^2\label{h}\tag{5}
と簡潔に表現できます。$(\ref{h})$の$E$が最小となる$\boldsymbol{W}$を求めればいいわけです。
$E$が最小となるとき、$\boldsymbol{W}$での偏微分は0になります。したがって、
\frac{\partial E}{\partial\boldsymbol{W}}=0\label{i}\tag{6}
を満たす$\boldsymbol{W}$を求めればいいことになります（そこが極小値となることを示す厳密な証明はカット）。
$(\ref{h})$、$(\ref{i})$より、
\begin{align}
\frac{\partial}{\partial\boldsymbol{W}}|\boldsymbol{Y-XW}|^2=0\label{ii}\tag{7}
\end{align}
を解いていき、$\boldsymbol{W}$を求めます。
$|\boldsymbol{Y-XW}|^2=(\boldsymbol{Y-XW})^T(\boldsymbol{Y-XW})$なので（ここで$T$は転置行列です）、$(\ref{ii})$は以下のように変換できます
\begin{align}
\frac{\partial}{\partial\boldsymbol{W}}(\boldsymbol{Y-XW})^T&(\boldsymbol{Y-XW})=0\label{j}\tag{8}\\
\end{align}
ここで、ベクトルの微分の公式として
\begin{align}
\frac{\partial}{\partial\boldsymbol{x}}\boldsymbol{x^TAx}=\boldsymbol{(A+A^T)x}\label{jj}\tag{9}\\
\end{align}
というのがあります（証明はこのサイト）。
$\boldsymbol{A=I}$（Iは単位行列）のとき、$(\ref{jj})$は、
\begin{align}
\frac{\partial}{\partial\boldsymbol{x}}\boldsymbol{x^Tx}=2\boldsymbol{x}\label{jjj}\tag{10}\\
\end{align}
が成り立ちます。
また、ベクトルの微分においては、以下の式が成り立ちます（証明は同様にこのサイト）。
\begin{align}
\frac{\partial}{\partial\boldsymbol{x}}\boldsymbol{Ax}=\boldsymbol{A^T}\label{jjjj}\tag{11}\\
\end{align}
$(\ref{jjj})$と$(\ref{jjjj})$を使うと、$(\ref{j})$は以下のように書くことができます、
\begin{align}
\frac{\partial}{\partial\boldsymbol{W}}(\boldsymbol{Y-XW})^T&(\boldsymbol{Y-XW})=0\\
\frac{\partial(Y-XW)}{\partial\boldsymbol{W}}\frac{\partial}{\partial\boldsymbol{(Y-XW)}}(\boldsymbol{Y-XW})^T&(\boldsymbol{Y-XW})=0\\
-X^T-2(\boldsymbol{Y-XW})&=0\\
\boldsymbol{X}^T(\boldsymbol{Y-XW})&=0\label{k}\tag{12}
\end{align}
$(\ref{k})$を展開して整理すると、
\begin{align}
\boldsymbol{X^TY}-\boldsymbol{X^TXW}&=0\\
\boldsymbol{X^TXW}&=\boldsymbol{X}^T\boldsymbol{Y}\\
\boldsymbol{(X^TX)}^{-1}\boldsymbol{X^TXW}&=\boldsymbol{(X^TX)}^{-1}\boldsymbol{X}^T\boldsymbol{Y}\\
\boldsymbol{W}&=\boldsymbol{(X^TX)}^{-1}\boldsymbol{X}^T\boldsymbol{Y}\label{l}\tag{13}
\end{align}
となります。$(\ref{l})$の$\boldsymbol{W}$が最小二乗法で求めるべき多項式の係数となります。
Pythonでの実装
ついでに、ノイズがのった三次関数を作成して、先ほど求めた最小二乗法をPythonにて実装し、係数を推定してみます。コードは以下のような感じになります。
import numpy as np
import matplotlib.pyplot as plt
import os
fig = plt.figure(figsize=(19.2*0.8,10.03),dpi=300)
ax1 = fig.add_subplot(1,1,1)
# 最小二乗法の記述
def LeastSquareMethod(x, y, degree):
x = x.reshape(-1, 1)
y = y.reshape(-1, 1)
X = np.empty((len(x),0))
for i in range(0, degree+1):
current_x = x ** i
X = np.hstack((X, current_x))
w = np.dot(np.dot(np.linalg.inv(np.dot(X.T, X)), X.T), y)
w = w.reshape(-1)
return w
# 多項式の係数と入力を受け取ってグラフを作成する関数
def f(x,w):
y = 0
for i, w in enumerate(w):
y += w * (x ** i)
return y
# サンプルデータ生成
x = np.linspace(-5, 5, 100)
# 例として次数3の多項式(y = -0.2x*x*x - 0.4x*x + 3x + 1)を用意
true_coefficients = np.array([1, 3, -0.4, -0.2])
y_true = f(x, true_coefficients)
# ノイズを付加
noise = np.random.normal(0, 1, len(x))
y_noisy = y_true + noise
# 最小二乗法を適用
w = LeastSquareMethod(x, y_noisy, 3)
# グラフ作成
ax1.scatter(x, y_noisy, s=100, label='Noisy Data')
ax1.plot(x, f(x,w), color='maroon', linewidth=5.0, label='regression curve')
ax1.set_xlabel('$x$', fontsize=40)
ax1.set_ylabel('$y$', fontsize=40)
ax1.tick_params(labelsize=30)
ax1.grid(linestyle='--')
ax1.legend(loc='best', fontsize=25, ncol=1, edgecolor='black', fancybox=False)
# 最小二乗法によって導出された係数
print(w)
# [ 0.85604308 2.96769492 -0.39790196 -0.19977282] (x^0, x^1, x^2, x^3 の順)
plt.tight_layout()
plt.savefig('least_square_method.png')
生成されたノイズと最小二乗法により推定されたカーブは以下のようになります。
データ点に重みを付けた最小二乗法
先ほどまでは、全てのデータ点を"等しく"扱い二乗誤差和が最小となる多項式の係数を求めていました。ここでは、データ点に"バイアスを付け"、重み付き二乗誤差和から多項式の係数を求める手法を軽くですが紹介してみようと思います。簡単にいうと、どのデータ点を回帰に重要視するかを決定することができるアルゴリズムとなっています。
以下のようなデータ点を考えてみます。
x = np.array([1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 5, 5, 6, 6, 6, 7, 7, 7, 8, 8, 8, 9, 9, 9])
y = np.array([1, 2, 3, 1, 2, 3, 1, 2, 3, 7, 8, 9, 7, 8, 9, 7, 8, 9, 4, 5, 6, 4, 5, 6, 4, 5, 6])
beta = np.array([10, 1, 1, 1, 10, 1, 1, 1, 10, 10, 1, 1, 10, 1, 1, 10, 1, 1, 10, 1, 1, 1, 10, 1, 1, 1, 10])
ここで、それぞれのデータ点にbetaとして重みが与えてられています。図示すると以下のようになります。betaの値が大きいほどそのデータの評価が高いことを示します。
これに、performence-weighted least square method, PWLSと呼ばれる最小二乗法を拡張した、データ点に"バイアス"を与え回帰曲線を求めることができるアルゴリズムを適用すると、以下のように回帰曲線が求まります。
Least Square Methodは最小二乗法で求まった回帰曲線です。最小二乗法と比較するとPWLSは、betaの大きい、つまり評価値の高いデータ点の近くを通り回帰曲線が引かれていることがわかります。
コードは以下の通りです。
def PWLS(x, y, beta, degree):
"""
This function calculates the coefficients of polynomial that minimizes the weighted square error.
x, y (NumPy array): These represent the data points. Please provide 1-dimensional NumPy arrays.
e.g., x = np.array([0.0, 1.0, 2.5]), y = np.array([1.0, 0.5, 3.5])
beta (NumPy array): Weights corresponding to each data point. Please provide a 1-dimensional NumPy array of the same length as x and y.
e.g., beta = np.array([1.0, 1.0, 2.0])
degree (int): This denotes the degree of the polynomial function. For example, setting degree=2 will calculate a quadratic curve.
"""
x = x.reshape(-1, 1)
Y = y.reshape(-1, 1)
beta = beta.reshape(-1, 1)
B = np.hstack([beta]*(degree+1))
X = np.empty((len(x),0))
for i in range(0, degree+1):
current_x = x ** i
X = np.hstack((X, current_x))
W = np.dot(np.dot(np.linalg.inv(np.dot((B*X).T, (B*X))), (B*X).T), beta*Y)
W = W.reshape(-1)
return W
上記で示したdef LeastSquareMethod(x, y, degree)と使い方はほぼ同じで、入力に重みのbetaが増えた形になっています。
以下の論文が原著となっています。
S Koseki, M Hayashib, and D Owaki. "Identifying essential factors for energy-efficient walking control across a wide range of velocities in reflex-based musculoskeletal systems." PLOS Computational Biology 20.1 (2024): e1011771.
https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1011771
ちなみに、このアルゴリズムの提案者は筆者です。
詳しくは論文内にも書いてあるのですが、後日時間のあるときに、このアルゴリズムの導出について日本語で説明する記事を作成しようと思っています。