sage/text/シミュレーション
をテンプレートにして作成
[
トップ
] [
新規
|
一覧
|
単語検索
|
最終更新
|
ヘルプ
]
開始行:
[[FrontPage]]
#contents
2011/06/20からのアクセス回数 &counter;
ここで紹介したSageワークシートは、以下のURLからダウンロー...
http://www15191ue.sakura.ne.jp:8000/home/pub/10/
また、Sageのサーバを公開しているサイト（http://www.sagenb...
アップロードし、実行したり、変更していろいろ動きを試すこ...
* シミュレーション [#k053556f]
シミュレーションは、自然現象や社会現象をモデルを使って再...
コンピュータ上で専用のアプリケーションを使って行われてい...
ここでは、Sageを使って参考文献１「MathematicaによるOR」4...
をシミュレーションする方法を紹介します。
モデルの作成には、論理的なモデルを使用する場合とデータに...
データから見つける方法があります。数学モデルの場合には、...
うまく表現しているかを調べるためにデータを教師用データと...
モデルの表現能力を確認します。
** 経済データ [#b0efec6d]
モデルを作成する前に、今回使用する経済データがどのような...
経済データには、1965年から1989年までの、所得、消費、投資...
います。
それでは、Sageを使って経済データを読み込み、どのような分...
*** データの読み込み [#d5382cb3]
sageでは、豊富なPythonのライブラリを使うことができるため、
データの読み込み処理がとても楽に行えます。
今回は、csvライブラリを使ってカンマ区切りのCSV形式のデータ
(Econometrics.txt)を読み込み、
Year, Ct, Yt, It, Gtの値をリスト変数にセットします。
sageへの入力:
#pre{{
# テストデータ（CSV形式）の読み込み
import csv
fileName = DATA + 'Econometrics.txt'
csvfile = open(fileName) # CSVファイルの読み込み
reader = csv.reader(csvfile) # CSVリーダーの取得
hdr = reader.next() # ヘッダを読む
# データを読む
Year = []; CtData =[]; YtData = []; ItData = []; GtData =...
for row in reader:
Year.append(N(row[0]))
CtData.append(N(row[1]))
YtData.append(N(row[2]))
ItData.append(N(row[3]))
GtData.append(N(row[4]))
csvfile.close()
# print hdr
# print Year, YtData, CtData, ItData, GtData
}}
*** データのプロット [#q29ddda3]
以下に読み込んだ政府支出、消費、投資、所得をプロットして...
これから、消費と所得の類似性や投資の変動が大きいことが見...
sageへの入力:
#pre{{
# データの関係
gtPlot = list_plot(GtData, rgbcolor='red', figsize=3)
ctPlot = list_plot(CtData, rgbcolor='red', figsize=3)
itPlot = list_plot(ItData, rgbcolor='red', figsize=3)
ytPlot = list_plot(YtData, rgbcolor='red', figsize=3)
# プロット表示
html.table([["政府支出", "消費"], [gtPlot, ctPlot], ["投...
}}
&ref(out1.png);
** モデルの作成 [#n8865b7e]
経済データを表現するモデルとして以下のような関係をもつ数...
&ref(model.png);
矩形はモデルの変数を表し、
Y:所得、C:消費、I:投資、G:政府支出です。
矢印の方向は影響を与える向きを表します。
消費は自分自身へ矢印があり、これは前年度の消費が今年度の
消費に影響を与えていることを表しています。
このモデルを連立方程式で表すと、以下のようになります。
$$
\left\{ \begin{array}{l c l}
C_t & = & a_0 + a_1 Y_t + a_2 C_{t-1} \\
I_t & = & b_0 + b_1 Y_t \\
Y_t & = & C_t + I_t + G_t \\
\end{array} \right.
$$
モデルの作成手順は、以下のようにします。
+ 所得、投資に含まれる係数$a_0, a_1, a_2, b_0, b_1$をデー...
+ 連立方程式を解いて、消費、投資、所得を算出する関数を求...
** 係数の決定 [#y61852f9]
連立方程式のCt, Itの係数$a_0, a_1, a_2, b_0, b_1$の値を回...
(find_fit関数)を使って求めます。
教師用データとして、CtList, YtList, GtListには、1966年か...
Ct1Listには、1年ずらしたCtの値をセットします。
sageへの入力:
#pre{{
# Ctの計算には前年度のCtであるCt-1を使うため1966年から198...
CtList = CtData[1:20]
Ct1List= CtData[0:19]
YtList = YtData[1:20]
GtList = GtData[1:20]
ItList = ItData[1:20]
}}
*** 消費の係数を求める [#w25581b5]
最初に、消費Ctの式$C_t = a_0 + a_1 Y_t + a_2 C_{t-1}$
の係数$a_0, a_1, a_2$を求めてみましょう。
find_fit関数の引数は以下のように与えます。
#pre{{
find_fit(data, model, options)
data: 変数1の値, 変数2の値, ... , 変数nの値, 関数の値をタ...
model: model(変数1の値, 変数2の値, ... , 変数nの値)の引数...
options: solution_dict=Trueを指定すると係数名をキーとする...
}}
最初に、変数Yt, Ct, Ct1と係数a0, a1, a2を定義します。
次に、消費のモデルとしてCtModelシンボリック関数を定義しま...
データの形式は、消費のモデルCtModelの引数がYt, Ct1であり、
モデルの表す値がCtですから、Yt, Ct1, Ctのタプルのリストを...
各年度のYt, Ct1, Ctのタプルは、zip関数を使って作成します。
sageへの入力:
#pre{{
# 各年度のYt, Ct1, Ctをタプルにまとめたリストを作成
data = zip(YtList, Ct1List, CtList)
# モデル関数と変数を定義
(Yt, Ct, Ct1, a0, a1, a2) = var('Yt Ct Ct1 a0 a1 a2')
CtModel(Yt, Ct1) = a0 + a1*Yt + a2*Ct1
# 最適な解のa0, a1, a2を求める
CtFit = find_fit(data, CtModel, solution_dict=True); CtFit
}}
sageの出力:
#pre{{
{a2: 0.59268191617908039, a1: 0.21857976971146675, a0: 75...
}}
*** 投資の係数を求める [#q20e3942]
消費と同様に、投資Itの式$I_t = b_0 + b_1 Y_t$
の係数$b_0, b_1$を求めます。
新たに変数Yt, b0, b1を定義し、投資のモデルとしてItModelシ...
データの形式は、投資のモデルItModelの引数がYtであり、
モデルの表す値がItですから、Yt, Itのタプルのリストを渡し...
sageへの入力:
#pre{{
# 同様にItに対するb0, b1を求める
data = zip(YtList, ItList)
(Yt, b0, b1) = var('Yt b0 b1')
ItModel(Yt) = b0 + b1*Yt
ItFit = find_fit(data, ItModel, solution_dict=True); ItFit
}}
sageの出力:
#pre{{
{b1: 0.15818655082721944, b0: 311.21923398160885}
}}
** 連立方程式を解く [#w1a41a89]
係数a0, a1, b2, b0, b1が決まったので、solve関数を使って連...
を解きます。
係数を入力して式を定義する代わりに、
#pre{{
CtEq = (Ct == CtModel(Yt, Ct1)).subs(CtFit)
}}
のようにsubs関数を使ってCt == a0 + a1*Yt + a2*Ct1のa0, a1...
これで、モデルを変更しても再計算が楽にできます。
解の結果、Ct, It, YtはGtとCt1からなる式に変換されています...
経済モデルから前年度の消費と今年度の政府支出によって、今...
sageへの入力:
#pre{{
# 連立方程式をCt, It, Ytに対して解く
(Gt, It) = var('Gt It')
CtEq = (Ct == CtModel(Yt, Ct1)).subs(CtFit)
ItEq = (It == ItModel(Yt)).subs(ItFit)
YtEq = (Yt == Ct + It + Gt)
eqn = [CtEq, ItEq, YtEq]; eqn
}}
sageの出力:
#pre{{
[Ct == 0.59268191617908039*Ct1 + 0.21857976971146675*Yt +...
It == 0.15818655082721944*Yt + 311.21923398160885, Yt == ...
}}
sageへの入力:
#pre{{
sol = solve(eqn, [Ct, It, Yt], solution_dict=True); sol
}}
sageの出力:
#pre{{
[{Ct: 327042806728/408524384645*Ct1 + 59574711/169864609*...
9519709947510127/927460765140, Yt: 388497958456/408524384...
272553640/169864609*Gt + 582572871301814/46373038257,
It: 61455151728/408524384645*Ct1 + 43114320/169864609*Gt +
710582492842051/309153588380}]
}}
** 計算結果と実データの比較 [#n723b432]
計算結果と実データを比較するために、zip関数を使って前年度...
sageへの入力:
#pre{{
# 入力データをセット
data = zip(Ct1List, GtList)
}}
*** 消費・投資・所得計算関数の定義 [#a82eb912]
連立方程式から消費Ct、投資It、所得Ytを計算する式を取り出...
前年度消費Ct1と今年度政府支出を引数とする関数ctFunc, itFu...
sageへの入力:
#pre{{
ctFunc(Ct1, Gt) = sol[0][Ct]; view(ctFunc) # 消費計算...
itFunc(Ct1, Gt) = sol[0][It]; view(itFunc) # 投資計算...
ytFunc(Ct1, Gt) = sol[0][Yt]; view(ytFunc) # 所得計算...
}}
&ref(out2.png);
*** 教師用データに対するシミュレーション [#q2829564]
求めた関数ctFunc, itFunc, ytFuncを使って教師用データの消...
してみましょう。
dataの各タプルに対し、リスト内包を使って関数から計算結果...
sageへの入力:
#pre{{
# 消費のシミュレーション結果
calCtList = [ctFunc(Ct1, Gt).n() for (Ct1, Gt) in data]
# 投資のシミュレーション結果
calItList = [itFunc(Ct1, Gt) for (Ct1, Gt) in data]
# 所得のシミュレーション結果
calYtList = [ytFunc(Ct1, Gt) for (Ct1, Gt) in data]
}}
*** 教師用データに対するシミュレーション結果 [#ya7a61c8]
list_plotを使って教師用データに対するシミュレーション結果...
単純なモデルの割には消費と所得に対しては、シミュレーショ...
投資については、所得Ytのみが影響するモデルとなっているた...
sageへの入力:
#pre{{
# 計算値（青）と実測値（赤）のプロット
gtPlot = list_plot(GtList, rgbcolor='red', figsize=3)
calCtPlot = list_plot(calCtList, plotjoined=true, figsize...
ctPlot = list_plot(CtList, rgbcolor='red', figsize=3)
calItPlot = list_plot(calItList, plotjoined=true, figsize...
itPlot = list_plot(ItList, rgbcolor='red', figsize=3)
calYtPlot = list_plot(calYtList, plotjoined=true, figsize...
ytPlot = list_plot(YtList, rgbcolor='red', figsize=3)
# プロット表示
html.table([["政府支出", "消費"], [gtPlot, (calCtPlot + c...
}}
&ref(out3.png);
** 検証用データとの比較 [#lfa4a75b]
同様に検証用データに対してシミュレーションを行ってみまし...
消費、所得については、モデルはよく表現できていますが、
投資については急激な増加をうまく表現できていません。
sageへの入力:
#pre{{
# 検証用データ
valCtList = [calCtList[18]]
valItList = [calItList[18]]
valYtList = [calYtList[18]]
Ct1List = CtData[0:24]
GtList = GtData[1:25]
CtList = CtData[1:25]
for i in range(19, 24):
valCtList.append(ctFunc(Ct1List[i], GtList[i]))
valItList.append(itFunc(Ct1List[i], GtList[i]))
valYtList.append(ytFunc(Ct1List[i], GtList[i]))
}}
sageへの入力:
#pre{{
# 教師用＋検証用の実測値（赤）
ctPlot = list_plot(CtData[1:24], rgbcolor='red', figsize=3)
itPlot = list_plot(ItData[1:25], rgbcolor='red', figsize=3)
ytPlot = list_plot(YtData[1:25], rgbcolor='red', figsize=3)
gtPlot = list_plot(GtData[1:25], rgbcolor='red', figsize=3)
# 計算値（青）、検証値（灰）
valCtPlot = list_plot(zip(range(18,24), valCtList), plotj...
valItPlot = list_plot(zip(range(18,24), valItList), plotj...
valYtPlot = list_plot(zip(range(18,24), valYtList), plotj...
# プロット
ctPlt = (calCtPlot + ctPlot +valCtPlot)
itPlt = (calItPlot + itPlot + valItPlot)
ytPlt = (calYtPlot + ytPlot + valYtPlot)
html.table([["政府支出", "消費"], [gtPlot, ctPlt], ["投資...
}}
&ref(out4.png);
** 予測 [#u69fd97a]
シミュレーションの結果から、モデルは消費をうまく表現して...
このことから、もし政府支出が予測できたら今年度度消費の値...
を計算することができます。
問題を簡単にするために、政府支出Gtを$G_t(x) = c_0 + c_1 x...
sageへの入力:
#pre{{
# 1985年から1989年のデータを予測
var('x c0 c1')
data = zip(range(20), GtList[0:20])
# 直線回帰から政府支出のモデルGtの係数を求める
GtModel(x) = c0 + c1*x
GtFit = find_fit(data, GtModel, solution_dict=True); GtFit
}}
sageの出力:
#pre{{
{c1: 1837.6999989361932, c0: 36476.850068506872}
}}
*** 政府支出の予測結果 [#w03f99d7]
直線回帰で求めた政府支出と実データをプロットすると以下の...
sageへの入力:
#pre{{
gtPlot = list_plot(GtList, rgbcolor='red', figsize=3)
gtFunc(x) = GtModel.subs(GtFit)
calGtPlot = plot(gtFunc, [x, 0, 19], figsize=3)
(calGtPlot + gtPlot).show()
}}
&ref(out5.png);
*** 1967年から1989年までの予測 [#g757862e]
ちょっと乱暴な遊びですが、
政府支出関数gtFuncを使って、1967年から1989年まで政府支出...
の値とgtFuncの値から消費、投資、所得を順番に計算してみま...
sageへの入力:
#pre{{
# 予測
predCt1List =[ Ct1List[0]] + range(23)
predCtList = [calCtList[0]]
predItList = [calItList[0]]
predYtList = [calYtList[0]]
for i in range(1, 24):
predCt1List[i] = ctFunc(predCt1List[i-1], gtFunc(i-1))
predCtList.append(ctFunc(predCt1List[i], gtFunc(i)))
predItList.append(itFunc(predCt1List[i], gtFunc(i)))
predYtList.append(ytFunc(predCt1List[i], gtFunc(i)))
}}
*** 予測結果 [#a7249c61]
政府支出の乱暴な予測にも関わらず、予測値(緑)は先のシミュ...
ほぼ同じ傾向を示しています。
このようにシミュレーションは、未来の予測にも活用できるこ...
sageへの入力:
#pre{{
# 計算値（青）、予測値（緑）と実測値（赤）
predCtPlot = list_plot(zip(range(24), predCtList), plotjo...
predItPlot = list_plot(zip(range(24), predItList), plotjo...
predYtPlot = list_plot(zip(range(24), predYtList), plotjo...
calGtPlot = plot(gtFunc, [x, 0, 23], figsize=3)
# プロット
ctPlt = (calCtPlot + predCtPlot + ctPlot + valCtPlot)
itPlt = (calItPlot + predItPlot + itPlot + valItPlot)
ytPlt = (calYtPlot + predYtPlot + ytPlot + valYtPlot)
html.table([["政府支出", "消費"], [(calGtPlot + gtPlot), ...
}}
&ref(out6.png);
** コメント [#e7a9ce66]
#vote(おもしろかった[3],そうでもない[0],わかりずらい[0])
皆様のご意見、ご希望をお待ちしております。
#comment_kcaptcha
終了行:
[[FrontPage]]
#contents
2011/06/20からのアクセス回数 &counter;
ここで紹介したSageワークシートは、以下のURLからダウンロー...
http://www15191ue.sakura.ne.jp:8000/home/pub/10/
また、Sageのサーバを公開しているサイト（http://www.sagenb...
アップロードし、実行したり、変更していろいろ動きを試すこ...
* シミュレーション [#k053556f]
シミュレーションは、自然現象や社会現象をモデルを使って再...
コンピュータ上で専用のアプリケーションを使って行われてい...
ここでは、Sageを使って参考文献１「MathematicaによるOR」4...
をシミュレーションする方法を紹介します。
モデルの作成には、論理的なモデルを使用する場合とデータに...
データから見つける方法があります。数学モデルの場合には、...
うまく表現しているかを調べるためにデータを教師用データと...
モデルの表現能力を確認します。
** 経済データ [#b0efec6d]
モデルを作成する前に、今回使用する経済データがどのような...
経済データには、1965年から1989年までの、所得、消費、投資...
います。
それでは、Sageを使って経済データを読み込み、どのような分...
*** データの読み込み [#d5382cb3]
sageでは、豊富なPythonのライブラリを使うことができるため、
データの読み込み処理がとても楽に行えます。
今回は、csvライブラリを使ってカンマ区切りのCSV形式のデータ
(Econometrics.txt)を読み込み、
Year, Ct, Yt, It, Gtの値をリスト変数にセットします。
sageへの入力:
#pre{{
# テストデータ（CSV形式）の読み込み
import csv
fileName = DATA + 'Econometrics.txt'
csvfile = open(fileName) # CSVファイルの読み込み
reader = csv.reader(csvfile) # CSVリーダーの取得
hdr = reader.next() # ヘッダを読む
# データを読む
Year = []; CtData =[]; YtData = []; ItData = []; GtData =...
for row in reader:
Year.append(N(row[0]))
CtData.append(N(row[1]))
YtData.append(N(row[2]))
ItData.append(N(row[3]))
GtData.append(N(row[4]))
csvfile.close()
# print hdr
# print Year, YtData, CtData, ItData, GtData
}}
*** データのプロット [#q29ddda3]
以下に読み込んだ政府支出、消費、投資、所得をプロットして...
これから、消費と所得の類似性や投資の変動が大きいことが見...
sageへの入力:
#pre{{
# データの関係
gtPlot = list_plot(GtData, rgbcolor='red', figsize=3)
ctPlot = list_plot(CtData, rgbcolor='red', figsize=3)
itPlot = list_plot(ItData, rgbcolor='red', figsize=3)
ytPlot = list_plot(YtData, rgbcolor='red', figsize=3)
# プロット表示
html.table([["政府支出", "消費"], [gtPlot, ctPlot], ["投...
}}
&ref(out1.png);
** モデルの作成 [#n8865b7e]
経済データを表現するモデルとして以下のような関係をもつ数...
&ref(model.png);
矩形はモデルの変数を表し、
Y:所得、C:消費、I:投資、G:政府支出です。
矢印の方向は影響を与える向きを表します。
消費は自分自身へ矢印があり、これは前年度の消費が今年度の
消費に影響を与えていることを表しています。
このモデルを連立方程式で表すと、以下のようになります。
$$
\left\{ \begin{array}{l c l}
C_t & = & a_0 + a_1 Y_t + a_2 C_{t-1} \\
I_t & = & b_0 + b_1 Y_t \\
Y_t & = & C_t + I_t + G_t \\
\end{array} \right.
$$
モデルの作成手順は、以下のようにします。
+ 所得、投資に含まれる係数$a_0, a_1, a_2, b_0, b_1$をデー...
+ 連立方程式を解いて、消費、投資、所得を算出する関数を求...
** 係数の決定 [#y61852f9]
連立方程式のCt, Itの係数$a_0, a_1, a_2, b_0, b_1$の値を回...
(find_fit関数)を使って求めます。
教師用データとして、CtList, YtList, GtListには、1966年か...
Ct1Listには、1年ずらしたCtの値をセットします。
sageへの入力:
#pre{{
# Ctの計算には前年度のCtであるCt-1を使うため1966年から198...
CtList = CtData[1:20]
Ct1List= CtData[0:19]
YtList = YtData[1:20]
GtList = GtData[1:20]
ItList = ItData[1:20]
}}
*** 消費の係数を求める [#w25581b5]
最初に、消費Ctの式$C_t = a_0 + a_1 Y_t + a_2 C_{t-1}$
の係数$a_0, a_1, a_2$を求めてみましょう。
find_fit関数の引数は以下のように与えます。
#pre{{
find_fit(data, model, options)
data: 変数1の値, 変数2の値, ... , 変数nの値, 関数の値をタ...
model: model(変数1の値, 変数2の値, ... , 変数nの値)の引数...
options: solution_dict=Trueを指定すると係数名をキーとする...
}}
最初に、変数Yt, Ct, Ct1と係数a0, a1, a2を定義します。
次に、消費のモデルとしてCtModelシンボリック関数を定義しま...
データの形式は、消費のモデルCtModelの引数がYt, Ct1であり、
モデルの表す値がCtですから、Yt, Ct1, Ctのタプルのリストを...
各年度のYt, Ct1, Ctのタプルは、zip関数を使って作成します。
sageへの入力:
#pre{{
# 各年度のYt, Ct1, Ctをタプルにまとめたリストを作成
data = zip(YtList, Ct1List, CtList)
# モデル関数と変数を定義
(Yt, Ct, Ct1, a0, a1, a2) = var('Yt Ct Ct1 a0 a1 a2')
CtModel(Yt, Ct1) = a0 + a1*Yt + a2*Ct1
# 最適な解のa0, a1, a2を求める
CtFit = find_fit(data, CtModel, solution_dict=True); CtFit
}}
sageの出力:
#pre{{
{a2: 0.59268191617908039, a1: 0.21857976971146675, a0: 75...
}}
*** 投資の係数を求める [#q20e3942]
消費と同様に、投資Itの式$I_t = b_0 + b_1 Y_t$
の係数$b_0, b_1$を求めます。
新たに変数Yt, b0, b1を定義し、投資のモデルとしてItModelシ...
データの形式は、投資のモデルItModelの引数がYtであり、
モデルの表す値がItですから、Yt, Itのタプルのリストを渡し...
sageへの入力:
#pre{{
# 同様にItに対するb0, b1を求める
data = zip(YtList, ItList)
(Yt, b0, b1) = var('Yt b0 b1')
ItModel(Yt) = b0 + b1*Yt
ItFit = find_fit(data, ItModel, solution_dict=True); ItFit
}}
sageの出力:
#pre{{
{b1: 0.15818655082721944, b0: 311.21923398160885}
}}
** 連立方程式を解く [#w1a41a89]
係数a0, a1, b2, b0, b1が決まったので、solve関数を使って連...
を解きます。
係数を入力して式を定義する代わりに、
#pre{{
CtEq = (Ct == CtModel(Yt, Ct1)).subs(CtFit)
}}
のようにsubs関数を使ってCt == a0 + a1*Yt + a2*Ct1のa0, a1...
これで、モデルを変更しても再計算が楽にできます。
解の結果、Ct, It, YtはGtとCt1からなる式に変換されています...
経済モデルから前年度の消費と今年度の政府支出によって、今...
sageへの入力:
#pre{{
# 連立方程式をCt, It, Ytに対して解く
(Gt, It) = var('Gt It')
CtEq = (Ct == CtModel(Yt, Ct1)).subs(CtFit)
ItEq = (It == ItModel(Yt)).subs(ItFit)
YtEq = (Yt == Ct + It + Gt)
eqn = [CtEq, ItEq, YtEq]; eqn
}}
sageの出力:
#pre{{
[Ct == 0.59268191617908039*Ct1 + 0.21857976971146675*Yt +...
It == 0.15818655082721944*Yt + 311.21923398160885, Yt == ...
}}
sageへの入力:
#pre{{
sol = solve(eqn, [Ct, It, Yt], solution_dict=True); sol
}}
sageの出力:
#pre{{
[{Ct: 327042806728/408524384645*Ct1 + 59574711/169864609*...
9519709947510127/927460765140, Yt: 388497958456/408524384...
272553640/169864609*Gt + 582572871301814/46373038257,
It: 61455151728/408524384645*Ct1 + 43114320/169864609*Gt +
710582492842051/309153588380}]
}}
** 計算結果と実データの比較 [#n723b432]
計算結果と実データを比較するために、zip関数を使って前年度...
sageへの入力:
#pre{{
# 入力データをセット
data = zip(Ct1List, GtList)
}}
*** 消費・投資・所得計算関数の定義 [#a82eb912]
連立方程式から消費Ct、投資It、所得Ytを計算する式を取り出...
前年度消費Ct1と今年度政府支出を引数とする関数ctFunc, itFu...
sageへの入力:
#pre{{
ctFunc(Ct1, Gt) = sol[0][Ct]; view(ctFunc) # 消費計算...
itFunc(Ct1, Gt) = sol[0][It]; view(itFunc) # 投資計算...
ytFunc(Ct1, Gt) = sol[0][Yt]; view(ytFunc) # 所得計算...
}}
&ref(out2.png);
*** 教師用データに対するシミュレーション [#q2829564]
求めた関数ctFunc, itFunc, ytFuncを使って教師用データの消...
してみましょう。
dataの各タプルに対し、リスト内包を使って関数から計算結果...
sageへの入力:
#pre{{
# 消費のシミュレーション結果
calCtList = [ctFunc(Ct1, Gt).n() for (Ct1, Gt) in data]
# 投資のシミュレーション結果
calItList = [itFunc(Ct1, Gt) for (Ct1, Gt) in data]
# 所得のシミュレーション結果
calYtList = [ytFunc(Ct1, Gt) for (Ct1, Gt) in data]
}}
*** 教師用データに対するシミュレーション結果 [#ya7a61c8]
list_plotを使って教師用データに対するシミュレーション結果...
単純なモデルの割には消費と所得に対しては、シミュレーショ...
投資については、所得Ytのみが影響するモデルとなっているた...
sageへの入力:
#pre{{
# 計算値（青）と実測値（赤）のプロット
gtPlot = list_plot(GtList, rgbcolor='red', figsize=3)
calCtPlot = list_plot(calCtList, plotjoined=true, figsize...
ctPlot = list_plot(CtList, rgbcolor='red', figsize=3)
calItPlot = list_plot(calItList, plotjoined=true, figsize...
itPlot = list_plot(ItList, rgbcolor='red', figsize=3)
calYtPlot = list_plot(calYtList, plotjoined=true, figsize...
ytPlot = list_plot(YtList, rgbcolor='red', figsize=3)
# プロット表示
html.table([["政府支出", "消費"], [gtPlot, (calCtPlot + c...
}}
&ref(out3.png);
** 検証用データとの比較 [#lfa4a75b]
同様に検証用データに対してシミュレーションを行ってみまし...
消費、所得については、モデルはよく表現できていますが、
投資については急激な増加をうまく表現できていません。
sageへの入力:
#pre{{
# 検証用データ
valCtList = [calCtList[18]]
valItList = [calItList[18]]
valYtList = [calYtList[18]]
Ct1List = CtData[0:24]
GtList = GtData[1:25]
CtList = CtData[1:25]
for i in range(19, 24):
valCtList.append(ctFunc(Ct1List[i], GtList[i]))
valItList.append(itFunc(Ct1List[i], GtList[i]))
valYtList.append(ytFunc(Ct1List[i], GtList[i]))
}}
sageへの入力:
#pre{{
# 教師用＋検証用の実測値（赤）
ctPlot = list_plot(CtData[1:24], rgbcolor='red', figsize=3)
itPlot = list_plot(ItData[1:25], rgbcolor='red', figsize=3)
ytPlot = list_plot(YtData[1:25], rgbcolor='red', figsize=3)
gtPlot = list_plot(GtData[1:25], rgbcolor='red', figsize=3)
# 計算値（青）、検証値（灰）
valCtPlot = list_plot(zip(range(18,24), valCtList), plotj...
valItPlot = list_plot(zip(range(18,24), valItList), plotj...
valYtPlot = list_plot(zip(range(18,24), valYtList), plotj...
# プロット
ctPlt = (calCtPlot + ctPlot +valCtPlot)
itPlt = (calItPlot + itPlot + valItPlot)
ytPlt = (calYtPlot + ytPlot + valYtPlot)
html.table([["政府支出", "消費"], [gtPlot, ctPlt], ["投資...
}}
&ref(out4.png);
** 予測 [#u69fd97a]
シミュレーションの結果から、モデルは消費をうまく表現して...
このことから、もし政府支出が予測できたら今年度度消費の値...
を計算することができます。
問題を簡単にするために、政府支出Gtを$G_t(x) = c_0 + c_1 x...
sageへの入力:
#pre{{
# 1985年から1989年のデータを予測
var('x c0 c1')
data = zip(range(20), GtList[0:20])
# 直線回帰から政府支出のモデルGtの係数を求める
GtModel(x) = c0 + c1*x
GtFit = find_fit(data, GtModel, solution_dict=True); GtFit
}}
sageの出力:
#pre{{
{c1: 1837.6999989361932, c0: 36476.850068506872}
}}
*** 政府支出の予測結果 [#w03f99d7]
直線回帰で求めた政府支出と実データをプロットすると以下の...
sageへの入力:
#pre{{
gtPlot = list_plot(GtList, rgbcolor='red', figsize=3)
gtFunc(x) = GtModel.subs(GtFit)
calGtPlot = plot(gtFunc, [x, 0, 19], figsize=3)
(calGtPlot + gtPlot).show()
}}
&ref(out5.png);
*** 1967年から1989年までの予測 [#g757862e]
ちょっと乱暴な遊びですが、
政府支出関数gtFuncを使って、1967年から1989年まで政府支出...
の値とgtFuncの値から消費、投資、所得を順番に計算してみま...
sageへの入力:
#pre{{
# 予測
predCt1List =[ Ct1List[0]] + range(23)
predCtList = [calCtList[0]]
predItList = [calItList[0]]
predYtList = [calYtList[0]]
for i in range(1, 24):
predCt1List[i] = ctFunc(predCt1List[i-1], gtFunc(i-1))
predCtList.append(ctFunc(predCt1List[i], gtFunc(i)))
predItList.append(itFunc(predCt1List[i], gtFunc(i)))
predYtList.append(ytFunc(predCt1List[i], gtFunc(i)))
}}
*** 予測結果 [#a7249c61]
政府支出の乱暴な予測にも関わらず、予測値(緑)は先のシミュ...
ほぼ同じ傾向を示しています。
このようにシミュレーションは、未来の予測にも活用できるこ...
sageへの入力:
#pre{{
# 計算値（青）、予測値（緑）と実測値（赤）
predCtPlot = list_plot(zip(range(24), predCtList), plotjo...
predItPlot = list_plot(zip(range(24), predItList), plotjo...
predYtPlot = list_plot(zip(range(24), predYtList), plotjo...
calGtPlot = plot(gtFunc, [x, 0, 23], figsize=3)
# プロット
ctPlt = (calCtPlot + predCtPlot + ctPlot + valCtPlot)
itPlt = (calItPlot + predItPlot + itPlot + valItPlot)
ytPlt = (calYtPlot + predYtPlot + ytPlot + valYtPlot)
html.table([["政府支出", "消費"], [(calGtPlot + gtPlot), ...
}}
&ref(out6.png);
** コメント [#e7a9ce66]
#vote(おもしろかった[3],そうでもない[0],わかりずらい[0])
皆様のご意見、ご希望をお待ちしております。
#comment_kcaptcha
ページ名:
SmartDoc
![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
