最小2乗法

最小２乗法は、ある変量の組(説明変数)とその変量に対する望みの結果(目的 変数、教師信号)が学習データとして与えられた時、説明変数から目的変数を 予測するモデルを構築するための統計的データ解析手法で、最も基本的で、最 も広く用いられています。ここでは、以下のような例題に対するプログラムを 作ってみましょう。

問題2.

15人の中学生のボール投げの記録($t$)と握力($x1$)、身長 ($x2$)、体重($x3$)のデータがあります。

表 1: 中学生のボール投げの記録($t$)と握力($x1$)、身長($x2$)、体重 ($x3$)のデータ

生徒番号	ボール投げ(m)	握力(kg)	身長(cm)	体重(kg)
1	22	28	146	34
2	36	46	169	57
3	24	39	160	48
4	22	25	156	38
5	27	34	161	47
6	29	29	168	50
7	26	38	154	54
8	23	23	153	40
9	31	42	160	62
10	24	27	152	39
11	23	35	155	46
12	27	39	154	54
13	31	38	157	57
14	25	32	162	53
15	23	25	142	32

このデータを用いて、握力($x1$)、身長($x2$)、体重($x3$)からボール投げの 記録($t$)を予測するための線形モデル

$$y(x1,x2,x3) = a_0 + a_1 x1 + a_2 x2 + a_3 x3$$ (9)

を求めなさい。

つまり、中学生のボール投げの記録($t$)と握力($x1$)、身長($x2$)、体重 ($x3$)のデータを $\{<t_l,x1_l,x2_l,x3_l>\vert l=1,\ldots,15\}$とすると、$l$ 番目の生徒のボール投げの記録を

$$y(x1_l,x2_l,x3_l) = a_0 + a_1 x1_l + a_2 x2_l + a_3 x3_l$$ (10)

で予測します。このモデルには、$a_0,a_1,a_2$および$a_3$の4個のパラメー タが含まれています。これらのパラメータを学習用のデータから決めなければ なりません。最小２乗法では、予測のための線形モデル良さの評価基準として、 望みの結果とモデルが予測した結果との誤差の2乗の期待値（平均２乗誤差） を用い、２乗誤差が最も小さくなるようなパラメータを探索します。今考えて いる例題では、説明変数の組 $<x1_l,x2_l,x3_l>$に対する望みの結果が $t_l$ で、モデルの出力が $y_l$ ですので、その誤差 $(t_l - y_l)$ の2乗 の平均（平均2乗誤差） $\varepsilon^2(a_0,a_1,a_2,a_3)$ は、

$$\varepsilon^2(a_0,a_1,a_2,a_3) = \frac{1}{15} \sum_{l=1}^{15} \varepsilon_l^2 = \frac{1}{15} \sum_{l=1}^{15} (t_l - y_l)^2$$

$$= \frac{1}{15} \sum_{l=1}^{15} \{t_l - (a_0 + a_1 x1_l + a_2 x2_l + a_3 x3_l)\}^2$$ (11)

のようになります。

この問題では、評価関数 $\varepsilon^2(a_0,a_1,a_2,a_3)$は、各パラメータ に関して2次の関数ですので、以下のような方法で最適なパラメータを求める ことが可能です。

最小２乗法の解法

式(11)のパラメータ$a_0$に関する微分 を求めると

$$\frac{\partial \varepsilon^2}{\partial a_0} = - 2 \{ (\frac{1}{15} \sum_{l=1}^{15} t_l) - a_0 (\frac{1}{15} \su... ...\frac{1}{15} \sum_{l=1}^{15} x2_l) - a_3 (\frac{1}{15} \sum_{l=1}^{15} x3_l) \}$$

$$= - 2 \{ \bar{t} - a_0 - a_1 \bar{x1} - a_2 \bar{x2} - a_3 \bar{x3} \}$$ (12)

のようになります。最適な解は、この値が$0$となることが必要ですから、 $\frac{\partial \varepsilon^2}{\partial a_0}=0$とおくと、

$$a_0 = \bar{t} - a_1 \bar{x1} - a_2 \bar{x2} - a_3 \bar{x3}$$ (13)

が得られます。ここで、$\bar{t}$, $\bar{x1}$, $\bar{x2}$および $\bar{x3}$は、それぞれ、$t$, $x1$, $x2$および$x3$の平均値で、

