pianofisica

Mathematics & Physics, Maxima, a bit Python & Wolfram, and Arts

Maximaで学ぶ最小二乗法

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

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

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

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

平均・分散・標準偏差

2つのデータからなる  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}

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

フィッティング

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

 \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;

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

データ点と1次関数を

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


最小二乗法(統計の言葉で)

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

 \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 について解いて求められたものが上で与えた1次関数です。

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

 \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;

より

 \displaystyle{\ X=\left(\begin{array}{c} A \\ B \end{array}\right)=\left(\begin{array}{c} \frac{338}{501} \\ \frac{751}{501} \end{array}\right)}

を得ます。これはもちろん、統計の言葉で表した場合に求めた結果と一致します。

Maximaでの行列の取り扱いについては過去の記事

pianofisica.hatenablog.com

pianofisica.hatenablog.com

を参考にしてみてください。



キーワードMaxima、最小二乗法、統計学、平均、分散、標準偏差、共分散