複数の点(データ)が与えられたときに、それらによく当てはまる関数を求めること(特に、パラメータを持つ関数について最適なパラメータを決定すること)を、統計学の用語で回帰(regression)と言います。与えられた点を通る必要はないため、前回と前々回の授業で扱った補間とは異なります。
$xy$ 平面上の $x$ 座標が相異なる $n$ 個の点 $(x_1,y_1),\ldots,(x_n,y_n)$ に対して、$m$ 個の関数 $g_1(x),\ldots,g_m(x)$ の線形結合として表される
$$ g(x)=\sum_{j=1}^m c_j g_j(x) $$の最適なパラメータ $c_1,\ldots,c_m\in\mathbb{R}$ を決定することを考えます。例えば
$$ g_j(x)=x^{j-1}\quad (j=1\ldots,m) $$の場合は $g(x)$ は $m-1$ 次の多項式、特に $m=2$ の場合は $g(x)$ は $1$ 次関数(直線)になり、実用上このように指定することが多いです。
(※最大次数の項の係数が $0$ になる可能性があるため、上記は厳密には「$m-1$ 次以下」です。)
あとは何をもって「パラメータが最適である」と言うかが問題となりますが、最小二乗法(least squares method)においては、残差平方和(residual sum of squares)
$$ \sum_{i=1}^n \left(g(x_i)-y_i\right)^2=\sum_{i=1}^n \left(\sum_{j=1}^m c_j g_j(x_i)-y_i\right)^2 $$が最小となるようにパラメータ $c_1,\ldots,c_m\in\mathbb{R}$ を決定します。これは、大雑把に言えば「関数値とデータのずれを最小にしよう」という考え方です。
特に、ある関数 $f$ 上の点 $(x_1,f(x_1)),\ldots,(x_n,f(x_n))$ が与えられた場合には
$$ \sum_{i=1}^n \left(g(x_i)-f(x_i)\right)^2=\sum_{i=1}^n \left(\sum_{j=1}^m c_j g_j(x_i)-f(x_i)\right)^2 $$が最小となるようにパラメータ $c_1,\ldots,c_m\in\mathbb{R}$ を決定することになり、このようにして得られる関数 $g(x)$ は $f(x)$ の近似と解釈することができます。
以下では、最小二乗法における最適なパラメータ $c_1,\ldots,c_m\in\mathbb{R}$ を理論的に求めます。
$n\times m$ 行列 $A$、$m$ 次元ベクトル $c$、$n$ 次元ベクトル $b$ を
$$ A= \begin{pmatrix} g_1(x_1) & g_2(x_1) & \cdots & g_m(x_1)\\ g_1(x_2) & g_2(x_2) & \cdots & g_m(x_2)\\ \vdots & \vdots & \vdots & \vdots\\ g_1(x_n) & g_2(x_n) & \cdots & g_m(x_n) \end{pmatrix} ,\quad c= \begin{pmatrix} c_1\\ c_2\\ \vdots\\ c_m \end{pmatrix} ,\quad b= \begin{pmatrix} y_1\\ y_2\\ \vdots\\ y_n \end{pmatrix} $$とおくと、残差平方和について
$$ \sum_{i=1}^n \left(g(x_i)-y_i\right)^2=\sum_{i=1}^n \left(\sum_{j=1}^m c_j g_j(x_i)-y_i\right)^2=(Ac-b)^T(Ac-b) $$が成り立ちます。そのため、関数 $r:\mathbb{R}^m\to\mathbb{R}$ を
$$ r(c)=(Ac-b)^T(Ac-b) $$により定めれば、最小二乗法では $r(c)$ を最小にするようなベクトル $c$ を求めればよいことが分かります。
ここで、$r(c)$ の勾配ベクトルを計算すると
$$ \nabla r(c)= \begin{pmatrix} \frac{\partial r(c)}{\partial c_1}\\ \vdots\\ \frac{\partial r(c)}{\partial c_m} \end{pmatrix} =2A^T(Ac-b)=2A^TAc-2A^Tb $$となります(証明は省略、$1$ 変数関数 $(ax-b)^2$ の微分が $2a(ax-b)$ であることを考えると理解しやすいです)。
$r(c)$ は凸関数であるため、$0_m$ を $m$ 次元の零ベクトルとして
$$ \nabla r(c)=0_m $$を満たす $c$ を求めれば、それが $r(c)$ を最小にするベクトルになります。
この条件を変形すると連立一次方程式
$$ A^TAc=A^Tb $$が得られ、これを正規方程式(normal equation)と呼びます。
正規方程式が唯一の解を持つかどうかは、行列 $A^TA$ の正則性に依存します。$A^TA$ を $A$ のグラム行列(Gram matrix)と呼び、次の二つは同値であることが知られています。
$A^TA$ が正則
$A$ の列ベクトルが線形独立
なお、$n\times m$ 行列 $A$ は $m$ 個の $n$ 次元ベクトルを列ベクトルとして持つため、もし $n<m$ であれば明らかに $A$ の列ベクトルは線形独立になりません。よって、最小二乗法を用いる上では、前提条件として少なくとも $n\geq m$ が必要であることに注意してください。
$n$ 個の点 $(x_1,y_1),\ldots,(x_n,y_n)$ に対して、最適な $1$ 次関数
$$ g(x)=c_1+c_2x $$を求めることを考えます。これを線形回帰(linear regression)と言い、求めた $1$ 次関数のグラフを回帰直線(regression line)と呼びます。
最小二乗法を用いるならば、前述の式で
$$ m=2,\quad g_1(x)=1,\quad g_2(x)=x $$の場合を考えることになるため、
$$ A= \begin{pmatrix} 1 & x_1\\ 1 & x_2\\ \vdots & \vdots\\ 1 & x_n \end{pmatrix} ,\quad c= \begin{pmatrix} c_1\\ c_2 \end{pmatrix} ,\quad b= \begin{pmatrix} y_1\\ y_2\\ \vdots\\ y_n \end{pmatrix} $$とおいて正規方程式
$$ A^TAc=A^Tb $$を解けばよいことが分かります。
次の $11$ 個の点に対して、正規方程式を解いて最小二乗法による線形回帰を行い、それらの点および求めた回帰直線(範囲は $0\leq x\leq 1$)を描画してください。
$$ (0,1.09),\quad (0.1,1.17),\quad (0.2,1.42),\quad (0.3,1.54), $$$$ (0.4,1.71),\quad (0.5,1.99),\quad (0.6,2.26),\quad (0.7,2.44), $$$$ (0.8,2.55),\quad (0.9,2.84),\quad (1,3.04) $$%演習1のコード
data = [0,1.09; 0.1,1.17; 0.2,1.42; 0.3,1.54; 0.4,1.71; 0.5,1.99; 0.6,2.26; 0.7,2.44; 0.8,2.55; 0.9,2.84; 1,3.04]; %点を並べた行列
n = size(data,1); %点の個数
%コードの最初の行に「%plot --format svg」と書くと、グラフがきれいに表示される
%点を描画するには、関数plotの3つ目の引数に「"."」などを指定すればよい
%行列の転置は「'」を使う
%連立一次方程式を解く際は「\」を使う
最小二乗法による線形回帰では、正規方程式を解かなくても $(x_1,y_1),\ldots,(x_n,y_n)$ に関する統計量を用いて
$$ c_2=\frac{s_{xy}}{s_x^2}, $$$$ c_1=\overline{y}-c_2\overline{x} $$と計算することができます。
ここで、$\overline{x}$ は $x_1,\ldots,x_n$ の平均、$\overline{y}$ は $y_1,\ldots,y_n$ の平均、$s_x^2$ は $x_1,\ldots,x_n$ の分散、$s_{xy}$ は $x_1,\ldots,x_n$ と $y_1,\ldots,y_n$ の共分散です。
演習1と同じ $11$ 個の点に対して、統計量を用いた計算で最小二乗法による線形回帰を行い、それらの点および求めた回帰直線(範囲は $0\leq x\leq 1$)を描画してください。
ただし、統計量の計算には平均はmean(ベクトル)
、分散はvar(ベクトル)
、共分散はcov(ベクトル1,ベクトル2)
を使用してよいです。
%演習2のコード
関数 $f(x)$ と $n$ 個の節点 $x_1,\ldots,x_n$ が与えられたとき、最小二乗法によって $f(x)$ を近似する $m-1$ 次以下の多項式
$$ g(x)=\sum_{j=1}^m c_j x^{j-1} $$を求めることができます(ただし、$n\geq m$)。前述のように、これは
$$ g_j(x)=x^{j-1}\quad (j=1\ldots,m) $$の場合に相当します。
節点
$$ x_i=-1+\frac{2(i-1)}{n-1}\quad (i=1,\ldots,n) $$を用いて関数
$$ f(x)=\frac{1}{1+25x^2} $$を $m-1$ 次以下の多項式
$$ g(x)=\sum_{j=1}^m c_j x^{j-1} $$で近似することを考えます。
$n=30,60$、$m=10,20$ を組み合わせた計 $4$ 個の場合について近似多項式を求め、元の関数 $f(x)$ とともに $-1\leq x\leq 1$ の範囲でグラフを描画してください。
%演習3のコード
演習1~演習3に取り組んでください。