線形回帰のプログラミング

線形回帰（Linear Regression）を理解するためにプログラム（Processingを使用）を組んでみました。相変わらずこのような数学的な内容はWikipediaの説明ではわかりにくいので（もうすでに理解している人のために書いている内容であって、分からない人が理解できるような説明にはなっていない）、プログラミングすれば理解が深まるかもしれないと思って試してみました。

センサーを使ってデータを集めた際に、以下のようなグラフになることがあります。これらのドット分布を簡単な数式（今回の場合は単純な直線の式：赤い直線）に置き換えることができないかということです。

直線の一次式は、

y = m*x + b

で表せるので、要はmの傾きとbの値（x=0の時のyの値）を求めることができればいいというわけです。
どうやら公式があるようで、

m = Σ(Xi - Xmean) * (Yi - Ymean) / Σ(Xi - Xmean) *(Xi - Xmean)
b = Ymean - m*Xmean

となるようです（Σはn=0からn-1までの合計値）。
Xiはある点のXの値、XmeanはすべてのXの点の平均値、Yについても同様。
ただし、この式だとわかりにくいので、Σを使わないで表せば、

それぞれのX座標値：X0、X1、X2・・・Xn-1
それぞれのY座標値：Y0、Y1、Y2・・・Yn-1

とすれば（n個の点があるという前提）、
プログラム上のfor文では、for(int i=0; i<n; i++)となるのでn-1まで、

(X0-X平均値)*(Y0-Y平均値)+(X1-X平均値)*(Y1-Y平均値)+(X1-X平均値)*(Y1-Y平均値)・・・
+(Xn-1-X平均値)*(Yn-1-Y平均値)

の合計と、

(X0-X平均値)*(X0-X平均値)+(X1-X平均値)*(X1-X平均値)+(X2-X平均値)*(X2-X平均値)・・・
+(Xn-1-X平均値)*(Xn-1-X平均値)

の合計を求めて、
上の合計を下の合計で割ればmが求まり、
bについては、

b = Y平均値 - m*X平均値

で直線の式がでてくるようです。
ある式をn=0からi=n-1まで計算していくには、プログラミングではfor文を使えばいいので簡単です。
ただし、方法としてはこれもまたよくでてくる最小二乗法というものです。この最小二乗法が解法のテクニックだと思うのですが、求めたい直線と各点のyの値とのズレをすべて合計して（さらにプラスマイナスの符号をなくすために二乗する）、その差を少なくすればいいというのは直感的に分かるので、そのあたりをプログラムを通して理解しようとしてみました。
複数ある点の平均値を求めるので、そこが点分布の重心となり、かならずこの重心を通るような線となるようです。

空白のグラフ（上）がありますが（上に画面が現れない場合や表示が変な場合はこちらをクリック）、任意の位置をクリックすると点が追加されていきます（Processing/Processing.jsで書いています）。
複数点を打っていくと、その分布する点に対する一次式の直線が計算により求められるというプログラムです。当然二点打てば、その二点を通る直線ができるわけですが、三点目を少しずらしたところに打てば、その三点以上の点分布に対する直線式が出てきます。異なる位置に点を追加していくことで、微妙に線の傾きが変化していきます。
直線は点分布の重心を必ず通ることから、二点打ったあとにその中点を打ち、その上下（y軸と平行に）に等間隔でさらに二点打てば、傾きは変化しないということになります。こうやって、点を追加していくことで、傾きが変化したりしなかったりという性質が分かるので理解が深まります。
＊リセット機能はプログラムに含まれていないので、ブラウザでリロードすれば初期化されます。

以下が、Processingによるコード。

//クリックで追加されていく点の座標の配列を用意
int[] x = new int[0]; 
int[] y = new int[0];

void setup(){
  size(400,400);
}

void draw(){
  background(51);

  //xとyの平均値を求める
  int xsum = 0; //xの合計値の変数
  int ysum = 0; //yの合計値の変数
  fill(200);
  stroke(200);
  for(int i = 0; i < x.length; i++){//点の数だけ合計値を求める
    xsum += x[i];
    ysum += y[i];
    ellipse(x[i],y[i],4,4);//各点の描画
  }
  float xmean = 1.0 * xsum / x.length; //xの平均値
  float ymean = 1.0 * ysum / y.length; //yの平均値

  //公式を使って計算する
  float xy = 0; //公式の分子の変数
  float xx = 0; //公式の分母の変数
  for(int i = 0; i < x.length; i++){
    xy += (x[i]-xmean)*(y[i]-ymean); //分子の加算
    xx += (x[i]-xmean)*(x[i]-xmean); //分母の加算
  }
  
  float m = xy / xx; //上記公式からmを求める
  float b = ymean - m * xmean; //公式で求めたmと各平均値からbを求める
  float x1 = 0; //グラフ左端x=0のとき
  float y1 = m * x1 + b; //x=0のときのyの値、結果y1= bとなる
  float x2 = width;  //グラフ右端x=画面幅のとき
  float y2 = m * x2 + b; //x=画面幅のときのyの値、mとbは求まっているのでy2も求まる
  stroke(255,0,0);
  line(x1,y1,x2,y2); //直線を描く
}

//マウスクリックするごとに配列に新たな座標を加える
void mousePressed(){ 
  x = append(x, mouseX); 
  y = append(y, mouseY);
}

プログラミングにおけるY座標は下向きに+ですが、mにマイナスをかけて、height-bとすれば、Y座標上向きを+にすることができます（今回は数値表現していないため省略）。

ネットなどに書いてある線形回帰の説明だと数式だけで数学的に記述されているので、わかりにくいのですが、こうやってプログラミングすることで、どのように計算しているのか理解が深まるという感じです。
もう一つ、勾配法による求め方もあるので、次回やってみようと思います。
以前、A*経路探索アルゴリズムも試してみましたが、このあたりのテクニック/アルゴリズムは、いろいろな場面で使えるので、知っていれば便利かと思って試している最中です。

続き：
線形回帰のプログラミング（その２）：勾配降下法