sage/text/データフィッティング
をテンプレートにして作成
[
トップ
] [
新規
|
一覧
|
単語検索
|
最終更新
|
ヘルプ
]
開始行:
[[FrontPage]]
#contents
2011/06/19からのアクセス回数 &counter;
ここで紹介したSageワークシートは、以下のURLからダウンロー...
http://www15191ue.sakura.ne.jp:8000/home/pub/8/
また、Sageのサーバを公開しているサイト（http://www.sagenb...
アップロードし、実行したり、変更していろいろ動きを試すこ...
** データフィッティング [#u39725eb]
与えられたデータに最もよくフィットする関数を求めるのが、...
今回は、\( sin(2 \pi x) \)曲線に正規分布のノイズを加えた...
近似しながら、Sageでのデータフィッティングの方法と手法の...
*** テストデータの生成 [#mb4e469a]
テストデータは、区間[0, 1]でランダムに抽出したxに対して、...
$$
y = sin(2 \pi x) + \mathcal{N}(0,0.3)
$$
変数X, Yにそれぞれ１０個のx座標、y座標のリストをセットし...
このままで計算するたびにX, Yの値が異なってしまうので、一...
load関数で読み込んでいます。
sageへの入力:
#pre{{
# テストデータ生成
#X = [random() for i in range (10)]
#Y = [sin(2*pi*x) + gauss(0, 0.3) for x in X]
# データの保存
#save(X, DATA+'X')
#save(Y, DATA+'Y')
# データのリストア
X = load(DATA+'X')
Y = load(DATA+'Y')
}}
生成されたデータをもとの\(sin(2 \pi x) \)曲線と一緒にプロ...
後で、少ないデータでのフィッティングに使うため、最初の３...
残りを青でプロットしています。
sageへの入力:
#pre{{
# データのプロット
x = var('x')
sin_plt = plot(sin(2*pi*x),[x, 0, 1], rgbcolor='green')
blue_plt = list_plot(zip(X[3:], Y[3:]))
red_plt = list_plot(zip(X[:3], Y[:3]), rgbcolor='red')
data_plt = list_plot(zip(X, Y));
(blue_plt + red_plt + sin_plt).show(xmin=0, xmax= 1, ymin...
}}
&ref(out1.png);
** フィッティング関数 [#y9724116]
Sageが提供している曲線のフィッティング関数がfind_fit関数...
#pre{{
find_fit(データ, モデル, オプション)
}}
データには、モデルで使用する変数の値（ここではx）リストと...
リストを渡し、モデルには変数を与えると予測値を計算するモ...
今回の例では、modelとして、以下の３次多項式を使用します。
$$
model(x) = w_0 + w_1 x + w_2 x^2 + w_3 x^3
$$
find_fitで求められた\(w_0, w_1, w_2, w_3\)を使って関数f_f...
辞書型で結果をもらうと便利です。solution_dict=Trueオプシ...
辞書型で返されます。
sageへの入力:
#pre{{
(w0, w1, w2, w3) = var('w0 w1 w2 w3')
model(x) = w0 + w1*x + w2*x^2 + w3*x^3
data = zip(X, Y)
fit = find_fit(data, model, solution_dict=True); view(fit)
}}
*** フィッティング曲線のプロット [#e33bc607]
モデルmodelにfind_fitの結果fitを代入して、フィッティング...
赤で表された曲線がデータとうまく合っていることが見て取れ...
sageへの入力:
#pre{{
f_fit(x) = model.subs(fit)
fit_plt = plot(f_fit, [x, 0, 1], rgbcolor='red')
(blue_plt + red_plt + sin_plt + fit_plt).show(xmin=0, xma...
}}
&ref(out2.png);
*** データ数が足りない場合 [#j70fcf4e]
モデルmodelでは、４個の未知数$w_0, w_1, w_2, w_3$を求めて...
データの個数が４個よりも少ない場合には、どのようになるの...
１０個のデータのうち最初の３個（赤でプロットされた点）を...
以下のようにエラーとなってしまいます。
sageへの入力:
#pre{{
# データ数が足りない場合には、エラーとなる
data = zip(X[:3], Y[:3]; print len(data)
fit = find_fit(data, model, solution_dict=True); print fit
f_fit(x) = model.subs(fit)
fit_plt = plot(f_fit, [x, 0, 1])
(blue_plt + red_plt + sin_plt + fit_plt).show()
}}
sageの出力:
#pre{{
Traceback (most recent call last):
...
SyntaxError: invalid syntax
}}
** 最小二乗法 [#w0471463]
多項式によるフィッティング関数yを以下のように定義し、
$$
y(x, w) = \sum^{M}_{j=0} w_j x^j
$$
各項の\(x^n\)を以下のように表すと、
$$
\phi_j(x) = x^j
$$
観測値YとXの関係を行列で表すと以下の式でもっと誤差の少な...
$$
Y \approx \Phi w
$$
ここで、\(\Phi\)の各要素は、以下の式で表されます。
$$
\Phi_{nj} = \phi_j(x_n)
$$
[[これなら分かる最適化数学>href="http://www.amazon.co.jp/...
によると、
ムーア・ベンローズの一般逆行列\(\Phi^{\dagger}\)を使うと...
$$
w = \Phi^{\dagger} Y
$$
$\Phi^{\dagger}$は、データの個数n、未知数wの個数mとすると...
$$
\Phi^{\dagger} = \left\{\begin{eqnarray}
( \Phi^T \Phi )^{-1} \Phi^T & m > n \\
\Phi^T (\Phi \Phi^T)^{-1} & m < n
\end{eqnarray}\right.
$$
ムーア・ベンローズの一般逆行列\(\Phi^{\dagger}\)を求める...
フィッティング関数を求めることができます。
sageへの入力:
#pre{{
# 次数のセット
M = 3
# Φ関数定義
def _phi(x, j):
return x^j
}}
*** ムーア・ベンローズの一般逆行列 [#b6cb9533]
最初にデータ数が１０個の場合（n > m）の場合の、ムーア・ベ...
\(\Phi^{\dagger}\)をほとんど定義と同じような形でSageで書...
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 * vector(Y)
}}
*** 多項式yの定義 [#g0e152a1]
次に、求まった重みWmlを使って多項式y(x)を以下のように定義...
sageへの入力:
#pre{{
# 出力関数yの定義
y = lambda x : sum(Wml[i]*x^i for i in (0..M));
}}
*** 多項式線形回帰の結果 [#w33766d3]
データ数N=10, 未知数M=4の時の、多項式回帰の結果（赤）をサ...
プロットします。
sageへの入力:
#pre{{
y_plt = plot(y, [x, 0, 1], rgbcolor='red');
(y_plt + data_plt + sin_plt).show(xmin=0, xmax= 1, ymin=-...
}}
&ref(out3.png);
** データが少ない場合のムーア・ベンローズの一般逆行列 [#x...
それでは、データ数N=3, 未知数M=4の時のムーア・ベンローズ...
その結果をプロットしてみましょう。
sageへの入力:
#pre{{
# 計画行列Φ
Phi = matrix([[ _phi(x,j) for j in range(0, (M+1))] for x...
Phi_t = Phi.transpose()
# ムーア・ベンローズの一般逆行列
PhiPhit = Phi*Phi_t
Phi_dag = Phi_t*PhiPhit.inverse(); view(Phi_dag)
}}
&ref(out4.png);
sageへの入力:
#pre{{
# 平均の重み
Wml = Phi_dag * vector(Y[:3]); Wml.apply_map(lambda x : n...
}}
sageの出力:
#pre{{
(0.100463412348611, 7.01017348792056, -5.25466361108248, ...
}}
与えられた３点を通る3次曲線が求まっていることがわかります。
sageへの入力:
#pre{{
x = var('x')
y_plt = plot(y, [x, -2, 1.5], rgbcolor='blue')
sin_plt = plot(sin(2*pi*x),[x, -2, 1.5], rgbcolor='green')
(y_plt + red_plt + sin_plt).show()
}}
&ref(out5.png);
** コメント [#f2be9b8b]
#vote(おもしろかった[2],そうでもない[0],わかりずらい[0])
皆様のご意見、ご希望をお待ちしております。
#comment_kcaptcha
終了行:
[[FrontPage]]
#contents
2011/06/19からのアクセス回数 &counter;
ここで紹介したSageワークシートは、以下のURLからダウンロー...
http://www15191ue.sakura.ne.jp:8000/home/pub/8/
また、Sageのサーバを公開しているサイト（http://www.sagenb...
アップロードし、実行したり、変更していろいろ動きを試すこ...
** データフィッティング [#u39725eb]
与えられたデータに最もよくフィットする関数を求めるのが、...
今回は、\( sin(2 \pi x) \)曲線に正規分布のノイズを加えた...
近似しながら、Sageでのデータフィッティングの方法と手法の...
*** テストデータの生成 [#mb4e469a]
テストデータは、区間[0, 1]でランダムに抽出したxに対して、...
$$
y = sin(2 \pi x) + \mathcal{N}(0,0.3)
$$
変数X, Yにそれぞれ１０個のx座標、y座標のリストをセットし...
このままで計算するたびにX, Yの値が異なってしまうので、一...
load関数で読み込んでいます。
sageへの入力:
#pre{{
# テストデータ生成
#X = [random() for i in range (10)]
#Y = [sin(2*pi*x) + gauss(0, 0.3) for x in X]
# データの保存
#save(X, DATA+'X')
#save(Y, DATA+'Y')
# データのリストア
X = load(DATA+'X')
Y = load(DATA+'Y')
}}
生成されたデータをもとの\(sin(2 \pi x) \)曲線と一緒にプロ...
後で、少ないデータでのフィッティングに使うため、最初の３...
残りを青でプロットしています。
sageへの入力:
#pre{{
# データのプロット
x = var('x')
sin_plt = plot(sin(2*pi*x),[x, 0, 1], rgbcolor='green')
blue_plt = list_plot(zip(X[3:], Y[3:]))
red_plt = list_plot(zip(X[:3], Y[:3]), rgbcolor='red')
data_plt = list_plot(zip(X, Y));
(blue_plt + red_plt + sin_plt).show(xmin=0, xmax= 1, ymin...
}}
&ref(out1.png);
** フィッティング関数 [#y9724116]
Sageが提供している曲線のフィッティング関数がfind_fit関数...
#pre{{
find_fit(データ, モデル, オプション)
}}
データには、モデルで使用する変数の値（ここではx）リストと...
リストを渡し、モデルには変数を与えると予測値を計算するモ...
今回の例では、modelとして、以下の３次多項式を使用します。
$$
model(x) = w_0 + w_1 x + w_2 x^2 + w_3 x^3
$$
find_fitで求められた\(w_0, w_1, w_2, w_3\)を使って関数f_f...
辞書型で結果をもらうと便利です。solution_dict=Trueオプシ...
辞書型で返されます。
sageへの入力:
#pre{{
(w0, w1, w2, w3) = var('w0 w1 w2 w3')
model(x) = w0 + w1*x + w2*x^2 + w3*x^3
data = zip(X, Y)
fit = find_fit(data, model, solution_dict=True); view(fit)
}}
*** フィッティング曲線のプロット [#e33bc607]
モデルmodelにfind_fitの結果fitを代入して、フィッティング...
赤で表された曲線がデータとうまく合っていることが見て取れ...
sageへの入力:
#pre{{
f_fit(x) = model.subs(fit)
fit_plt = plot(f_fit, [x, 0, 1], rgbcolor='red')
(blue_plt + red_plt + sin_plt + fit_plt).show(xmin=0, xma...
}}
&ref(out2.png);
*** データ数が足りない場合 [#j70fcf4e]
モデルmodelでは、４個の未知数$w_0, w_1, w_2, w_3$を求めて...
データの個数が４個よりも少ない場合には、どのようになるの...
１０個のデータのうち最初の３個（赤でプロットされた点）を...
以下のようにエラーとなってしまいます。
sageへの入力:
#pre{{
# データ数が足りない場合には、エラーとなる
data = zip(X[:3], Y[:3]; print len(data)
fit = find_fit(data, model, solution_dict=True); print fit
f_fit(x) = model.subs(fit)
fit_plt = plot(f_fit, [x, 0, 1])
(blue_plt + red_plt + sin_plt + fit_plt).show()
}}
sageの出力:
#pre{{
Traceback (most recent call last):
...
SyntaxError: invalid syntax
}}
** 最小二乗法 [#w0471463]
多項式によるフィッティング関数yを以下のように定義し、
$$
y(x, w) = \sum^{M}_{j=0} w_j x^j
$$
各項の\(x^n\)を以下のように表すと、
$$
\phi_j(x) = x^j
$$
観測値YとXの関係を行列で表すと以下の式でもっと誤差の少な...
$$
Y \approx \Phi w
$$
ここで、\(\Phi\)の各要素は、以下の式で表されます。
$$
\Phi_{nj} = \phi_j(x_n)
$$
[[これなら分かる最適化数学>href="http://www.amazon.co.jp/...
によると、
ムーア・ベンローズの一般逆行列\(\Phi^{\dagger}\)を使うと...
$$
w = \Phi^{\dagger} Y
$$
$\Phi^{\dagger}$は、データの個数n、未知数wの個数mとすると...
$$
\Phi^{\dagger} = \left\{\begin{eqnarray}
( \Phi^T \Phi )^{-1} \Phi^T & m > n \\
\Phi^T (\Phi \Phi^T)^{-1} & m < n
\end{eqnarray}\right.
$$
ムーア・ベンローズの一般逆行列\(\Phi^{\dagger}\)を求める...
フィッティング関数を求めることができます。
sageへの入力:
#pre{{
# 次数のセット
M = 3
# Φ関数定義
def _phi(x, j):
return x^j
}}
*** ムーア・ベンローズの一般逆行列 [#b6cb9533]
最初にデータ数が１０個の場合（n > m）の場合の、ムーア・ベ...
\(\Phi^{\dagger}\)をほとんど定義と同じような形でSageで書...
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 * vector(Y)
}}
*** 多項式yの定義 [#g0e152a1]
次に、求まった重みWmlを使って多項式y(x)を以下のように定義...
sageへの入力:
#pre{{
# 出力関数yの定義
y = lambda x : sum(Wml[i]*x^i for i in (0..M));
}}
*** 多項式線形回帰の結果 [#w33766d3]
データ数N=10, 未知数M=4の時の、多項式回帰の結果（赤）をサ...
プロットします。
sageへの入力:
#pre{{
y_plt = plot(y, [x, 0, 1], rgbcolor='red');
(y_plt + data_plt + sin_plt).show(xmin=0, xmax= 1, ymin=-...
}}
&ref(out3.png);
** データが少ない場合のムーア・ベンローズの一般逆行列 [#x...
それでは、データ数N=3, 未知数M=4の時のムーア・ベンローズ...
その結果をプロットしてみましょう。
sageへの入力:
#pre{{
# 計画行列Φ
Phi = matrix([[ _phi(x,j) for j in range(0, (M+1))] for x...
Phi_t = Phi.transpose()
# ムーア・ベンローズの一般逆行列
PhiPhit = Phi*Phi_t
Phi_dag = Phi_t*PhiPhit.inverse(); view(Phi_dag)
}}
&ref(out4.png);
sageへの入力:
#pre{{
# 平均の重み
Wml = Phi_dag * vector(Y[:3]); Wml.apply_map(lambda x : n...
}}
sageの出力:
#pre{{
(0.100463412348611, 7.01017348792056, -5.25466361108248, ...
}}
与えられた３点を通る3次曲線が求まっていることがわかります。
sageへの入力:
#pre{{
x = var('x')
y_plt = plot(y, [x, -2, 1.5], rgbcolor='blue')
sin_plt = plot(sin(2*pi*x),[x, -2, 1.5], rgbcolor='green')
(y_plt + red_plt + sin_plt).show()
}}
&ref(out5.png);
** コメント [#f2be9b8b]
#vote(おもしろかった[2],そうでもない[0],わかりずらい[0])
皆様のご意見、ご希望をお待ちしております。
#comment_kcaptcha
ページ名:
SmartDoc
![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
