最小二乗法による係数の推定

本日の講義概要

線形最小二乗法とは
MATLABの演算子

線形最小二乗法とは

問題の考え方

ある系列の２つの測定値が得られる場合に，この２つの測定値の関係を求める問題を考える．

例：

ある学級の生徒の身長と体重の値
実験実習で行った，金属柱の長さ及び直径の測定値と質量

系列	1	2	3	4	…	n
x	x1	x2	x3	x4	…	xn
y	y1	y2	y3	y4	…	yn

このときに，２つの系列の測定値の関係が，線形の関係にあると考える．

ｙを表現すると仮定する関数は基準周波数の整数倍の周波数を持つ正弦波でも良いし，より自由度の高い自由曲線を考えることもできる．

身長体重の単純な例では，

で表現されていると考えることができる．

すなわち，最小二乗法は，「測定されたｘとｙの系列から，ある条件を満たす係数列ａを求める問題」と考えることができる．

用語の定義

誤差: 測定値と真の値（永遠に不明）の差，従って具体的な誤差値も永遠に不明となる．通常の測定では目盛の10分の1を目視で読むので，この最小の桁に誤差が含まれる，として扱う．従って，誤差値に基づいて何かを行うことは出来ない．
例：一目盛が1mmの定規で長さを測定する場合は，0.1mmまで測定して，測定値が4.5mmとする．このとき真の値は，4.45mm～4.54mmの間にあることになり，誤差の範囲は，±0.05mmとなる．
残差: 推定値と測定値との差．この値は計算することが出来る．->最小二乗法で用いることができる．

線形最小二乗法の定義

最小二乗法では残差の二乗和が最小になるように各係数を求める方式である．具体的には，残差二乗和を各係数で偏微分した式を=0とした連立方程式を立てて解く．すなわち残差二乗和Sを

として

これらをa0及びa1の連立方程式として解く．すると

aは求めるべき係数，X及びｙは実験値から求まる，この２つからaを求めることが出来る．上の式のaにかかる行列を消去してｙ側に移せれば良い．Matlabでは，「左からの商」と呼ばれる演算子（\:日本語環境では円マーク，英語環境ではバックスラッシュ）が存在し，これを使うことによってaの値をダイレクトに求めることができる．

この計算を行うMTALABスクリプトは以下のようになる．MATLABらしいプログラムになっているので，プログラムのコメントと後ろの解説を良く読んでほしい．

注1：WWWブラウザからカットアンドペーストを行うと全角スペースが挿入されてシンタクスエラーとなることがあるので注意すること．
注2：作成した関数をセーブする際には，関数名と同じファイル名とすること！

% function lsm1()

% function for least square method 1

% Input: none

% Output: Graph out

% By your_name, your class, build date

% function y=lsm1()

%x,y をそれぞれ15行1列として定義．アポストロフィは転地行列を求める演算子

x=[73 72 56 77 64 71 81 32 47 79 33 92 83 67 55]'

y=[82 70 50 82 66 75 89 40 58 70 25 89 87 73 78]'

n=length(x) 　　　　　　 %次の行で，x行列と同じ行数で1列，かつ内容が全て1の行列を作るためにxの行数を求める

X=[x ones(n,1)] 　　　　　　 %15行2列の行列で1列目がｘ，2列目に全て1を入れる

a=X\y　　　　　　　　　　 % a=X\yは，X*a=yの解，即ち左からの商

yy=a(1).*x+a(2)　　 %求めたaと測定値ｘからyの推定値yyを求める

res=y-yy　　　　　　　　　　　%測定値yと推定値yyの差＝残差

plot(x,y,'o',x,yy,'-'), axis equal

練習問題：

各変数の内容を確認すること

補足

MATLAB特有の演算子と演算結果を示しておく．本来\は，バックスラッシュ記号だがSJISのコード体系中では￥記号で表示されるので注意すること．x=[1 2 3 4]; y=[5 6 7 8];として

*は行列の乗算	z1=x*y' ※ｘ，ｙのままでは乗算できないのでｙを転置している．	z1=70
\は行列の左除算	z2=x\y	z2 = 0 0 0 0 0 0 0 0 0 0 0 0 1.2500 1.5000 1.7500 2.0000

.*は，要素ごとの積．ｘとｙの両方が配列の時は，両方が同じ大きさでなくてはならない．	z1=x..*y	z1 =5 12 21 32
\ は，左からの商．従ってz2は，x*z2=y の解	z2=x.\y	z2 =5.0000 3.0000 2.3333 2.0000
.^は要素ごとのべき乗	z3=x.^y	z3 =1 64 2187 65536
	z4=x.^2	z4 =1 4 9 16
	z5=2.^[x y]	z5 =2 4 8 16 32 64 128 256

関数を2次多項式として考えると

係数は3個になるので，３元連立方程式を解くことになる．

練習問題　

では，1次関数の最小二乗法のプログラムを修正して，2次関数の最小二乗法を行ってみよう．

１次関数の場合と異なっている部分がどこか？
どのような演算子を使えばよいか？ヒントは，前節の補足に

単純に拡張しただけではグラフの出力がおかしくなるので，連続的な直線になるようにプログラムを修正せよ．ヒント：plot関数に渡す値を修正する必要がある．

数値解析２のページへ