情報工学実験I データ解析・直線によるデータの当てはめ（宮里）

情報工学実験 I
データ解析・直線によるデータの当てはめ（宮里）
学籍番号：065762A
氏名：与儀涼子
作成：19.07.13 提出：19.07.13
1
補足プログラムの動作
1
３つの補足プログラムの動作を確認せよ。
まず、プログラムの未完成の部分をコメントアウトしてコンパイルし、動作を確認した。
実行結果
----------------------------% ./lsf-demo
Curve fit data is created using the quadratic
y(x) = c(1) + c(2)*x + c(3)*x^2
Enter the coefficients:
c(1) = 1
c(2) = 2
c(3) = 3
Enter estimated error bar: 100
Fit parameters:
a(1) = 0 +/- 0
a(2) = 0 +/- 0
Chi square = 0; N-M = 48
----------------------------コンパイルに成功し、プログラムを実行することもできた。しかしまだ部分的に未完成であるた
め、データを入力しても結果として 0 が出力される。
線形回帰，最小二乗法について調べて簡潔にまとめよ。
2
2.1
線形回帰について
線形回帰とは、与えられた二次元のデータを一次関数で近似することである。全てのパラメータか
ら最も近い直線を選ぶための手法の代表として、以下の最小二乗法がある。
2.2
最小二乗法について
データを近似する手段の一つ。二次元データを一次関数で表す場合、推定値の直線と全てのデータ
との差分の二乗平均を求め、その値が最小になる直線が最も近似な直線である。
3
カイ二乗関数について調べて簡潔にまとめよ
データを最小二乗法を使って近似する場合、誤差範囲の大きい点ほど加重が減るような基準を用
いるとよい。つまり、推定値の直線から近い点が 1 ずれることは大きな差だが、推定値の直線から十
分遠い点が 1 ずれることの影響は小さい。
2
この基準を実現するために、推定値と実際のパラメータの差の二乗を、距離から求めた値で割ること
にする。この当てはめに使われるのがカイ二乗関数で、以下のように定義される。
χ ({aj }) =
2
)2
N (
∑
∆i
i=1
σi
=
N
∑
[Y (xi ; {aj }) − yi ]2
i=1
σi2
線形回帰当てはめプログラム
4
線形回帰当てはめプログラム (lin-reg.cpp) を完成させよ．
以下に lin reg.cpp のソースコードを示す。
// 線形回帰プログラム
// Inputs
//
x
独立変数
従属変数
//
y
//
sigma
// Outputs
y における推定誤差
//
//
a_fit
sig_a
a(1):切片, a(2):傾き
パラメータ a() における推定誤差
//
//
yy
chisqr
推定値
カイ二乗関数
#include "NumMeth.h"
#include "math.h"
void lin_reg(Matrix x, Matrix y, Matrix sigma,
Matrix &a_fit, Matrix &sig_a, Matrix &yy, double &chisqr) {
//* Evaluate various sigma sums
int i, nData = x.nRow();
double sigmaTerm;
double s = 0.0, sx = 0.0, sy = 0.0, sxy = 0.0, sxx = 0.0;
for( i=1; i<=nData; i++ ) {
//
式 (5.10) の計算を行う
sigmaTerm = 1/(sigma(i)*sigma(i));
s += sigmaTerm;
sx += x(i)*sigmaTerm;
sy += y(i)*sigmaTerm;
sxy += (x(i)*y(i))*sigmaTerm;
sxx += (x(i)*x(i))*sigmaTerm;
}
3
//* 切片 a_fit(1) と傾き a_fit(2) の計算
// ここに式 (5.11) に対応する部分を書く
//a_fit(1) = ・
・
・
・
a_fit(1) = ((sy*sxx) - (sx*sxy)) / ((s*sxx) - (sx*sx));
a_fit(2) = ((s*sxy) - (sy*sx)) / ((s*sxx) - (sx*sx));
//* 切片と傾きの誤差範囲を計算
//sig_a(1) = ・
・
・
・
sig_a(1) = sqrt(sxx / (s*sxx-(sx*sx)));
sig_a(2) = sqrt(s / (s*sxx-(sx*sx)));
//* 推定値とカイ二乗関数の計算
chisqr = 0.0;
//
ここに書く
for(i=1;i<=nData;i++){
chisqr += sigmaTerm * (a_fit(1) + (a_fit(2) * x(i)) - y(i))\\
* (a_fit(1) + (a_fit(2) * x(i)) - y(i));
yy(i) = a_fit(1) + (a_fit(2) * x(i));
}
}
4
5
octave を用いた結果のグラフ
octave を用いて結果をプロットし考察せよ．なお 3 種類のデータを発生させ，それぞれグラフ化
する事．
----------------------------% ./lsf-demo
Curve fit data is created using the quadratic
y(x) = c(1) + c(2)*x + c(3)*x^2
Enter the coefficients:
c(1) = 1
c(2) = 1
c(3) = 1
Enter estimated error bar: 1
Fit parameters:
a(1) = -441.301 +/- 0.287139
a(2) = 52.0032 +/- 0.00979992
Chi square = 1.73594e+06; N-M = 48
-----------------------------
緑色の直線が yy.txt、赤色の点で示されているのが y.txt である。赤色の点は、結ぶと二次曲線に
なるような分布である。
5
----------------------------% ./lsf-demo
Curve fit data is created using the quadratic
y(x) = c(1) + c(2)*x + c(3)*x^2
Enter the coefficients:
c(1) = 1
c(2) = 2
c(3) = 3
Enter estimated error bar: 300
Fit parameters:
a(1) = -1415.21 +/- 86.1418
a(2) = 155.948 +/- 2.93998
Chi square = 253.124; N-M = 48
-----------------------------
error bar(誤差の範囲) を大きくしてみた。今度は赤色の点の分布が二次曲線からかなりばらつい
たのがわかる。
6
----------------------------% ./lsf-demo
Curve fit data is created using the quadratic
y(x) = c(1) + c(2)*x + c(3)*x^2
Enter the coefficients:
c(1) = 1
c(2) = 2
c(3) = -3
Enter estimated error bar: 100
Fit parameters:
a(1) = 1296.93 +/- 28.7139
a(2) = -150.684 +/- 0.979992
Chi square = 1509.24; N-M = 48
-----------------------------
c(3) の値をマイナスにしてみた。上に凸の二次関数の点がばらけたような分布になる。
参考文献
[1] wikipedia http://ja.wikipedia.org/
[2] 授業配布のプリント
7

Download Report