10. 一次の最小二乗法

10. 一次の最小二乗法
2014/12/14
1
最小二乗法
実験やシミュレーションによって、データが得られた際に、そのデータが従う
式を推定することが多くあります。最小二乗法は、そのような方式の代表的なも
のです。
二次元のデータ点の列 {xk , yk }(k = 0, · · · , n − 1) があるとし、その関数形を
y = f (x) と予想したとします。関数には、一般にパラメタが含まれていることに
注意します。そのときに、データと予想した関数との二乗誤差
S=
n−1
∑
[yk − f (xk )]2
(1.1)
k=0
を最小化することで、データに合う曲線 y = f (x) を求めるのが最小二乗法です。
その最も簡単な形式が、予想した関数として一次関数 y = ax + b を仮定するも
のです。このときの二乗誤差は
S=
n−1
∑
[yk − axk − b]2
(1.2)
k=0
となります。S を最小化するような係数 a と b を求めることを考えましょう。
S は係数 a と b の下に凸な二次関数ですから、a と b を変化させた時の極値、つ
まり傾きゼロの点は、極小値となります。従って
n−1
∑
∂S
= −2
xk [yk − axk − b] = 0
∂a
k=0
(1.3)
n−1
∑
∂S
= −2
[yk − axk − b] = 0
∂b
k=0
(1.4)
の解、a と b が二乗誤差を最小化します。この式を整理すると
a
n−1
∑
x2k
k=0
n−1
∑
a
k=0
+b
n−1
∑
xk =
k=0
n−1
∑
xk + b
k=0
1
1 =
n−1
∑
k=0
n−1
∑
k=0
xk yk
(1.5)
yk
(1.6)
という連立方程式となります。ここで、
∑
1 n−1
xk
n k=0
(1.7)
=
∑
1 n−1
x2
n k=0 k
(1.8)
hyi =
∑
1 n−1
yk
n k=0
(1.9)
hxi =
D
x2
E
hxyi =
∑
1 n−1
xk yk
n k=0
(1.10)
を用いると、
D
E
a x2 + b hxi = hxyi
a hxi + b = hyi
(1.11)
(1.12)
と整理できます。従って、解として
1
[hxyi − hxi hyi]
σx2
]
1 [D 2 E
b =
x
hyi
−
hxi
hxyi
σx2
a =
D
σx2 =
E
x2 − hxi2
(1.13)
(1.14)
(1.15)
を得ることができます。
片対数と両対数
2
9
8
7
5
y
y
6
4
3
2
1
0
0
2
4
6
8
10
1010
100
10-10
10-20
10-30
10-40
10-50
10-60
10-70
10-80
10-90
0
x
10 20 30 40 50 60 70 80 90 100
x
図 1: データの例:1
データの従う関数として一次関数を採用した場合の最小二乗法は、簡単に以下
の二つの場合に拡張することができます。
2
図 1 左のようなデータがあったとしましょう。これは、当然直線に従うもので
はありません。その大きな特徴は、x が大きくなると、y が急速に小さくなること
ですそこで、y 軸を対数で表示してみましょう。それが図 1 右です。直線的に変化
している様子が見えます。y 軸の対数をとったものが、x の一次関数に従うので、
ln y = −ax + b
(2.1)
のような振る舞いです。ここで、ln y = loge y は自然対数を表していることに注意
します。式 (2.1) の両辺の指数関数をとると
y = e−ax+b = eb e−ax
(2.2)
101
10
9
8
7
6
5
4
3
2
1
0
100
10-1
10-2
y
y
つまり、y は x の指数関数に従うことが分かります。
10-3
10-4
10
-5
10
-6
10-7
0
2
4
6
8
10-8
100
10
101
x
102
x
103
104
図 2: データの例:2
次に図 2 左の図を見ましょう。これも、当然直線に従うものではありません。そ
の大きな特徴は、x が大きくなると、y が急速に小さくなることですが、図 1 より
も、減少の仕方が緩く見えます。そこで、x と y の両方を対数をとって表示したも
のが図 2 右です。直線になりました。
ln y = −a ln x + b
(2.3)
という振る舞いです。式 (2.3) の両辺の指数関数をとると
y = eb x−a
(2.4)
となります。つまり、y は x のべき関数ということが分かります。
3
準備
それでは、プログラムを作成する準備をしましょう。プロジェクト LeastSquare
を作成します。さらに、以下の URL から Java ソースファイルを取得し、プロジェ
クト LeastSquare のソースファイルの下に置きます。また、ライブラリ MyLib も
設定しましょう。
3
http://http://aoba.cc.saga-u.ac.jp/lecture
/ModelingAndSimulation/javasrc/LeastSquare/src.zip
4
データの準備
プロジェクト LeastSquare のソースファイルには、二つのパッケージがありま
す。最初のパッケージ data は、最小二乗法で合わせることができるデータを生成
します。クラス DataGenerator を実行すると、三つのデータファイルを生成しま
す。直線で合わせることができるデータ out-linear.txt、指数関数で合わせるこ
とができるデータ out-exp.txt、そしてべき関数で合わせることができるデータ
out-pow.txt の三つです。まずは、gnuplot をもちいて、データを図示しましょう。
5
データの関数形の推計
もう一つのパッケージ leastSquare は、データに対して最小二乗法で関数形を
推計するクラスが入っています。まず、データをファイルから読む必要がありま
す。クラス FileRead は、そのようなメソッドもつライブラリです。
メソッド openReader は、ファイル名を指定して、BufferedReader クラスのイ
ンスタンスを開きます。メソッド readData は、ファイル名を指定すると、内部で
メソッド openReader を使って BufferedReader クラスのインスタンスを開き、そ
こから、double 型数値がスペース区切りで一行に二つあると仮定して、データを
読み出します。各行のデータは、xy-座標のクラス Point2D.Double のインスタン
スとしてリストになって戻ります。
クラス LeastSquare は、ファイル名を指定してデータを読み出し、直線を使っ
た最小二乗法を実行します。このクラスは、前述の指数関数、べき関数にも対応
します。その区別は
static public enum Type {
Normal, SemiLog, LogLog
}
という enum 型で定義され、メソッド fit に引数として指定します。
最小二乗法の本体はプライベートメソッド fitSub です (Program 5.1)。引数は xy座標となったデータ点のリストです。指数関数やべき関数での推計が対象の場合に
は、すでに対数をとる処理がなされています。クラス Point2D.Double のインスタ
ンスを p とすると、p.x が x-座標の値、p.y が y-座標の値になります。Program 5.1
を完成させましょう。
クラス LeastSquare を実行すると、先に作った out-linear.txt、out-exp.txt、
及び out-pow.txt の三つを読み込み、それぞれの関数形を推計して画面に表示し
ます。データとともに、それらの関数形を図示しましょう。
4
private Point2D.Double fitSub(List<Point2D.Double> d) {
double ax = 0.;//x の平均
double ax2 = 0;//x^2 の平均
double ay = 0;//y の平均
double axy = 0.;//x*y の平均
int n = 0;//データ点の数
for (Point2D.Double p : d) {//全てのデータに対して
//x、x^2、y、xy の和を計算
}
ax /= n;
ax2 /= n;
ay /= n;
axy /= n;
double sigma2 = ax2 - ax * ax;
double a = (axy - ax * ay) / sigma2;
double b = (ay * ax2 - axy * ax) / sigma2;
return new Point2D.Double(a, b);
}
Program 5.1: fitSub メソッド
5

Download Report