Maximaで学ぶ最小二乗法

数学の具体的な計算にMaximaを使って、数学もMaximaも同時に学んでしまいましょう。

今回はMaximaを使って統計学における最小二乗法をみてみたいと思います。

最小二乗法とは、複数のデータの組を適当な関数でフィッティングする方法のことです。

本記事では簡単のために１次関数によるフィッティングを例にとりたいと思います。

平均・分散・標準偏差
共分散・相関係数
フィッティング
最小二乗法（統計の言葉で）
最小二乗法（行列の言葉で）

平均・分散・標準偏差

２つのデータからなる $N$ 組のデータ $(x_i,y_i)$ （ $i=1,2,\dots,N$ ）があるとき、

それぞれの成分の平均値は

$\displaystyle{\ {\rm Avr}_x=\frac{\sum_{i=1}^N x_i}{N}, \qquad {\rm Avr}_y=\frac{\sum_{i=1}^N y_i}{N}}$

で定義されます。例えば10組のデータ

$\displaystyle{(1,3),~(2,2),~(2,3),~(3,4),~(4,3),~(4,5),~(6,4),~(6,6),~(7,6),~(8,8)}$

に対しては：

kill(all)$
N:10$
xy:[[1,3],[2,2],[2,3],[3,4],[4,3],[4,5],[6,4],[6,6],[7,6],[8,8]]$
Avr_x:sum(xy[i][1],i,1,N)/N;
Avr_y:sum(xy[i][2],i,1,N)/N;

より

$\displaystyle{\ {\rm Avr}_x=\frac{43}{10}, \qquad {\rm Avr}_y=\frac{22}{5}}$

と求まります。ファイルから読み込んだデータに本記事の解析を適用しているのが

以下の記事になります。本記事の続編として、あわせて読んでみてください：

pianofisica.hatenablog.com

分散は

$\displaystyle{\ {\rm Var}_x=\frac{\sum_{i=1}^N (x_i-{\rm Avr}_x)^2}{N}=\frac{\sum_{i=1}^N x_i^2}{N}-{\rm Avr}_x^2 \\ \ {\rm Var}_y=\frac{\sum_{i=1}^N (y_i-{\rm Avr}_y)^2}{N}=\frac{\sum_{i=1}^N y_i^2}{N}-{\rm Avr}_y^2}$

で定義されます。データ例に対しては

Var_x:sum((xy[i][1]-Avr_x)^2,i,1,N)/N;
Var_y:sum((xy[i][2]-Avr_y)^2,i,1,N)/N;

より

$\displaystyle{\ {\rm Var}_x=\frac{501}{100}, \qquad {\rm Var}_y=\frac{76}{25}}$

と計算されます。

標準偏差は分散の二乗根で定義されます：

$\displaystyle{\ \sigma_x=\sqrt{{\rm Var}_x}, \qquad \sigma_y=\sqrt{{\rm Var}_y}}$

データ例に対しては

sigma_x:sqrt(Var_x);
sigma_y:sqrt(Var_y);

より

$\displaystyle{\ {\rm Var}_x=\frac{\sqrt{501}}{10}, \qquad {\rm Var}_y=\frac{2\sqrt{19}}{5}}$

と求まります。

共分散・相関係数

共分散は

$\displaystyle{\ {\rm Cov}(X,Y)=\frac{\sum_{i=1}^N(x_i-{\rm Avr}_x)(y_i-{\rm Avr}_y)}{N}=\frac{\sum_{i=1}^N x_i y_i}{N}-{\rm Avr}_x{\rm Avr}_y}$

で定義されます。データ例に対しては

Cov:sum((xy[i][1]-Avr_x)*(xy[i][2]-Avr_y),i,1,N)/N;

より

$\displaystyle{\ {\rm Cov}(X,Y)=\frac{169}{50}}$

と計算されます。

相関係数は

$\displaystyle{\ \rho=\frac{{\rm Cov}(X,Y)}{\sigma_x\sigma_y}}$

で定義されます。データ例に対しては

rho:Cov/(sigma_x*sigma_y);
bfloat(rho);

より

$\displaystyle{\ \rho\simeq0.866}$

と求まります。相関係数は１に近いほど確率変数は線形な関係になります。

フィッティング

フィッティングして得られる１次関数は

$\displaystyle{\ y=Ax+B; \qquad A=\frac{{\rm Cov(x,y)}}{\sigma_x^2}, \qquad B={\rm Avr}_y-\frac{{\rm Cov(x,y)}{\rm Avr}_x}{\sigma_x^2}}$

によって与えられます。なぜこのような係数になるかについては下で説明しています。

データ例に対しては

A:Cov/sigma_x^2;
B:Avr_y-Cov*Avr_x/sigma_x^2;

より次の１次関数を得ます： $\displaystyle{y=\frac{338}{501}x+\frac{751}{501}}$

データ点と１次関数を

A:Cov/sigma_x^2;
B:Avr_y-Cov*Avr_x/sigma_x^2;
plot2d([ [discrete, xy], A*x+B], [x, 0,10],
[style, points, lines],[color, red, blue],
[point_type, asterisk],[legend, "data", "fitting"])$

によってプロットしたのが以下の図です：

f:id:pianofisica:20190519041941p:plain

最小二乗法（統計の言葉で）

さて、上で求めた１次関数は、１次関数とデータ点とのずれ

$\displaystyle{\ \Delta y_i = y-y_i= Ax_i+B-y_i}$

の二乗和

$\displaystyle{\ \Delta(A,B)=\sum_{i=1}^N\left(\Delta y_i\right)^2}$

を最小化することによって求められたものです。

すなわち、最小化するための必要条件

