sage/PRML-線形回帰
をテンプレートにして作成
[
トップ
] [
新規
|
一覧
|
単語検索
|
最終更新
|
ヘルプ
]
開始行:
[[FrontPage]]
#contents
2011/04/08からのアクセス回数 &counter;
ここで紹介したSageワークシートは、以下のURLからダウンロー...
http://sage.math.canterbury.ac.nz/home/pub/86/
また、Sageのサーバを公開しているサイト（http://sage.math....
アップロードし、実行したり、変更していろいろ動きを試すこ...
* 第1章-Sageを使って線形回帰を試してみる [#x22ef249]
[[パターン認識と機械学習>http://www.amazon.co.jp/dp/44311...
は、機械学習に関するとても優れた教科書です。ここでは、Sag...
第1章の例題として、\( sin(2 \pi x) \)の回帰問題について見...
** データについて [#bcfa72da]
上巻の付録Aにあるように http://research.microsoft.com/~cm...
例題のデータをダウンロードすることができます。
これから計算に使うSin曲線は、
$$
y = sin(2 \pi x) + \mathcal{N}(0,0.3)
$$
で与えられたデータです。ここでは、curve fitting dataで提...
最初に、データを座標Xと目的値tにセットし、sin曲線と一緒に...
sageへの入力:
#pre{{
# PRMLのsin曲線のデータ
data = matrix([
[0.000000, 0.349486],
[0.111111, 0.830839],
[0.222222, 1.007332],
[0.333333, 0.971507],
[0.444444, 0.133066],
[0.555556, 0.166823],
[0.666667, -0.848307],
[0.777778, -0.445686],
[0.888889, -0.563567],
[1.000000, 0.261502],
]);
X = data.column(0)
t = data.column(1)
M = 3;
# データのプロット
sin_plt = plot(sin(2*pi*x),[x, 0, 1], rgbcolor='green');
data_plt = list_plot(zip(X, t)); data_plt
(data_plt + sin_plt).show()
}}
&ref(sage0-0.png);
** 最小自乗法の解 [#q4d22cc3]
3章の式(3.15)によると最小自乗法の解は、
$$
W_{ML} = ( \Phi^T \Phi )^{-1} \Phi^T t
$$
で与えられます。ここで、\(\Phi\) は計画行列(design matrix...
与えらます。
また、行列 \( \Phi^{\dagger} \)は、ムーア＿ベンローズの疑...
$$
\Phi^{\dagger} = \left( \Phi^T \Phi \right)^{-1} \Phi^T
$$
*** 多項式フィッティング [#db96ba18]
多項式フィッティングの式(1.1)では、\(y(x, w)\)を以下のよ...
$$
y(x, w) = \sum^{M}_{j=0} w_j x^j
$$
そこで、\(\phi_j(x)\)を以下のように定義します。
$$
\phi_j(x) = x^j
$$
sageへの入力:
#pre{{
# Φ関数定義
def _phi(x, j):
return x^j
}}
*** 行列のセットと解の計算 [#r72c35af]
定義に従って計画行列\(\Phi\)、計画行列の転置行列\(\Phi^T\...
平均の重み\(W_{ML}\)を計算します。
sageへの入力:
#pre{{
# 計画行列Φ
Phi = matrix([[ _phi(x,j) for j in range(0, (M+1))] for x...
Phi_t = Phi.transpose();
# ムーア＿ベンローズの疑似逆行列
Phi_dag = (Phi_t * Phi).inverse() * Phi_t;
# 平均の重み
Wml = Phi_dag * t; Wml
}}
#pre{{
(0.313702727439780, 7.98537103199157, -25.4261022423404, ...
}}
*** 多項式yの定義 [#b562770f]
sageでは、多項式y(x)を以下のように定義します。
sageへの入力:
#pre{{
# 出力関数yの定義
y = lambda x : sum(Wml[i]*x^i for i in range(0, (M+1)));
}}
*** 多項式線形回帰の結果 [#d0d628cf]
\(M=3\)の時の、多項式回帰の結果（赤）をサンプリング（青）...
プロットします。
sageへの入力:
#pre{{
y_plt = plot(y, [x, 0, 1], rgbcolor='red');
(y_plt + data_plt + sin_plt).show();
}}
&ref(sage0-1.png);
*** 様々なMの結果 [#de724aca]
図1.3に相当する図を表示します。スライダーでMの値を変えて...
\(M=9\)ではすべての点を通る曲線になりますが、元のsin曲線...
これは、過学習（over fitting）と呼ばれる現象で、機械学習...
sageへの入力:
#pre{{
# PRMLの図1.3
@interact
def _(M=(0..9)):
# 計画行列Φ
Phi = matrix([[ _phi(x,j) for j in range(0, (M+1))] f...
Phi_t = Phi.transpose();
# ムーア＿ベンローズの疑似逆行列
Phi_dag = (Phi_t * Phi).inverse() * Phi_t;
# 平均の重み
Wml = Phi_dag * t;
f = lambda x : sum(Wml[i]*x^i for i in range(0, (M+1)...
y_plt = plot(f, [x, 0, 1], rgbcolor='red');
(y_plt + data_plt + sin_plt).show(ymin=-1.5, ymax=1.5);
}}
&ref(sage0-2.png);
*** 100個の検証用データを生成 [#p026eb42]
サンプルを100個に増やした検証用データを生成します。
sageへの入力:
#pre{{
# 100個の検証用データを生成する
X100 = vector([random() for i in range(100)]);
t100 = vector([(sin(2*pi*x) + +gauss(0, 0.3)).n() for x i...
# 100個のデータをプロット
lst100_plt = list_plot(zip(X100, t100), rgbcolor='gray');
(data_plt + lst100_plt + sin_plt).show()
}}
&ref(sage0-3.png);
*** サンプルを増やした場合 [#se6959a3]
サンプルを100個に増やし、\(M=9\)の多項式フィッティングを...
サンプルを増やせば、\(M=9\)でも元のsin曲線に近い形になり...
sageへの入力:
#pre{{
M = 9;
# 計画行列Φ
Phi = matrix([[ _phi(x,j) for j in range(0, (M+1))] for x...
Phi_t = Phi.transpose();
# ムーア＿ベンローズの疑似逆行列
Phi_dag = (Phi_t * Phi).inverse() * Phi_t;
# 平均の重み
Wml = Phi_dag * t100; Wml
f = lambda x : sum(Wml[i]*x^i for i in range(0, (M+1)));
y_plt = plot(f, [x, 0, 1], rgbcolor='red');
(y_plt + lst100_plt + sin_plt).show(ymin=-1.5, ymax=1.5);
}}
&ref(sage0-4.png);
*** 平均自乗平方根誤差 [#qa20890c]
原書の1.1では訓練データとテスト用データの平均自乗平方根誤...
$$
E(w) = \frac{1}{2} \sum^N_{n=1} \{ y(x_n, w) - t_n \}^2
$$
を使って最適なMの値を求める手法を紹介しています。
以下に訓練用データの\(E_{RMS}\)（青）、テスト用データの\(...
値の変化を示します。M=3で訓練、テスト共に低い値となること...
sageへの入力:
#pre{{
# 平均自乗平方根誤差
Erms_t = [];
Erms_t100 = [];
for M in range(10):
# 計画行列Φ
Phi = matrix([[ _phi(x,j) for j in range(0, (M+1))] f...
Phi_t = Phi.transpose();
# ムーア＿ベンローズの疑似逆行列
Phi_dag = (Phi_t * Phi).inverse() * Phi_t;
# 平均の重み
Wml = Phi_dag * t;
f = lambda x : sum(Wml[i]*x^i for i in range(0, (M+1)...
Erms_t += [sqrt(sum((t[i] - f(X[i]))^2 for i in range...
Erms_t100 += [sqrt(sum((t100[i] - f(X100[i]))^2 for i...
Erms_t_plt = list_plot(Erms_t);
Erms_t100_plt = list_plot(Erms_t100, rgbcolor = 'red');
(Erms_t_plt + Erms_t100_plt).show();
# グラフの傾向はPRMLと同じですが、値が若干ずれている？
}}
&ref(sage0-5.png);
** リッジ回帰 [#ze18b2d1]
誤差関数にペナルティ項を追加し、過学習を防ぐ例として、リ...
$$
\tilde{E}^(w) = \frac{1}{2} \sum^N_{n=1} \{ y(x_n, w) - t...
$$
を使って正規化します。
以下にM=9の練習用データにリッジ回帰を行った結果を示します。
sageへの入力:
#pre{{
# リッジ回帰
M = 9;
N = len(t);
lam = n(e^-18);
Phi = matrix([[ _phi(x,j) for j in range(0, (M+1))] for x...
Phi_t = Phi.transpose();
# ムーア＿ベンローズの疑似逆行列
Phi_dag = (lam*matrix((M+1),(M+1),1) + Phi_t * Phi).inver...
# 平均の重み
Wml = Phi_dag * t;
f = lambda x : sum(Wml[i]*x^i for i in range(0, (M+1)));
y_plt = plot(f, [x, 0, 1], rgbcolor='red');
(y_plt + data_plt + sin_plt).show(ymin=-1.5, ymax=1.5);
}}
&ref(sage0-6.png);
*** リッジ回帰に対する平均自乗平方根誤差 [#l010da42]
リッジ回帰に対する訓練用データの\(E_{RMS}\)（青）、テスト...
値の変化を示します。
sageへの入力:
#pre{{
# Ermsを取って過学習の度合いを見る
# 平均自乗平方根誤差
Erms_t = [];
Erms_t100 = [];
for ln_lam in range(-38, 1):
lam = n(e^ln_lam);
Phi = matrix([[ _phi(x,j) for j in range(0, (M+1))] f...
Phi_t = Phi.transpose();
# ムーア＿ベンローズの疑似逆行列
Phi_dag = (lam*matrix((M+1),(M+1),1) + Phi_t * Phi).i...
# 平均の重み
Wml = Phi_dag * t;
f = lambda x : sum(Wml[i]*x^i for i in range(0, (M+1)...
Erms_t += [[ln_lam, sqrt(sum((t[i] - f(X[i]))^2 for i...
Erms_t100 += [[ln_lam, sqrt(sum((t100[i] - f(X100[i])...
Erms_t_plt = list_plot(Erms_t);
Erms_t100_plt = list_plot(Erms_t100, rgbcolor = 'red');
(Erms_t_plt + Erms_t100_plt).show();
}}
&ref(sage0-7.png);
** ベイズ的なフィッティング [#n134ab6d]
もっとも良い結果を示すのが、ベイズ的なフィッティングの結...
$$
m_N = \beta S_N \Phi^T t
$$
$$
S^{-1}_{N} = \alpha I + \beta \Phi^T \Phi
$$
で与えられます。（\(m_N\)は、直線回帰の\(W_{ML}\)に相当し...
以下に\(\alpha = 5 \times 10^{-3}, \beta = 11.1\)の練習用...
sageへの入力:
#pre{{
# ベイズ的なフィッティング
# α=5*10^-3, β=11.1を使用
alpha = 5*10^-3;
beta = 11.1;
Phi = matrix([[ _phi(x,j) for j in range(0, (M+1))] for x...
Phi_t = Phi.transpose();
# ムーア＿ベンローズの疑似逆行列
Phi_dag = (alpha*matrix((M+1),(M+1),1) + beta*Phi_t * Phi...
# 平均の重み
Wml = beta*Phi_dag * t;
f = lambda x : sum(Wml[i]*x^i for i in range(0, (M+1)));
# 分散
def s(x):
phi_x = vector([x^i for i in range(M+1)]);
s_sqr = 1/beta + phi_x * Phi_dag * phi_x;
return sqrt(s_sqr);
s_u_plt = plot(lambda x : f(x) + s(x), [x, 0, 1], rgbcol...
s_d_plt = plot(lambda x : f(x) - s(x), [x, 0, 1], rgbcol...
y_plt = plot(f, [x, 0, 1], rgbcolor='red');
(y_plt + data_plt + sin_plt + s_u_plt + s_d_plt).show(ymi...
}}
&ref(sage0-8.png);
** lasso回帰 [#jcca8eb4]
鈴木由宇さんから、最近はlasso回帰も注目されていると聞き、
[[smlyさんのブログ「線形回帰モデルとか」>http://d.hatena....
から、cvx_optを使った方法を使わせて頂き、Sageでプロットし...
lasso回帰は、リッジ回帰のペナルティ項が２次から絶対値にな...
$$
\tilde{E}(w) = \frac{1}{2} \sum^N_{n=1} \{ y(x_n, w) - t_...
$$
sageへの入力:
#pre{{
# lasso回帰
# http://d.hatena.ne.jp/smly/20100630/1277904761
# を参考にlasso回帰を実施
# cvx_optの宣言
from cvxopt import solvers
from cvxopt.base import matrix as _matrix
import numpy as np
#
M = 9
lim = 30.0
# 計画行列Φ
phi = _matrix([[ float(_phi(x,j).n()) for j in range(0, (...
P = _matrix(float(0.0), (2*(M+1), 2*(M+1)))
P[:M+1, :M+1] = phi.T * phi
q = _matrix(float(0.0), (2*(M+1), 1))
t = _matrix(np.matrix(t)).T
q[:M+1] = -phi.T * t
Ident = _matrix(np.identity(M+1, float))
G = _matrix([[Ident, -Ident, _matrix(float(0.0), (1,M+1))...
h = _matrix(float(lim), (2*(M+1)+1,1))
# constraint (PRML ex.3.5, eq3.30)
x = solvers.qp(P, q, G, h)['x'][:M+1]
Wml = np.array(x).reshape(M+1)
print Wml
}}
sageからの出力:
#pre{{
pcost dcost gap pres dres
0: -7.3129e-01 -1.5730e+04 2e+04 5e-17 8e+01
1: -1.1908e+00 -2.3198e+02 2e+02 2e-16 1e+00
2: -1.6934e+00 -1.1863e+01 1e+01 2e-16 5e-02
3: -1.8817e+00 -2.4425e+00 6e-01 2e-16 2e-03
4: -1.8890e+00 -1.9346e+00 5e-02 2e-16 2e-04
5: -1.8949e+00 -1.9100e+00 2e-02 2e-16 4e-05
6: -1.8990e+00 -1.9027e+00 4e-03 3e-16 4e-14
7: -1.8998e+00 -1.9009e+00 1e-03 2e-16 3e-14
8: -1.9003e+00 -1.9005e+00 2e-04 3e-16 4e-14
9: -1.9004e+00 -1.9005e+00 1e-04 3e-16 3e-14
10: -1.9004e+00 -1.9004e+00 1e-05 5e-16 4e-14
11: -1.9004e+00 -1.9004e+00 2e-07 2e-16 3e-14
Optimal solution found.
[ 3.50668337e-01 4.72355549e+00 -5.66598952e-01 -3.3...
5.28640500e-04 5.46111228e+01 2.04434686e+01 -9.7...
-1.30982583e+02 8.49936038e+01]
}}
sageへの入力:
#pre{{
var('x')
f = lambda x : sum(Wml[i]*_phi(x,i) for i in range(0, (M+...
y_plt = plot(f(x), [x, 0, 1], rgbcolor='red');
(y_plt + data_plt + sin_plt).show(ymin=-1.5, ymax=1.5);
}}
&ref(sage0-9.png);
** コメント [#y99ad97a]
#vote(おもしろかった[7],そうでもない[0],わかりずらい[0])
皆様のご意見、ご希望をお待ちしております。
- lasso回帰を追加 -- [[竹本 浩]] &new{2011-05-31 (火) 23...
- Bishop の本の式からコードを書くのに困っていた。参考にさ...
- 姉妹書は、Netlab Algorithms for Pattern Recognition ISB...
- Netlab は Neural Networks for Pattern Recognition ,C.M...
- sage -ipython notebook で実行してみました。変換する時の...
- losso回帰 http://d.hatena.ne.jp/smly/20100630/127790476...
- https://github.com/PRML/PRML PRML toolbox を見つけまし...
- ysatoさま、貴重な情報をありがとうございます。 -- [[竹本...
- PRML toolbox は説明が無いものでした。それに対してこのペ...
- ysatoさま、うれしいコメントありがとうございます。 -- [[...
#comment_kcaptcha
終了行:
[[FrontPage]]
#contents
2011/04/08からのアクセス回数 &counter;
ここで紹介したSageワークシートは、以下のURLからダウンロー...
http://sage.math.canterbury.ac.nz/home/pub/86/
また、Sageのサーバを公開しているサイト（http://sage.math....
アップロードし、実行したり、変更していろいろ動きを試すこ...
* 第1章-Sageを使って線形回帰を試してみる [#x22ef249]
[[パターン認識と機械学習>http://www.amazon.co.jp/dp/44311...
は、機械学習に関するとても優れた教科書です。ここでは、Sag...
第1章の例題として、\( sin(2 \pi x) \)の回帰問題について見...
** データについて [#bcfa72da]
上巻の付録Aにあるように http://research.microsoft.com/~cm...
例題のデータをダウンロードすることができます。
これから計算に使うSin曲線は、
$$
y = sin(2 \pi x) + \mathcal{N}(0,0.3)
$$
で与えられたデータです。ここでは、curve fitting dataで提...
最初に、データを座標Xと目的値tにセットし、sin曲線と一緒に...
sageへの入力:
#pre{{
# PRMLのsin曲線のデータ
data = matrix([
[0.000000, 0.349486],
[0.111111, 0.830839],
[0.222222, 1.007332],
[0.333333, 0.971507],
[0.444444, 0.133066],
[0.555556, 0.166823],
[0.666667, -0.848307],
[0.777778, -0.445686],
[0.888889, -0.563567],
[1.000000, 0.261502],
]);
X = data.column(0)
t = data.column(1)
M = 3;
# データのプロット
sin_plt = plot(sin(2*pi*x),[x, 0, 1], rgbcolor='green');
data_plt = list_plot(zip(X, t)); data_plt
(data_plt + sin_plt).show()
}}
&ref(sage0-0.png);
** 最小自乗法の解 [#q4d22cc3]
3章の式(3.15)によると最小自乗法の解は、
$$
W_{ML} = ( \Phi^T \Phi )^{-1} \Phi^T t
$$
で与えられます。ここで、\(\Phi\) は計画行列(design matrix...
与えらます。
また、行列 \( \Phi^{\dagger} \)は、ムーア＿ベンローズの疑...
$$
\Phi^{\dagger} = \left( \Phi^T \Phi \right)^{-1} \Phi^T
$$
*** 多項式フィッティング [#db96ba18]
多項式フィッティングの式(1.1)では、\(y(x, w)\)を以下のよ...
$$
y(x, w) = \sum^{M}_{j=0} w_j x^j
$$
そこで、\(\phi_j(x)\)を以下のように定義します。
$$
\phi_j(x) = x^j
$$
sageへの入力:
#pre{{
# Φ関数定義
def _phi(x, j):
return x^j
}}
*** 行列のセットと解の計算 [#r72c35af]
定義に従って計画行列\(\Phi\)、計画行列の転置行列\(\Phi^T\...
平均の重み\(W_{ML}\)を計算します。
sageへの入力:
#pre{{
# 計画行列Φ
Phi = matrix([[ _phi(x,j) for j in range(0, (M+1))] for x...
Phi_t = Phi.transpose();
# ムーア＿ベンローズの疑似逆行列
Phi_dag = (Phi_t * Phi).inverse() * Phi_t;
# 平均の重み
Wml = Phi_dag * t; Wml
}}
#pre{{
(0.313702727439780, 7.98537103199157, -25.4261022423404, ...
}}
*** 多項式yの定義 [#b562770f]
sageでは、多項式y(x)を以下のように定義します。
sageへの入力:
#pre{{
# 出力関数yの定義
y = lambda x : sum(Wml[i]*x^i for i in range(0, (M+1)));
}}
*** 多項式線形回帰の結果 [#d0d628cf]
\(M=3\)の時の、多項式回帰の結果（赤）をサンプリング（青）...
プロットします。
sageへの入力:
#pre{{
y_plt = plot(y, [x, 0, 1], rgbcolor='red');
(y_plt + data_plt + sin_plt).show();
}}
&ref(sage0-1.png);
*** 様々なMの結果 [#de724aca]
図1.3に相当する図を表示します。スライダーでMの値を変えて...
\(M=9\)ではすべての点を通る曲線になりますが、元のsin曲線...
これは、過学習（over fitting）と呼ばれる現象で、機械学習...
sageへの入力:
#pre{{
# PRMLの図1.3
@interact
def _(M=(0..9)):
# 計画行列Φ
Phi = matrix([[ _phi(x,j) for j in range(0, (M+1))] f...
Phi_t = Phi.transpose();
# ムーア＿ベンローズの疑似逆行列
Phi_dag = (Phi_t * Phi).inverse() * Phi_t;
# 平均の重み
Wml = Phi_dag * t;
f = lambda x : sum(Wml[i]*x^i for i in range(0, (M+1)...
y_plt = plot(f, [x, 0, 1], rgbcolor='red');
(y_plt + data_plt + sin_plt).show(ymin=-1.5, ymax=1.5);
}}
&ref(sage0-2.png);
*** 100個の検証用データを生成 [#p026eb42]
サンプルを100個に増やした検証用データを生成します。
sageへの入力:
#pre{{
# 100個の検証用データを生成する
X100 = vector([random() for i in range(100)]);
t100 = vector([(sin(2*pi*x) + +gauss(0, 0.3)).n() for x i...
# 100個のデータをプロット
lst100_plt = list_plot(zip(X100, t100), rgbcolor='gray');
(data_plt + lst100_plt + sin_plt).show()
}}
&ref(sage0-3.png);
*** サンプルを増やした場合 [#se6959a3]
サンプルを100個に増やし、\(M=9\)の多項式フィッティングを...
サンプルを増やせば、\(M=9\)でも元のsin曲線に近い形になり...
sageへの入力:
#pre{{
M = 9;
# 計画行列Φ
Phi = matrix([[ _phi(x,j) for j in range(0, (M+1))] for x...
Phi_t = Phi.transpose();
# ムーア＿ベンローズの疑似逆行列
Phi_dag = (Phi_t * Phi).inverse() * Phi_t;
# 平均の重み
Wml = Phi_dag * t100; Wml
f = lambda x : sum(Wml[i]*x^i for i in range(0, (M+1)));
y_plt = plot(f, [x, 0, 1], rgbcolor='red');
(y_plt + lst100_plt + sin_plt).show(ymin=-1.5, ymax=1.5);
}}
&ref(sage0-4.png);
*** 平均自乗平方根誤差 [#qa20890c]
原書の1.1では訓練データとテスト用データの平均自乗平方根誤...
$$
E(w) = \frac{1}{2} \sum^N_{n=1} \{ y(x_n, w) - t_n \}^2
$$
を使って最適なMの値を求める手法を紹介しています。
以下に訓練用データの\(E_{RMS}\)（青）、テスト用データの\(...
値の変化を示します。M=3で訓練、テスト共に低い値となること...
sageへの入力:
#pre{{
# 平均自乗平方根誤差
Erms_t = [];
Erms_t100 = [];
for M in range(10):
# 計画行列Φ
Phi = matrix([[ _phi(x,j) for j in range(0, (M+1))] f...
Phi_t = Phi.transpose();
# ムーア＿ベンローズの疑似逆行列
Phi_dag = (Phi_t * Phi).inverse() * Phi_t;
# 平均の重み
Wml = Phi_dag * t;
f = lambda x : sum(Wml[i]*x^i for i in range(0, (M+1)...
Erms_t += [sqrt(sum((t[i] - f(X[i]))^2 for i in range...
Erms_t100 += [sqrt(sum((t100[i] - f(X100[i]))^2 for i...
Erms_t_plt = list_plot(Erms_t);
Erms_t100_plt = list_plot(Erms_t100, rgbcolor = 'red');
(Erms_t_plt + Erms_t100_plt).show();
# グラフの傾向はPRMLと同じですが、値が若干ずれている？
}}
&ref(sage0-5.png);
** リッジ回帰 [#ze18b2d1]
誤差関数にペナルティ項を追加し、過学習を防ぐ例として、リ...
$$
\tilde{E}^(w) = \frac{1}{2} \sum^N_{n=1} \{ y(x_n, w) - t...
$$
を使って正規化します。
以下にM=9の練習用データにリッジ回帰を行った結果を示します。
sageへの入力:
#pre{{
# リッジ回帰
M = 9;
N = len(t);
lam = n(e^-18);
Phi = matrix([[ _phi(x,j) for j in range(0, (M+1))] for x...
Phi_t = Phi.transpose();
# ムーア＿ベンローズの疑似逆行列
Phi_dag = (lam*matrix((M+1),(M+1),1) + Phi_t * Phi).inver...
# 平均の重み
Wml = Phi_dag * t;
f = lambda x : sum(Wml[i]*x^i for i in range(0, (M+1)));
y_plt = plot(f, [x, 0, 1], rgbcolor='red');
(y_plt + data_plt + sin_plt).show(ymin=-1.5, ymax=1.5);
}}
&ref(sage0-6.png);
*** リッジ回帰に対する平均自乗平方根誤差 [#l010da42]
リッジ回帰に対する訓練用データの\(E_{RMS}\)（青）、テスト...
値の変化を示します。
sageへの入力:
#pre{{
# Ermsを取って過学習の度合いを見る
# 平均自乗平方根誤差
Erms_t = [];
Erms_t100 = [];
for ln_lam in range(-38, 1):
lam = n(e^ln_lam);
Phi = matrix([[ _phi(x,j) for j in range(0, (M+1))] f...
Phi_t = Phi.transpose();
# ムーア＿ベンローズの疑似逆行列
Phi_dag = (lam*matrix((M+1),(M+1),1) + Phi_t * Phi).i...
# 平均の重み
Wml = Phi_dag * t;
f = lambda x : sum(Wml[i]*x^i for i in range(0, (M+1)...
Erms_t += [[ln_lam, sqrt(sum((t[i] - f(X[i]))^2 for i...
Erms_t100 += [[ln_lam, sqrt(sum((t100[i] - f(X100[i])...
Erms_t_plt = list_plot(Erms_t);
Erms_t100_plt = list_plot(Erms_t100, rgbcolor = 'red');
(Erms_t_plt + Erms_t100_plt).show();
}}
&ref(sage0-7.png);
** ベイズ的なフィッティング [#n134ab6d]
もっとも良い結果を示すのが、ベイズ的なフィッティングの結...
$$
m_N = \beta S_N \Phi^T t
$$
$$
S^{-1}_{N} = \alpha I + \beta \Phi^T \Phi
$$
で与えられます。（\(m_N\)は、直線回帰の\(W_{ML}\)に相当し...
以下に\(\alpha = 5 \times 10^{-3}, \beta = 11.1\)の練習用...
sageへの入力:
#pre{{
# ベイズ的なフィッティング
# α=5*10^-3, β=11.1を使用
alpha = 5*10^-3;
beta = 11.1;
Phi = matrix([[ _phi(x,j) for j in range(0, (M+1))] for x...
Phi_t = Phi.transpose();
# ムーア＿ベンローズの疑似逆行列
Phi_dag = (alpha*matrix((M+1),(M+1),1) + beta*Phi_t * Phi...
# 平均の重み
Wml = beta*Phi_dag * t;
f = lambda x : sum(Wml[i]*x^i for i in range(0, (M+1)));
# 分散
def s(x):
phi_x = vector([x^i for i in range(M+1)]);
s_sqr = 1/beta + phi_x * Phi_dag * phi_x;
return sqrt(s_sqr);
s_u_plt = plot(lambda x : f(x) + s(x), [x, 0, 1], rgbcol...
s_d_plt = plot(lambda x : f(x) - s(x), [x, 0, 1], rgbcol...
y_plt = plot(f, [x, 0, 1], rgbcolor='red');
(y_plt + data_plt + sin_plt + s_u_plt + s_d_plt).show(ymi...
}}
&ref(sage0-8.png);
** lasso回帰 [#jcca8eb4]
鈴木由宇さんから、最近はlasso回帰も注目されていると聞き、
[[smlyさんのブログ「線形回帰モデルとか」>http://d.hatena....
から、cvx_optを使った方法を使わせて頂き、Sageでプロットし...
lasso回帰は、リッジ回帰のペナルティ項が２次から絶対値にな...
$$
\tilde{E}(w) = \frac{1}{2} \sum^N_{n=1} \{ y(x_n, w) - t_...
$$
sageへの入力:
#pre{{
# lasso回帰
# http://d.hatena.ne.jp/smly/20100630/1277904761
# を参考にlasso回帰を実施
# cvx_optの宣言
from cvxopt import solvers
from cvxopt.base import matrix as _matrix
import numpy as np
#
M = 9
lim = 30.0
# 計画行列Φ
phi = _matrix([[ float(_phi(x,j).n()) for j in range(0, (...
P = _matrix(float(0.0), (2*(M+1), 2*(M+1)))
P[:M+1, :M+1] = phi.T * phi
q = _matrix(float(0.0), (2*(M+1), 1))
t = _matrix(np.matrix(t)).T
q[:M+1] = -phi.T * t
Ident = _matrix(np.identity(M+1, float))
G = _matrix([[Ident, -Ident, _matrix(float(0.0), (1,M+1))...
h = _matrix(float(lim), (2*(M+1)+1,1))
# constraint (PRML ex.3.5, eq3.30)
x = solvers.qp(P, q, G, h)['x'][:M+1]
Wml = np.array(x).reshape(M+1)
print Wml
}}
sageからの出力:
#pre{{
pcost dcost gap pres dres
0: -7.3129e-01 -1.5730e+04 2e+04 5e-17 8e+01
1: -1.1908e+00 -2.3198e+02 2e+02 2e-16 1e+00
2: -1.6934e+00 -1.1863e+01 1e+01 2e-16 5e-02
3: -1.8817e+00 -2.4425e+00 6e-01 2e-16 2e-03
4: -1.8890e+00 -1.9346e+00 5e-02 2e-16 2e-04
5: -1.8949e+00 -1.9100e+00 2e-02 2e-16 4e-05
6: -1.8990e+00 -1.9027e+00 4e-03 3e-16 4e-14
7: -1.8998e+00 -1.9009e+00 1e-03 2e-16 3e-14
8: -1.9003e+00 -1.9005e+00 2e-04 3e-16 4e-14
9: -1.9004e+00 -1.9005e+00 1e-04 3e-16 3e-14
10: -1.9004e+00 -1.9004e+00 1e-05 5e-16 4e-14
11: -1.9004e+00 -1.9004e+00 2e-07 2e-16 3e-14
Optimal solution found.
[ 3.50668337e-01 4.72355549e+00 -5.66598952e-01 -3.3...
5.28640500e-04 5.46111228e+01 2.04434686e+01 -9.7...
-1.30982583e+02 8.49936038e+01]
}}
sageへの入力:
#pre{{
var('x')
f = lambda x : sum(Wml[i]*_phi(x,i) for i in range(0, (M+...
y_plt = plot(f(x), [x, 0, 1], rgbcolor='red');
(y_plt + data_plt + sin_plt).show(ymin=-1.5, ymax=1.5);
}}
&ref(sage0-9.png);
** コメント [#y99ad97a]
#vote(おもしろかった[7],そうでもない[0],わかりずらい[0])
皆様のご意見、ご希望をお待ちしております。
- lasso回帰を追加 -- [[竹本 浩]] &new{2011-05-31 (火) 23...
- Bishop の本の式からコードを書くのに困っていた。参考にさ...
- 姉妹書は、Netlab Algorithms for Pattern Recognition ISB...
- Netlab は Neural Networks for Pattern Recognition ,C.M...
- sage -ipython notebook で実行してみました。変換する時の...
- losso回帰 http://d.hatena.ne.jp/smly/20100630/127790476...
- https://github.com/PRML/PRML PRML toolbox を見つけまし...
- ysatoさま、貴重な情報をありがとうございます。 -- [[竹本...
- PRML toolbox は説明が無いものでした。それに対してこのペ...
- ysatoさま、うれしいコメントありがとうございます。 -- [[...
#comment_kcaptcha
ページ名:
SmartDoc
![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
