最小自乗法

回帰分析における最小自乗法の原理、回帰係数の導出過

程、C言語による実装

【学習内容の構成】

1. 回帰分析：説明変数と目的変数の関係をモデル化する

統計的手法

2. 最小自乗法：残差の自乗和を最小化して回帰係数a, bを

決定する手法

3. 数式導出：偏微分による連立方程式の解法で回帰係数

を算出

4. プログラム実装：C言語によるファイル入力と回帰係数

計算

• 前提：微分の基礎、C言語のファイル操作

• 意義：データから予測モデルを構築する能力の習得

回帰分析

• 説明変数と目的変数の関係のモデルを推

定する統計的手法

– ある変量(説明変数)とその変量に対する望み

の結果(目的変数)の値がいくつか与えられる

• 説明変数が2つ以上ある場合には，重回

帰分析と言う

回帰分析

• 説明変数( x )と目的変数( y )の関係の実測値がいく

つか与えられたとする（下図）

• x に対する y の予測値 y は，y = ax + b で与えられる

• この直線の最適な係数 a, b を求めることによって，変

数 x と y の関係が予測する

• 係数 a, b を回帰係数と呼ぶ

x

y

実測値

y=ax+b

^

最小自乗法

• 回帰係数を求める手法

• 実測値 y

i

と予測値 y

i

(=ax

i

+b) の差(残差)の自乗

和を最小にするような回帰係数 a,b を求める

– 実測値のデータの個数：N

– x

i

： i 番目の x の実測値 (i = 1・・・N)

– y

i

： x

i

に対する y の実測値とする (i = 1・・・N)

– a：直線の傾き

– b：切片

＾

最小自乗法

1月

2月

3月

4月

5月

110

150

120

130

？

5月の売上は？

月売上

実測値

月X

売上Y

Y＝3X＋120

最小自乗法

X

Y

y=a+bx+誤差

誤差の自乗和（図の正方形の和）を最小

にするようにして未知係数aとbを求める

ので最小自乗法と呼ぶ

最小自乗法

• 残差の自乗和: Q

Q = Σ(y

i

- y

i

)

2

= Σ(y

i

- (ax

i

+b))

2

(0)

i=1

N

i=1

N

＾

最小自乗法

• Qを最小とする a, b は， = 0 かつ = 0

を満たす

= -2Σx

i

(y

i

- (ax

i

+b)) = 0 - (1)

= -2Σ(y

i

- (ax

i

+b)) = 0 - (2)

(1)式を整理すると

Σx

i

(y

i

- (ax

i

+ b)) = 0

Σx

i

y

i

- Σ(ax

i

x

i

+ bx

i

)) = 0

Σx

i

y

i

- aΣx

i

x

i

= bΣx

i

- (3)

ρQ

ρa

ρQ

ρb

i=1

N

ρQ

ρa

ρQ

ρb

i=1

N

i=1

N

i=1

N

i=1

N

i=1

N

i=1

N

i=1

N

最小自乗法

• 前ページ(2)式を整理すると

Σy

i

–Σax

i

–Σb = 0

Σy

i

– aΣx

i

– Nb = 0

Σy

i

– aΣx

i

= Nb

b = 1/N (Σy

i

– aΣx

i

) - (4)

i=1

N

i=1

N

i=1

N

i=1

N

i=1

N

i=1

N

i=1

N

i=1

N

i=1

N

(4)式を(3)式に代入して

以上のようにして，回帰係数を決定

この時のa,bが予測値と実測値の差の自乗を最小にするような，

予測値の式の係数

Σxiyi-aΣxixi= Σxi

i=1

N

N

i=1

i=1

N

(Σyi-aΣxi)

i=1 i=1

N N

N

1

NΣxiyi-aNΣxixi = Σxi

i=1

N N

i=1 i=1

N

ΣxiΣyi-aΣxi

i=1 i=1

N N

i=1

N

a(ΣxiΣxi-NΣxixi) = ΣxiΣyi-NΣxiyi

i=1

N

i=1

N

i=1

N

i=1

N

i=1

N

i=1

N

ΣxiΣxi-NΣxixi

i=1

N

i=1

N

i=1

N

ΣxiΣyi-NΣxiyi

i=1

N

i=1

N

i=1

N

a =

(5)

#include<stdio.h>

main()

{

FILE *infile;

int i, n;

double a, b, x, y, x_sum, y_sum, xy_sum, xx_sum;

if ((infile=fopen("input.dat","r"))==NULL) {

printf("can't open file ¥n");

exit;

}

n=0;

xy_sum=0;

x_sum=0;

y_sum=0;

xx_sum=0;

while( fscanf(infile, "%lf%lf¥n", &x, &y) != EOF ) {

xy_sum=xy_sum+x*y; /* x*y の和 */

x_sum=x_sum+x; /* x の和 */

y_sum=y_sum+y; /* y の和 */

xx_sum=xx_sum+x*x; /* x*x の和 */

n=n+1; /* データの数 n */

}

a=(x_sum*y_sum-n*xy_sum) / (x_sum*x_sum-n*xx_sum);

b=(y_sum-a*x_sum)/n;

printf("y = ax + b¥n");

printf("a = %lf¥nb = %lf¥n", a, b);

fclose(infile);

}