元のページ

ステージ2-2: 線形回帰、グラフ描画の練習 (linear_model.LinearRegression, pl.scatter, pl.plot) (情報工学実験 3 : データマイニング班)

目次

想定環境

データセットの用意

# 1次元サンプルを -5〜+5 の範囲でn分割して用意
#ここでは、分割数 = 事例数。
import numpy as np
n = 50
x_base = np.linspace(-5, 5, n)

# 1次式(y = a*x +b) で出力を用意。
# 今回は y = 0.8*x - 4 とする。
a = 0.8
b = -4
y_base = np.polyval([a, b], x_base)

# 生成したデータセット(x_base, y_base) を散布図出力して確認
# pl.scatter(): s=点のサイズ
import pylab as pl
pl.figure()
pl.scatter(x_base, y_base, s=2)
pl.show()

# y1_base にノイズを付与したデータセットを用意
# 　0.5倍してるのはガウシアン分布のままだと、
# 　元の値(y_base)に対してノイズの割合が大きすぎるため。
y_noise = y_base + np.random.normal(size=n) * 0.5
pl.figure()
pl.scatter(x_base, y_noise, s=2)
pl.show()

# y_noise と実モデル(y=0.8*x-4)を同時に描画
# この実モデルを y_base や y_noise から推定したい。
pl.figure()
pl.plot(x_base, y_base)
pl.scatter(x_base, y_noise, s=2)
pl.show()






最小二乗法

# 上記リンクを参考に、
# ノイズがない場合 (x_base, y_base) と、
# ノイズがある場合 (x_base, y_noise) における a,b を算出せよ。






線形回帰で推定、グラフ描画の練習 (linear_model.LinearRegression)

# x_base は1サンプル=1数値として列挙されている。
# 1サンプルを1リストとして表現するため、型変換。
# さらに、このままでは一部の機能が使えないため、np.ndarray() に変換。
# オブジェクトの方は type() で確認可能。
data_list = [[x_base[i]] for i in range(len(x_base))]
data = np.array(data_list)
target = y_base

# 線形回帰モデルを準備
from sklearn import linear_model
regr = linear_model.LinearRegression()

# 手動で fitting してみる
#   LinearRegression(): coef_に傾き、intercept_に切片が得られる。
np.random.seed(0)
indices = np.random.permutation(len(data))
half = len(data)/2
data_train = data[indices[:half]]
data_test = data[indices[half:]]
target_train = target[indices[:half]]
target_test = target[indices[half:]]

regr.fit(data_train, target_train)
regr.score(data_test, target_test)
print "(LinearRegression:base) coef_=%s, intercept=%s" % (regr.coef_, regr.intercept_)

# ノイズデータでも fitting してみる
target = y_noise
target_train = target[indices[:half]]
target_test = target[indices[half:]]
regr.fit(data_train, target_train)
regr.score(data_test, target_test)
print "(LinearRegression:noise) coef_=%s, intercept=%s" % (regr.coef_, regr.intercept_)


# ノイズデータに対してシンプルに交差検定してみる
from sklearn import cross_validation
np.random.seed(0)
kfold = cross_validation.KFold(n=len(data), n_folds=2, shuffle=True)
scores = cross_validation.cross_val_score(regr, data, target, cv=kfold, n_jobs=-1)
print "scores =", scores
print "max = %s, ave = %s, min = %s" % (scores.max(), scores.mean(), scores.min())

# 交差検定しながら fitting 状況を描画
pl.figure()
pl.scatter(data, target, s=2)
for train_indices, test_indices in kfold:
    train_data = [data[i] for i in train_indices]
    train_target = [target[i] for i in train_indices]
    regr.fit(train_data, train_target)
    test_indices.sort()
    test_data = [data[i] for i in test_indices]
    test_target = [target[i] for i in test_indices]
    estimated = regr.predict(test_data)
    pl.plot(test_data, estimated)

pl.show()



参考サイト一覧