4月 282010
 

回帰分析(最小二乗法)を理解しているか?

データ解析の初歩的な話をしておこう.工学部の学生や卒業生であれば,最小二乗法でデータから直線を求めるなんてことは大学1回生のころから実験データの整理などでやっているはずなのだが,まともに理解しているのはほとんどいない.データ解析の産業応用を幅広く手掛けていることもあって,多くの企業で出張講義(?)をさせていただくのだが,企業のエンジニアも理解していない.考えてみれば,それもそのはずで,大学では最小二乗法(または回帰分析)をまともに教えてもらう機会がないわけだ.いきなり,もっと難しい話に飛んでしまう.

回帰分析(最小二乗法)で直線の傾きを求める

さて,説明変数(入力変数) x と目的変数(出力変数) y の測定データが与えられたとき, y=ax+b という直線でデータを近似することを考えよう.例えば,100人の身長 x と体重 y を計測して,身長から体重を計算する式 y=ax+b を構築することを考える.普通の教科書だと,ここで b をどうやって求めるかという話が来るのだが,そんなものは x と y 両方を平均 0 にしてしまえば無視できる.つまり, b=0 となる.ついでに, x と y 両方の標準偏差も 1 にしてしまおう.つまり, x と y 両方について,平均 0 かつ標準偏差(および分散) 1 とする.この操作は, x と y それぞれから平均を引いて,さらに標準偏差で割ればいいだけだ.以下では,説明変数(入力変数) x と目的変数(出力変数) y は共に平均 0 かつ標準偏差(および分散) 1 とする.このとき,問題は y=ax という直線の傾き a を求めることになる.

子供から大人まで様々な年齢の人達の身長 x と体重 y を測定したとしよう.そうすれば,身長 x と体重 y は強い正の相関を持つと予想される.つまり,身長が高ければ,体重は重いだろうというわけだ.このとき,回帰分析(最小二乗法)で傾き a を求めるとどうなるだろうか.身長 x と体重 y は共に平均 0 かつ標準偏差(および分散) 1 であるため, x-y 平面に測定データをプロットして散布図を描けば,データは斜め45度(傾き1)の直線の周りに分布する.ここまではいいだろうか.

多くの人はここで,回帰分析(最小二乗法)で求まる傾きは a=1 だと答える.自信を持って.なぜなら,データが斜め45度(傾き1)の直線の周りに分布するからだ.

しかし現実には,そうはならない.

簡単な計算で確かめることができるが,回帰分析(最小二乗法)で求まる傾きは, x と y の相関係数に等しくなる.したがって,傾きが1となるのは,すべてのデータが綺麗に一直線上にのったときだけだ.身長と体重の測定データにおいては,そんなことはありえない.散布図の名にふさわしく,データはばらつく.相関係数が正であることは間違いないだろうが,相関係数は1にはならない.つまり,現実に求まる直線は,妄想した直線よりも,もっと寝た(x軸側に偏った)直線にしかならない.

傾きが1にならないことを確かめる

このことを実際に数値例で確認してみよう. y = x + noise という式を使って,100点の測定データを発生させる.説明変数 x は正規分布に従うとし, noise も正規分布に従うとする.ここで, noise の大きさを色々と変えてみる.4種類の大きさの noise について実際にデータを発生させて,回帰分析(最小二乗法)を用いて直線を求めた結果を下図に示す.赤線が回帰分析(最小二乗法)の結果だ.

バラツキが大きくなると,明らかに直線は寝てしまう.傾きは相関係数に一致するのだから当然だ.図中の青線が求まると思っていた人はいないか?

回帰分析(最小二乗法)と主成分分析の傾きの違い
回帰分析(最小二乗法)と主成分分析の傾きの違い

傾き1の直線は求められるか?

さて,ここで注目して欲しいのは,図中の青線は傾き1の直線を描いたわけではなくて,100点の測定データから統計的手法によって求めた直線だということだ.回帰分析(最小二乗法)よりも,ずっとデータの分布を表現しているように思える.この青線はどうしたら求められるのだろうか.

実は,主成分分析を使っている.つまり,青線は第一主成分だ.主成分分析を知っている人なら,あぁ,そうか,と納得できるだろう.知らない人は仕方がない.主成分分析のテキストでも読んで勉強してみよう.損はしない.

産業界が抱える問題

問題は,このようなことも知らないで,産業界の実務で回帰分析(最小二乗法)が利用されていることだ.そのような状況で,まともな結論が導出できるとも思えない.もちろん回帰分析(最小二乗法)だけの所為ではないが,結果として,「データ解析は使えない」とか,「データ解析なんてただの数字遊びだ」とか,そんな戯れ言が罷り通るようになってきたのではないだろうか.いやいや,データ解析が悪いのではなくて,使い方が悪いんでしょと言いたくもなる.ところが,これを自信を持って言えるエンジニアが産業界に少ない.本当に少ない.最初の話に戻るが,教えられていないのだから当然とも言える.

ベイズだカーネルだと色々あるけれども,品質不良と戦う現場では,品質改善を目指して,膨大な要因を相手に,散布図を描いたり,相関係数を求めたりしていることだろう.その現場に,正しい知識を持ち込む必要がある.取り組むべき課題の1つだと思われる.

付録

今回の話は,このブログに書くような内容ではない気もするが,気になっていることでもあるので,書いてみた.私が担当する「プロセスデータ解析学」を受講すれば,こういう初歩的なことも習うわけだが,これも大学院向けの講義だ.学部では習わない.最小二乗法なんて,工学リテラシー(?)として必須だと思うのだが...

それはともかく,自分で確認してみたいという人のために,以下にMATLABのプログラムを記載しておこう.上図を作成するために使用したプログラムだ.簡単な短いプログラムなので,見ればわかるだろう(と期待する).


ndata = 100;

% explanatory variable
x = randn(ndata,1);
x = (x-repmat(mean(x),size(x)))./repmat(std(x),size(x));

% objective variable
noise = randn(ndata,1)*3;
y = x + noise;
y = (y-repmat(mean(y),size(y)))./repmat(std(y),size(y));

% multiple regression analysis
C_MRA = (x'*x)\(x'*y) % slope calculated by MRA

% principal component analysis
[U,S,V] = svd([x y]);
C_PCA = V(2,1)/V(1,1) % slope calculated by PCA

figure(1)
plot(x,y,'k.'),grid,hold on
y_MRA = C_MRA*[-4 4];
y_PCA = C_PCA*[-4 4];
plot([-4 4],y_MRA,'r-',[-4 4],y_PCA,'b-'),hold off
xlabel('x','fontsize',12)
ylabel('y','fontsize',12)
axis([-4 4 -4 4])
set(gca,'fontsize',12)
legend('data','MRA','PCA')

ちなみに,最初に触れた出張講義(?)は随時受け付けています.スケジュールが確保できればですが...