$$\bar{t} = \frac{1}{15} \sum_{l=1}^{15} t_l$$

$$\bar{x1} = \frac{1}{15} \sum_{l=1}^{15} x1_l$$

$$\bar{x2} = \frac{1}{15} \sum_{l=1}^{15} x2_l$$

$$\bar{x3} = \frac{1}{15} \sum_{l=1}^{15} x3_l$$ (14)

のように定義されます。今、これを、式(10)に代入すると

$$y(x1_l,x2_l,x3_l) = \bar{t} + a_1 (x1_l - \bar{x1}) + a_2 (x2_l - \bar{x2}) + a_3 (x3_l - \bar{x3})$$ (15)

のようになります。従って、平均２乗誤差は、パラメータ$a_1$, $a_2$ およ び $a_3$の関数として、

$$\varepsilon^2(a_1,a_2,a_3) = \frac{1}{15} \sum_{l=1}^{15} (t_l - y_l)^2$$

$$= \frac{1}{15} \sum_{l=1}^{15} \{t_l - \bar{t} - a_1 (x1_l - \bar{x1}) - a_2 (x2_l - \bar{x2}) - a_3 (x3_l - \bar{x3}) \}^2$$ (16)

のように書けます。今、この平均２乗誤差をパラメータ $a_1$ で微分すると

$$\frac{\partial \varepsilon^2}{\partial a_1} = - \frac{1}{15} \sum_{l=1}^{15} \{t_l - \bar{t} - a_1 (x1_l - \bar{x1}) - a_2 (x2_l - \bar{x2}) - a_3 (x3_l - \bar{x3})\}\{x1_l - \bar{x1} \}$$

$$= - \{\sigma_{t1} - a_1 \sigma_{11} - a_2 \sigma_{21} - a_3 \sigma_{31}\}$$ (17)

のようになります。同様に、平均２乗誤差をパラメータ $a_2$ および $a_3$ で微分すると、

$$\frac{\partial \varepsilon^2}{\partial a_2} = -\{ \sigma_{t2} - a_1 \sigma_{12} - a_2 \sigma_{22} - a_3 \sigma_{32} \}$$ (18)

$$\frac{\partial \varepsilon^2}{\partial a_3} = -\{ \sigma_{t3} - a_1 \sigma_{13} - a_2 \sigma_{23} - a_3 \sigma_{33} \}$$ (19)

のようになります。ここで、

$$\sigma_{11} = \frac{1}{15} \sum_{l=1}^{15} (x1_l - \bar{x1})(x1_l - \bar{x1}),$$

$$\sigma_{12} = \frac{1}{15} \sum_{l=1}^{15} (x1_l - \bar{x1})(x2_l - \bar{x2}),$$

$$\sigma_{13} = \frac{1}{15} \sum_{l=1}^{15} (x1_l - \bar{x1})(x3_l - \bar{x3}),$$

$$\sigma_{21} = \frac{1}{15} \sum_{l=1}^{15} (x2_l - \bar{x2})(x1_l - \bar{x1}),$$

$$\sigma_{22} = \frac{1}{15} \sum_{l=1}^{15} (x2_l - \bar{x2})(x2_l - \bar{x2}),$$

$$\sigma_{23} = \frac{1}{15} \sum_{l=1}^{15} (x2_l - \bar{x2})(x3_l - \bar{x3}),$$

$$\sigma_{31} = \frac{1}{15} \sum_{l=1}^{15} (x3_l - \bar{x3})(x1_l - \bar{x1}),$$

$$\sigma_{32} = \frac{1}{15} \sum_{l=1}^{15} (x3_l - \bar{x3})(x2_l - \bar{x2}),$$

$$\sigma_{33} = \frac{1}{15} \sum_{l=1}^{15} (x3_l - \bar{x3})(x3_l - \bar{x3}),$$

$$\sigma_{t1} = \frac{1}{15} \sum_{l=1}^{15} (t_l - \bar{t})(x1_l - \bar{x1}),$$

$$\sigma_{t2} = \frac{1}{15} \sum_{l=1}^{15} (t_l - \bar{t})(x2_l - \bar{x2}),$$

$$\sigma_{t3} = \frac{1}{15} \sum_{l=1}^{15} (t_l - \bar{t})(x3_l - \bar{x3})$$ (20)

です。これらは、分散あるいは共分散と呼ばれています。また、これらの定義 から、