$\displaystyle{\ \frac{\partial\Delta(A,B)}{\partial A}=0, \qquad \frac{\partial\Delta(A,B)}{\partial B}=0}$

を $A, \ B$ について解いて求められたものが上で与えた１次関数です。

実際、これらの条件から得られる $A, \ B$ に対する２元連立１次方程式

$\displaystyle{\ 0={\rm Cov}(X,Y)+{\rm Avr}_x{\rm Avr}_y-A\left({\rm Var}_x+{\rm Avr}_x^2\right)-B{\rm Avr}_x,\\ \ 0={\rm Avr}_y-A{\rm Avr}_x-B}$

を例えば

solve([0=c+ax*ay-a*(vx+ax^2)-b*ax, 0=ay-a*ax-b],[a,b]);

によって解くことで

$\displaystyle{A=\frac{{\rm Cov(x,y)}}{{\rm Var}_x}, \qquad B={\rm Avr}_y-\frac{{\rm Cov(x,y)}{\rm Avr}_x}{{\rm Var}_x}}$

が得られます。これがフィッティングの式に使われていた係数です。

最小二乗法（行列の言葉で）

さて、ずれの二乗和を行列で表してみましょう：

$\displaystyle{\ \left(\begin{array}{c} \Delta y_1\\ \Delta y_2\\ \vdots \\ \Delta y_N \end{array}\right)=\left(\begin{array}{c} Ax_1+B-y_1\\ Ax_2+B-y_2\\ \vdots \\ A_N+B-y_N \end{array}\right)=\left(\begin{array}{cc} x_1 & 1 \\ x_2 & 1 \\ \vdots \\ x_N & 1 \end{array}\right)\left(\begin{array}{c} A \\ B \end{array}\right)-\left(\begin{array}{c} y_1\\ y_2\\ \vdots \\ y_N \end{array}\right)}$

ここで行列

$\displaystyle{\ M=\left(\begin{array}{cc} x_1 & 1 \\ x_2 & 1 \\ \vdots \\ x_N & 1 \end{array}\right), \qquad X=\left(\begin{array}{c} A \\ B \end{array}\right), \qquad Y=\left(\begin{array}{c} y_1\\ y_2\\ \vdots \\ y_N \end{array}\right)}$

を導入すると、ずれの二乗和が

$\displaystyle{\ \Delta(A,B)=\sum_{i=1}^N\left(\Delta y_i\right)^2=(MX-Y)^{\rm T}(MX-Y)}$

と表せます。ここで一般に行列 $K$ に対して $K^{\rm T}$ はその転置行列を表すものとします。

ずれの二乗和 $\Delta(A,B)$ を最小化する必要条件として、

$\ X$ の各成分についての偏微分係数がゼロになることを課すことにより

$\displaystyle{\ 0=M^{\rm T}MX-M^{\rm T}Y}$

を得ます。これは以下のようにして示されます。まず

$\displaystyle{\ \Delta(A,B)=(X^{\rm T}M^{\rm T}-Y^{\rm T})(MX-Y)\\ \quad\quad\quad\ \ \ =X^{\rm T}M^{\rm T}MX-X^{\rm T}M^{\rm T}Y-Y^{\rm T}MX+Y^{\rm T}Y\\ \quad\quad\quad\ \ \ = XM^{\rm T}MX-2X^{\rm T}M^{\rm T}Y+Y^{\rm T}Y}$

です。ここで

$\displaystyle{\ -X^{\rm T}M^{\rm T}Y=-X_i(M^{\rm T})_{ij}Y_j=-X_iM_{ji}Y_j=-Y_jM_{ji}X_i=-Y^{\rm T}MX}$

を使いました。ただし添字についてはアインシュタインの規約に従うものとします。

$\Delta(A,B)$ を $X_k$ について微分すると

$\displaystyle{\ \frac{\partial\Delta(A,B)}{\partial X_k} =\frac{\partial}{\partial X_k} ( X_i(M^{\rm T}M)_{ij}X_j -2X_i (M^{\rm T} )_{ij} Y_j+ Y_iY_i )\\ \qquad\qquad \ \ \ = (M^{\rm T}M)_{kj}X_j+X_i(M^{\rm T}M)_{ik}-2 (M^{\rm T} )_{kj} Y_j}$

です。さらに

$\displaystyle{\ (M^{\rm T}M)_{ik}=(M^{\rm T})_{il}M_{lk}=M_{li}M_{lk}=(M^{\rm T})_{kl}M_{li}=(M^{\rm T}M)_{ki}}$

より

$\displaystyle{\ \frac{\partial\Delta(A,B)}{\partial X_k} = 2(M^{\rm T}M)_{kj}X_j-2 (M^{\rm T} )_{kj} Y_j}$

したがって条件式

$\displaystyle{\ 0=M^{\rm T}MX-M^{\rm T}Y}$

を得ます。

もし行列 $M^{\rm T}M$ が正則であればこの条件式を $X$ について解くことができます：

$\displaystyle{\ X=\left(M^{\rm T}M\right)^{-1}M^{\rm T}Y}$

これによって $A, ~ B$ の値を求めることができるというわけです。

データ例からこれらの行列を計算してみると

M:matrix([xy[1][1],1],[xy[2][1],1],[xy[3][1],1],[xy[4][1],1],[xy[5][1],1],[xy[6][1],1],
[xy[7][1],1],[xy[8][1],1],[xy[9][1],1],[xy[10][1],1])$
Y:matrix([xy[1][2]],[xy[2][2]],[xy[3][2]],[xy[4][2]],[xy[5][2]],[xy[6][2]],
[xy[7][2]],[xy[8][2]],[xy[9][2]],[xy[10][2]])$
invert(transpose(M).M).transpose(M).Y;