2012年 12月 1日掲載
「僕と契約して，AWK使いになってよ」
今回は，最小二乗法によるデータの平滑化をあつかってみる．
この記事は，昔何かのテキストで読んだ手法を思い出しながら書いたものなのだが，ざっとネットで調べたところでは，データの微分値を求める Savitzky-Golay法という手法を，ずっと簡便にしたものになるようだ．
まずは，サンプルデータの用意であるが，これは前回作った data1.txt を使う． 一応作り方を再掲しておく．
$ awk 'BEGIN{for(x=0;x<=6.0;x+=0.25)print 2*sin(x)}' > data0.txt
$ awk '{print $1+rand()-0.5}' data0.txt > data1.txt
data1.txt を，これまた前回作った plot1pbm.awk で図化すると，こんな感じ．
$ awk -f plot1pbm.awk data1.txt > data1a.pbm
さて，このギクシャクしたデータ列を，最小二乗法を用いて平滑化する方法を考えてみる．
まずは，隣接する 3点のデータを使って，1次近似式を作る方法を考えよう．
隣接する 3点を (xj-1, yj-1), (xj, yj), (xj+1, yj+1) とする． ただし，xj (j = 1, 2, ..., N) は，等間隔で，xj+1 - xj = Δx (j = 1, 2, ... , N-1) とする．
1次近似式を y~ = fj(x) = a (x - xj) + b とすると，最小二乗法を使って，次のようにして a, b を計算できる．
ここで，
であるから，
となる．
結果，近似式 fj(x) は，
となり，この式で，xj の近似値 y~j を計算すると，
となる． j = 2, 3, ... , N-1 の y~j は，3点の単純移動平均になるのだ．
fj(x) を使えば，単純移動平均では計算できなかった，データ列両端の y~1, y~N も計算できる． すなわち，
これを，AWK で書いたスクリプト smoothing1.awk は以下のようになる．
NR<3{
y[NR] = $1
}
NR==3{
y[NR] = $1
# y~_1
print (5 * y[1] + 2 * y[2] - y[3]) / 6
}
NR>3{
# y~_j, j = 2, ... , N-2
print (y[1] + y[2] + y[3]) / 3
y[1] = y[2]
y[2] = y[3]
y[3] = $1
}
END{
if(NR < 3) exit
# y~_N-1
print (y[1] + y[2] + y[3]) / 3
# y~_N
print (- y[1] + 2 * y[2] + 5 * y[3]) / 6
}
これを使って，data1.txt を平滑化してみる．
$ awk -f smoothing1.awk data1.txt > data3.txt $ awk -f plot1pbm.awk data3.txt > data3.pbm
データ数が減っていないこともたしかめておこうか．
$ wc -l data1.txt data3.txt
25 data1.txt
25 data3.txt
50 total
この手法は，容易に高次の近似式を使った方法に拡張できる． たとえば，隣接する 5点を使って，2次式で近似する方法を考えてみる．
近似式を，
として， x が等間隔の 5点のデータ (xj-2, yj-2), (xj-1, yj-1), (xj, yj), (xj+1, yj+1), (xj+2, yj+2) を使って，係数 a, b, c を最小二乗法で解くと，
となるので，
となる． 数字が多くなるので，まちがえないように計算するのが大変だが，3点 1次式のときと同じように計算式が求まる． これを AWK で書いた smoothing2.awk は，以下のとおり．
NR<5{
y[NR] = $1
}
NR==5{
y[NR] = $1
# y~_1
print (31 * y[1] + 9 * y[2] - 3 * y[3] - 5 * y[4] + 3 * y[5]) / 35
# y~_2
print (9 * y[1] + 13 * y[2] + 12 * y[3] + 6 * y[4] - 5 * y[5]) / 35
}
NR>5{
# y~_j, j = 3, ... , N-3
print (- 3 * y[1] + 12 * y[2] + 17 * y[3] + 12 * y[4] - 3 * y[5]) / 35
for(i = 1; i < 5; i++){
y[i] = y[i + 1]
}
y[5] = $1
}
END{
if(NR < 5) exit
# y~_N-2
print (- 3 * y[1] + 12 * y[2] + 17 * y[3] + 12 * y[4] - 3 * y[5]) / 35
# y~_N-1
print (- 5 * y[1] + 6 * y[2] + 12 * y[3] + 13 * y[4] + 9 * y[5]) / 35
# y~_N
print (3 * y[1] - 5 * y[2] - 3 * y[3] + 9 * y[4] + 31 * y[5]) / 35
}
計算が大変だったわりには，あっさりしたスクリプトである． data1.txt を，このスクリプトで処理すると，こうなる．
$ awk -f smoothing2.awk data1.txt > data4.txt $ awk -f plot1pbm.awk data4.txt > data4.pbm
え？移動平均（smoothing1.awk）の結果 data3.txt と，あんまり変わってないって？ たしかにグラフを見た限りでは，あまり変わったようには見えないが，最大値・最小値を見てみると（その1 で作った max.awk と min.awk を使う），
$ awk -f max.awk data3.txt 2.04519 $ awk -f max.awk data4.txt 2.16425 $ awk -f min.awk data3.txt -2.26189 $ awk -f min.awk data4.txt -2.33594
data4.txt の方が，data3.txt よりも，最大値が大きく，最小値が小さいことがわかる． smoothing2.awk の方が，2次式で近似している分，でっぱった部分のなまり方が，1次式で近似している smoothing1.awk よりも少ないようだ．
今回の AWKスクリプト：
過去の記事：
[ KDO ホームにもどる]
(C) 2012 KDO