$$\sigma_{12} = \sigma_{21}, \hspace*{2mm} \sigma_{13} = \sigma_{31}, \hspace*{2mm} \sigma_{23} = \sigma_{32}$$ (21)

のような関係が成り立つことが分かります。

最適なパラメータでは、平均２乗誤差のパラメータに関する微分が $0$ となる はずですので、それらを $0$ とおくと、結局、

$$a_1 \sigma_{11} + a_2 \sigma_{12} + a_3 \sigma_{13} = \sigma_{t1}$$

$$a_1 \sigma_{21} + a_2 \sigma_{22} + a_3 \sigma_{23} = \sigma_{t2}$$

$$a_1 \sigma_{31} + a_2 \sigma_{32} + a_3 \sigma_{33} = \sigma_{t3}$$ (22)

のような連立方程式が得られます。従って、この連立方程式を解くためのサブ ルーチンがあれば、最適なパラメータが求まることになります。この連立方程 式を行列とベクトルを使って表現すると、

$$\Sigma \mbox{\boldmath$a$} = \mbox{\boldmath$\sigma$}$$ (23)

のようになります。ここで、行列 $\Sigma$、およびベクトル $\mbox{\boldmath$a$}$, $\mbox{\boldmath$\sigma$}$は、

$$\Sigma = \left( \begin{array}{ccc} \sigma_{11} & \sigma_{... ...ma_{t1} \\ \sigma_{t2} \\ \sigma_{t3} \end{array} \right)$$ (24)

のように定義されます。行列 $\Sigma$ は、分散共分散行列と呼ばれています。 もし、この行列 $\Sigma$ が正則で逆行列$\Sigma^{-1}$が存在するなら、上 の連立方程式の両辺に左から逆行列$\Sigma^{-1}$をかけて、

$$\mbox{\boldmath$a$} = \Sigma^{-1} \mbox{\boldmath$\sigma$}$$ (25)

のように最適なパラメータが求まります。今考えている問題では、分散共分散 行列 $\Sigma$ は、$3 \times 3$ の行列で、その逆行列は、逆行列の公式から

$$\Sigma^{-1} = \frac{1}{\vert\Sigma\vert} \left( \begin{... ... \sigma_{11} \sigma_{22} - \sigma_{12}^2 \end{array} \right)$$ (26)

のように求まります。ここで、$\vert\Sigma\vert$ は、行列 $\Sigma$ の行列式で、 $\vert\Sigma\vert = \sigma_{11} \sigma_{22} \sigma_{33} - \sigma_{11} \sigma_{23... ...\sigma_{33} - \sigma_{13}^2 \sigma_{22} + 2 \sigma_{12} \sigma_{13} \sigma_{23}$ です。この行列式$\vert\Sigma\vert$が $0$でない場合に逆行列が存在します。つまり、これが、最適なパラメータが 計算できる条件になります。

この式に従って、ボール投げの記録を予測するための最適なパラメータを具体 的に計算すると、

$$a_0 = -13.21730$$

$$a_1 = 0.20138$$

$$a_2 = 0.17103$$

$$a_3 = 0.12494$$ (27)

のようになります。以上ののようにして学習データから最適なパラメータが求 まれば、握力($x1=30$)、身長($x2=165$)、体重($x3=55$)のデータからの学生 のボール投げの記録を

$$y = -13.21730 + 0.20138 x 30 + 0.17103 x 165 + 0.12494 x 55 = 27.91575$$ (28)

のように予測することができるようになります。この問題の場合は、行列 $\Sigma$ が $3 \times 3$ でしたので、手計算でも最適解が求まりましたが、 一般には、連立方程式を解くサブルーチンを用いるなどして、最適解を求める 必要があります。

最急降下法によるパラメータの推定

次に、最小2乗法の最適なパラメー タを線形方程式を解かないで、最急降下法により求めるプログラムについて考 えてみます。最急降下法は、適当な初期パラメータからはじめて、パラメータ の値を微分値と逆の方向にちょっとだけ変化させて徐々に最適なパラメータに 近づけて行く方法ですので、まずは、評価関数(最小２乗法では、最小２乗誤 差)の各パラメータでの微分を計算する必要があります。

式(11)の最小２乗誤差のパラメータ$a_0$に関する微分は、

$$\frac{\partial \varepsilon^2}{\partial a_0} = - 2 \frac{1}... ... = - 2 \frac{1}{15} \sum_{l=1}^{15} (t_l - y(x1_l,x2_l,x3_l))$$ (29)

のように書けます。

同様に、最小２乗誤差のパラメータ$a_1$, $a_2$および$a_3$に関する微分は、

$$\frac{\partial \varepsilon^2}{\partial a_1} = - 2 \frac{1}{15} \sum_{l=1}^{15} \{ \varepsilon_l \frac{\partial ... ...psilon_l x1_l = - 2 \frac{1}{15} \sum_{l=1}^{15} (t_l - y(x1_l,x2_l,x3_l)) x1_l$$

$$\frac{\partial \varepsilon^2}{\partial a_2} = - 2 \frac{1}{15} \sum_{l=1}^{15} \{ \varepsilon_l \frac{\partial ... ...psilon_l x2_l = - 2 \frac{1}{15} \sum_{l=1}^{15} (t_l - y(x1_l,x2_l,x3_l)) x2_l$$

$$\frac{\partial \varepsilon^2}{\partial a_3} = - 2 \frac{1}{15} \sum_{l=1}^{15} \{ \varepsilon_l \frac{\partial ... ...psilon_l x3_l = - 2 \frac{1}{15} \sum_{l=1}^{15} (t_l - y(x1_l,x2_l,x3_l)) x3_l$$ (30)

のようになります。

従って、最急降下法による各パラメータの更新式は、

$$a_0^{(k+1)} = a_0^{(k)} - \alpha \frac{\partial \varepsilon^2}{\partial a_0}\ve... ...} = a_0^{(k)} + 2 \alpha \frac{1}{15} \sum_{l=1}^{15} (t_l - y(x1_l,x2_l,x3_l))$$

$$a_1^{(k+1)} = a_1^{(k)} - \alpha \frac{\partial \varepsilon^2}{\partial a_1}\ve... ..._1^{(k)} + 2 \alpha \frac{1}{15} \sum_{l=1}^{15} (t_l - y(x1_l,x2_l,x3_l)) x1_l$$

$$a_2^{(k+1)} = a_2^{(k)} - \alpha \frac{\partial \varepsilon^2}{\partial a_2}\ve... ..._2^{(k)} + 2 \alpha \frac{1}{15} \sum_{l=1}^{15} (t_l - y(x1_l,x2_l,x3_l)) x2_l$$

$$a_3^{(k+1)} = a_3^{(k)} - \alpha \frac{\partial \varepsilon^2}{\partial a_3}\ve... ..._3^{(k)} + 2 \alpha \frac{1}{15} \sum_{l=1}^{15} (t_l - y(x1_l,x2_l,x3_l)) x3_l$$ (31)

のようになります。

では、この式を用いて、握力($x1$)、身長($x2$)、体重($x3$)のデータからボー ル投げの記録($t$)を予測するための線形モデルのパラメータを最急降下法で 求めるプログラムを作成してみます。ただし、学習を安定化させるため各変数 ($x1$,$x2$,$x3$)の値を $100$ で割って、

$$x1' = \frac{x1}{100}, \hspace*{2mm} x2' = \frac{x2}{100}, \hspace*{2mm} x3' = \frac{x3}{100}$$ (32)

のように正規化してから利用します。これにより、学習係数をパラメータ毎に 変える必要がなくなると思います。

以下がそのプログラムの例です。

#include <stdio.h>

#define NSAMPLE 15
#define XDIM 3

main() {
  FILE *fp;
  double t[NSAMPLE];
  double x[NSAMPLE][XDIM];
  double a[XDIM+1];
  int   i, j, l;
  double y, err, mse;
  double derivatives[XDIM+1];
  double alpha = 0.2; /* Learning Rate */

  /* Open Data File */
  if ((fp = fopen("ball.dat","r")) == NULL) {
    fprintf(stderr,"File Open Fail\n");
    exit(1);
  }

  /* Read Data */
  for (l = 0; l < NSAMPLE; l++) {
    /* Teacher Signal (Ball) */
    fscanf(fp,"%lf", &(t[l]));
    /* Input input vectors */
    for (j = 0; j < XDIM; j++) {
      fscanf(fp,"%lf",&(x[l][j]));
    }
  }

  /* Close Data File */
  fclose(fp);

  /* Print the data */
  for (l = 0; l < NSAMPLE; l++) {
    printf("%3d : %8.2f ", l, t[l]);
    for (j = 0; j < XDIM; j++) {
      printf("%8.2f ", x[l][j]);
    }
    printf("\n");
  }

  /* scaling the data */
  for (l = 0; l < NSAMPLE; l++) {
    /*    t[l] = t[l] / tmean;*/
    for (j = 0; j < XDIM; j++) {
      x[l][j] = x[l][j] / 100.0;
    }
  }

  /* Initialize the parameters by random number */
  for (j = 0; j < XDIM+1; j++) {
    a[j] = (drand48() - 0.5);
  }

  /* Open output file */
  fp = fopen("mse.out","w");

  /* Learning the parameters */
  for (i = 1; i < 20000; i++) { /* Learning Loop */

    /* Compute derivatives */
    /* Initialize derivatives */
    for (j = 0; j < XDIM+1; j++) {
      derivatives[j] = 0.0;
    }

    /* update derivatives */
    for (l = 0; l < NSAMPLE; l++) {

      /* prediction */
      y = a[0];
      for (j = 1; j < XDIM+1; j++) {
	y += a[j] * x[l][j-1];
      }

      /* error */
      err = t[l] - y;
      /*      printf("err[%d] = %f\n", l, err);*/

      /* update derivatives */
      derivatives[0] += err;
      for (j = 1; j < XDIM+1; j++) {
	derivatives[j] += err * x[l][j-1];
      }

    }
    for (j = 0; j < XDIM+1; j++) {
      derivatives[j] = -2.0 * derivatives[j] / (double)NSAMPLE;
    }

    /* update parameters */
    for (j = 0; j < XDIM+1; j++) {
      a[j] = a[j] - alpha * derivatives[j];
    }

    /* Compute Mean Squared Error */
    mse = 0.0;
    for (l = 0; l < NSAMPLE; l++) {
      /* prediction */
      y = a[0];
      for (j = 1; j < XDIM+1; j++) {
	y += a[j] * x[l][j-1];
      }

      /* error */
      err = t[l] - y;

      mse += err * err;
    }
    mse = mse / (double)NSAMPLE;

    printf("%d : Mean Squared Error is %f\n", i, mse);
    fprintf(fp, "%f\n", mse);

  }

  fclose(fp);

  /* Print Estmated Parameters */
  for (j = 0; j < XDIM+1; j++) {
    printf("a[%d]=%f, ",j, a[j]);
  }
  printf("\n");

  /* Prediction and Errors */
  for (l = 0; l < NSAMPLE; l++) {
    /* prediction */
    y = a[0];
    for (j = 1; j < XDIM+1; j++) {
      y += a[j] * x[l][j-1];
    }
    
    /* error */
    err = t[l] - y;

    printf("%3d : t = %f, y = %f (err = %f)\n", l, t[l], y, err);
  }

}

このプログラムを実行させると、

a[0]=-13.156891, a[1]=20.077345, a[2]=17.056200, a[3]=12.562173,
  0 : t = 22.000000, y = 21.637957 (err = 0.362043)
  1 : t = 36.000000, y = 32.064105 (err = 3.935895)
  2 : t = 24.000000, y = 27.993037 (err = -3.993037)
  3 : t = 22.000000, y = 23.243744 (err = -1.243744)
  4 : t = 27.000000, y = 27.034110 (err = -0.034110)
  5 : t = 29.000000, y = 27.601042 (err = 1.398958)
  6 : t = 26.000000, y = 27.522622 (err = -1.522622)
  7 : t = 23.000000, y = 22.581754 (err = 0.418246)
  8 : t = 31.000000, y = 30.354062 (err = 0.645938)
  9 : t = 24.000000, y = 23.088664 (err = 0.911336)
 10 : t = 23.000000, y = 26.085890 (err = -3.085890)
 11 : t = 27.000000, y = 27.723396 (err = -0.723396)
 12 : t = 31.000000, y = 28.411174 (err = 2.588826)
 13 : t = 25.000000, y = 27.556856 (err = -2.556856)
 14 : t = 23.000000, y = 20.102145 (err = 2.897855)

のようなパラメータが求まり、学習に用いたデータに対して、求まったパラメー タを用いたモデルで予測した結果が表示されます。パラメータが先に手計算で 求めた最適なパラメータと近い値に収束し、ボール投げの記録がだいたい予測 できるようになっていることがわかります。