最小二乗法

MATLABでは、マウスやカーソルから(x、y)のデータポイントを入力することができます。MATLABコードをワードプレスに変換するオンラインツールがありました。http://cs.unm.edu/~drew/matlab_converter/ これを使ってみることにします。

 [x,y] = ginput(10])

グラフが現れるので、ポチッ、ポチッ、っと10回クリックすると、x、yに10個点が入ります。なんとなく直線に乗るように打ってみました。

x =[

    0.5158
    0.5972
    0.7040
    0.8502
    0.7055
    0.4597
    0.2802
    0.1055
    0.0447
    0.8217
]
y =[

    0.5346
    0.6309
    0.7475
    0.9434
    0.7323
    0.5853
    0.3218
    0.2019
    0.0887
    0.9029]

グラフで確認。だいたい直線上。

plot(x,y,'o')

データポイントy_i – (a*x_i+b)の2乗の和を最小にする条件を求める。係数a,b でそれぞれ偏微分したものがゼロになるという条件を使うと、いわゆる正規方程式が導かれる。その正規方程式を行列を用いてかけば、

A*a=v (変数の記号や中身は、下のMATLABコードを参照。ややこしいがa=[うえのカーブフィットの直線の式の係数a;係数 b])

MATLABをつかうと、係数を求めるのは一瞬でできてしまう。

A=[sum(x.^2) sum(x); sum(x) 1];
v=[sum(x.*y);sum(y)];

% A*[a b]=v  を解けばよい。それにはAの逆行列を左からかけて、
inv(A)*v
ans =

    1.1228
   -0.0198

a=inv(A)*v  ややこしいがa=[係数a;係数 b]のつもり。

直線で近似できたかグラフを描いて確かめてみる。

a=inv(A)*v;
t=a(1)*s+a(2);
plot(s,t);
hold on;
plot(x,y,'o');
axis([min(x) max(x) min(y) max(y)]);

できた。

簡単な計算でも符号を間違ったり、簡単なMATLABコードなのに変数を逆にしてたり足し算と掛け算を間違えたり、いろいろやらかしてしまったので、期待される結果のグラフが描けたときには爽快感がありました。データにフィットする直線を描いたり、その式を出すのは、エクセルをつかうとボタン一つなわけですが、数学的な中身を理解するというのは、単純に気持ちのいいものです。それにしても、数式混じりの本は、読むのと、やるのとは大違いです。

参考図書:これならわかる応用数学教室 第1章 最小二乗法

この本は岡山大学工学部の数学未履修の学生向けの講義がもとになってできた教科書なので、数学に苦手意識がある人にお勧め。こまごましたことは大胆にはしょって、数学的な面白さはしっかり教えてくれる本。最終目標が画像処理やウェーブレットの入門あたりになっており、それに必要な線形代数の基礎事項が説明されています。学生と先生との対話がふんだんに盛り込まれていて、面白いです。学生が、「先輩に、情報系だから数学は要らないといわれたので解析学の授業はとらなかったのですが。」みたいな寝ぼけた発言をしたりして、先生はそんな学生にも根気強く数学を説明しています。