C言語で線形回帰（単回帰）・最小二乗法

こんにちは。

今回は、線形回帰（単回帰）を最小二乗法を使って行う方法について書いていきます。

単回帰は回帰の中でも、簡単なので覚えるといいと思います。

それでは書いていきます。

線形回帰（単回帰）の考え方

前提

まず、以下のように実験などでデータが取れたとします。

この時に、なんとなく線が引けそうなので、線が引きたいと思います。

単回帰では、その線の式を以下のように考えます。

$$y=ax+b$$

見てわかるように、単純な一次関数ですね。

式を見ると、$a$と$b$が求まれば、線の式がわかりそうです。

最小二乗法

最小二乗法では、上の式の$a$と$b$を各点との誤差が小さくなるように求めます。

例えば、下の図の緑矢印の大きさが誤差にあたります。

今回は上図のようにy軸方向の誤差を最小とするように考えます。

ここで、赤点の位置を $(x,y)=(x_k,y_k)$ とすると、
その点と $y=ax+b$の誤差の大きさは次のようになります。

$$y_k-ax_k+b$$

このままでは、この値が正の時と負の時があって扱いずらいので、二乗して正のみになるようにします。

$$(y_k-ax_k+b)^2$$

また、これは一点のみについての話なので、仮に赤点が全部でN点あったとすると、それらの誤差の2乗の和は以下のようになります。

$$\sum_{i=1}^{N}{(y_i-ax_i+b)^2}$$

結局は、この式の値を最小化するということになります。

上の式を展開すると、

$$\begin{aligned} &\sum\left(y_{i}-a x_{i}-b\right)^{2} \\ &=a^{2} \sum x_{i}^{2}+n b^{2}+\sum y_{i}^{2} -2 a \sum x_{i} y_{i}-2 b \sum y_{i}+2 a b \sum x_{i} \end{aligned}$$

となります。

この式が最小になるとき、 $a$, $b$ での偏微分の値が0となるので、
$a$, $b$ で偏微分して、

$a$ で偏微分 : $2 a \sum x_{i}^{2}-2 \sum x_{i} y_{i}+2 b \sum x_{i}=0$
$b$ で偏微分 : $2 n b-2 \sum y_{i}+2 a \sum x_{i}=0$

これを、連立方程式とみなして、解くと以下のようになります。

$$a=\frac{\sum_{i=1}^{n}\left\{\left(x_{i}-\bar{x}\right)\left(y_{i}-\bar{y}\right)\right\}}{\sum_{i=1}^{n}\left(x_{i}-\bar{x}\right)^{2}}$$ $$b=\bar{y}-a \bar{x}$$

ここで、 $\bar{x}$, $\bar{y}$ はそれぞれ $x$, $y$ の平均値としています。

小まとめ

• $(x,y)$ のように点がたくさんあって、それに直線を引きたい
• 直線を引くには、誤差の二乗の和を最小化する（最小二乗法）
• 誤差の二乗の和を最小化するには、偏微分をする
• 結果として、以下のようになる

$$a=\frac{\sum_{i=1}^{n}\left\{\left(x_{i}-\bar{x}\right)\left(y_{i}-\bar{y}\right)\right\}}{\sum_{i=1}^{n}\left(x_{i}-\bar{x}\right)^{2}}$$ $$b=\bar{y}-a \bar{x}$$

C言語でコードを書く

今回は以下のような点を与えます。

x	1	2	3	4	5	6	7	8	9	10
y	1.16	5.07	5.22	6.59	10.92	11.52	13.11	17.09	16.66	20.6

Excelでみると、以下のようになります。
Excelの近似曲線で入れた式も出しています。

コードは以下のように書きました。

#include <stdio.h>

int main(void){
  // 入力点
  double xs[10] = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10};
  double ys[10] = {1.16, 5.07, 5.22, 6.59, 10.92, 11.52, 13.11, 17.09, 16.66, 20.6};

  // xとyの平均を求める
  double x_sum = 0.0;
  double y_sum = 0.0;
  for(int i=0; i<10; i++){
    x_sum = x_sum + xs[i];
    y_sum = y_sum + ys[i];
  }
  double x_ave = x_sum / 10.0;
  double y_ave = y_sum / 10.0;

  // aの値を求める
  double a_top = 0.0;
  double a_bottom = 0.0;
  for(int i=0; i<10; i++){
    a_top = a_top + (xs[i]-x_ave)*(ys[i]-y_ave);
    a_bottom = a_bottom + (xs[i]-x_ave)*(xs[i]-x_ave);
  }
  double a = a_top / a_bottom;

  // bの値を求める
  double b = y_ave - a*x_ave;

  // aとbの結果を表示する
  printf("a: %lf\n", a);
  printf("b: %lf\n", b);
}

これを実行すると、以下のような結果が得られます。

a: 2.033939
b: -0.392667

Excelで実行した結果とあっていることが確認できました。

今回はこれで終わります。