sage/sageの紹介
をテンプレートにして作成
[
トップ
] [
新規
|
一覧
|
単語検索
|
最終更新
|
ヘルプ
]
開始行:
[[FrontPage]]
#contents
2009/12/02からのアクセス回数 &counter;
ここで紹介したSageワークシートは、以下のURLからダウンロー...
http://sage.math.canterbury.ac.nz/home/pub/89/
また、Sageのサーバを公開しているサイト（http://sage.math....
アップロードし、実行したり、変更していろいろ動きを試すこ...
* sageとは [#ra3b3502]
sageは、mathematicaのような数式処理を行うオープンソースの...
&ref(notebook.jpg);
sageの目標は、商用のMagma, Maple, MathematicaそしてMatlab...
** sageの誕生 [#lb35db09]
sageは、ウィリアム・スタイン (William Stein)によって2005...
最新のバージョンは、4.2.1（2009-11-14）です。
** sageに含まれるオープンソースのコンポーネント [#q5ca65e7]
sageに含まれるオープンソースのコンポーネントは、以下のURL...
http://www.sagemath.org/links-components.html
主なコンポーネントは、以下の通りです。((William Steinの20...
|計算処理|GMP, MPFR|
|可換環論|Singular (libcf, libfactor y)|
|暗号処理|OpenSSL, PyOpenSSL, PyCrypto|
|組み合わせ|GAP|
|グラフ理論|NetworkX|
|数論|PARI, NTL|
|数値計算|GSL, Numpy|
|計算機代数|Maxima|
|統計処理|R|
|特殊数値計算|多くの C/C++個別プログラム|
|コマンドライン処理|IPython|
|グラフィックインタフェース|Notebook, jsmath, Moin wiki|
|プロット処理|Matplotlib, Tachyon, libgd|
|ネットワーク処理|Twisted|
|データベース|ZODB, Python Pickles|
|プログラミング言語|Python, SageX (compiled python)|
* 基礎編 [#xe37bb90]
** インストールなしで使えるsage [#acf91637]
Sageのおもしろいところは、sageをインストールしなくてもオ...
Sageのホームページ(http://www.sagemath.org/)のTry Sage On...
ログインが完了すると以下のようなNotebook画面になります。
&ref(online.jpg);
** ノートブックの使い方 [#eb354044]
** 表記方法 [#e706e24d]
*** 変数の宣言と値の代入 [#l6a544f0]
sageの式で変数として認識させるには変数をvar関数で宣言しな...
変数の宣言は、
var('変数名')
複数の変数を宣言する場合には、スペースを空けて指定します...
x, y = var('x y')
のようにpython変数x, yに宣言したsage変数を代入します。
変数の値を指定して、関数の返す値を調べるには、
f1(x=1.5)
のように式を表す変数にカッコを付けて、x=1.5と値を指定する...
sage入力：
#pre{{
x = var('x')
f1 = (x - 1)*(x^2 - 1)
f1(x=1.5)
}}
sage出力：
#pre{{
0.625000000000000
}}
変数をx, yと複数存在する場合の例をです。f5にf5(x, y)=x^2+...
sage入力：
#pre{{
var('x y w')
f5=x^2 + y
f5(x=2, y=3)
}}
sage出力：
#pre{{
7
}}
*** 規則の代入 [#kcfd0935]
規則の代入には、subs_exprを使います。
sage入力：
#pre{{
f6=x^2 + 1; f6
}}
sage出力：
#pre{{
x^2 + 1
}}
sage入力：
#pre{{
f6.subs_expr(x^2 == w)
}}
sage出力：
#pre{{
w + 1
}}
*** 関数の定義 [#p480044a]
sageで関数を定義する方法が２つあります。
- pythonの関数を定義する
- lambda式を使って関数を定義する
です。
*** python関数の定義 [#lbff504c]
手続き的な関数を定義する場合には、python関数を使うと便利...
pythonでは
#pre{{
def 関数名(引数のリスト):
関数の本体
}}
のように定義します。
以下にx^3を返す関数cubeをpythonの関数を使って定義します。
sage入力：
#pre{{
def cube(x):
return x^3
cube(2)
}}
sage出力：
#pre{{
8
}}
sage入力：
#pre{{
x, y = var('x y')
cube(x+y)
}}
sage出力：
#pre{{
(x + y)^3
}}
*** lambda関数の定義 [#od48df9e]
つぎにlambda式を使って関数を定義する方法を示します。 x^3...
#pre{{
lambda 変数のリスト: 式
のように定義します。 先ほどと同じx3を返す関数の例を以下に...
}}
sage入力：
#pre{{
f = lambda x: x^3
f(2)
}}
sage出力：
#pre{{
f(x+y)
}}
** 基本的な計算 [#i2fc0a51]
** グラフ表示 [#w338af36]
*** グラフの使い方 [#ibafe2d2]
sageでのグラフの使い方について、説明します。
レファレンスマニュアル (( http://sagemath.com/doc/referen...
*** 基本図形 [#b11b16c5]
計算結果の表示の他に、補足説明などのために基本図形を表示...
- 円: circle
- 文字列: text
- 線: line
- 点: point
- ポリゴン: poligon
を示します。
*** 円 [#l0500bd1]
円は以下のように表示します。 circle((座標), 半径) circle...
sage入力：
#pre{{
circle((0,0), 1)
}}
sage出力：
&ref(basic-graph-1.png);
*** 文字列 [#c6a68586]
次に文字列textです。表示する文字列には、$で囲んでlatexの...
textは以下の形式で使用します。文字列の中心が指定した座標...
text(文字列, (座標))
textの例を以下に示します。原点近くで日本語文字が化けてい...
sage入力：
#pre{{
text('test', (1, 1)) + text('$f(x) = x^2 + 1$', (0.5,0.5)...
}}
sage出力：
&ref(basic-graph-2.png);
*** 線 [#baa25291]
線（line）は、指定された座標のリストを線で結びます。 line...
line([(開始座標), (終了座標)])
*** 点 [#cccbe72f]
これまで、説明しないで使ってきた点（point）です。
point((座標), オプション属性（pointsize, rgccolor, facet...
以下に例を示します。日本語が使えないため、タイトルをhtml...
sage入力：
#pre{{
html('<center>テスト</center>')
c = circle((0.5,0.5), 1)
l = line([(0,0), (1, 1)])
pt = point((0.5, 0.5), rgbcolor='white', pointsize=30, fa...
(c + l + pt).show(xmin=-1, xmax=2)
}}
sage出力：
&ref(basic-graph-3.png);
*** ポリゴンの塗りつぶし [#ga802d84]
polygon関数を使うとリストで指定した座標の図形を塗りつぶし...
sage入力：
#pre{{
polygon([(0,0), (1,1), (0,1)])
}}
sage出力：
&ref(basic-graph-4.png);
*** 2次元グラフ [#u9239084]
sageでは、簡単に2次元グラフを表示することができます。
2次元グラフには、plot関数を使用します。 plot関数のもっと...
plot(関数, 最小値, 最大値)
です。
例として、\( y=cos(x) \) のグラフを\(−2 \pi \) から\(2 \p...
sage入力：
#pre{{
plot(cos, -2*pi, 2*pi)
}}
sage出力：
&ref(basic-graph-5.png);
*** グラフの基本 [#zc8a0692]
もう少しplotの使い方を調べてみましょう。
plotの呼び出しは、以下の形式で覚えると便利です。
plot(関数, [変数名, 最小, 最大], オプション)
オプションは、省略可能です。plotのオプションは、plot.opti...
- 描画範囲指定のxmin, xmax, ymin, ymax
- グラフの比率を指定するaspect_ratio
- 線の色指定のrgbcolor
などです。
プロットのオプションは、以下のように知ることができます。
sage入力：
#pre{{
# plotのオプションを知る
plot.options
}}
sage出力：
#pre{{
{'fillalpha': 0.5, 'detect_poles': False, 'plot_points': ...
'adaptive_tolerance': 0.01, 'fillcolor': 'automatic','alp...
'rgbcolor': (0, 0, 1), 'fill':None}
}}
sage入力：
#pre{{
# Graphicsのオプションを知る
Graphics.show?
}}
sage出力：
&ref(basic-graph-6.png);
*** 重ね書き [#w8b71a83]
重ね書きの例として、直線に○をプロットしてみます。
関数 \( f(x)=\frac{(x2−1)}{(x−1)} \)を-1から3の範囲で表...
描画範囲は、(0, 0)から(3, 4)とします。
sage入力：
#pre{{
# 重ね書きの例
p = plot((x^2 - 1)/(x - 1), x, -1, 3)
pt = point((1, 2), rgbcolor='white', pointsize=30, facete...
g = p+pt
g.show(xmin=0, ymin=0, xmax=3, ymax=4)
}}
sage出力：
&ref(basic-graph-7.png);
*** plotはlimitを計算する [#n6f1a188]
グラフをみて、おやっと思われた方もいると思いますが、\(f(x...
sage入力：
#pre{{
limit((x^2 - 1)/(x - 1), x=1)
}}
sage出力：
#pre{{
2
}}
もう一つ0で不連続な関数\( \frac{sin(x)}{x} \)をグラフにし...
sage入力：
#pre{{
plot(sin(x)/x, x, -100, 100)
}}
sage出力：
&ref(basic-graph-8.png);
*** 場合分けのグラフ [#v47bd963]
以下のような場合分けのグラフを表示する例を示します。
$$
f(x)=\left\{
\begin{array}{l l}
x^2, & 0 \le x \le 1 \\
2 -x , & 1 \lt x \le 2 \\
x^2-3x+2& 2 \lt x \le 3 \\
\end{array}
\right.
$$
sage入力：
#pre{{
p1 = plot(x^2, x, 0, 1)
p2 = plot(-x+2, x, 1, 2)
p3 = plot(x^2-3*x+2, x, 2, 3)
pt1 = point((0, 0), rgbcolor='black', pointsize=30)
pt2 = point((3, 2), rgbcolor='black', pointsize=30)
(p1+p2+p3+pt1+pt2).show(xmin=0, xmax=3, ymin=0, ymax=2)
}}
sage出力：
&ref(basic-graph-9.png);
*** パラメトリック方程式の可視化 [#j7f615f0]
高校で習う運動力学で円運動がありますが、これを時間tをパラ...
$$
x=cos(t), y=sin(t)
$$
となります。xyを時間tをパラメータでので、媒介変数表示と呼...
媒介変数表示を行う関数が、parametric_plot関数です。
parametric_plot([座標を示す関数リスト], (パラメータ名, ...
の形式で使います。 例）単位円上の円運動を媒介変数表示で表...
sage入力：
#pre{{
var('t')
parametric_plot([cos(t),sin(t)],[t,0,2*pi], aspect_ratio=1)
}}
sage出力：
&ref(basic-graph-10.png);
*** リストプロット [#g5dabb05]
リスト (( []で括られた値のリスト )) をプロットするlist_pl...
list_plotの使い方は、以下の通りです。
list_plot(プロットするリスト)
sage入力：
#pre{{
list_plot([1, 2, 4, 3, 6])
}}
sage出力：
&ref(basic-graph-11.png);
*** 関係式の値をプロット [#tb4cc34f]
これまでのプロットは、プロットする値を明示的に示すもので...
sage入力：
#pre{{
implicit_plot(x^2 + y^2 == 1, [x, -1, 1], [y, -1, 1], asp...
}}
sage出力：
&ref(basic-graph-12.png);
*** 3次元グラフ [#wf265d38]
2次元グラフと同様に3次元のグラフを表示することができます...
sage入力：
#pre{{
x, y = var('x y')
plot3d(sin(x*y),(x,-pi,pi),(y,-pi,pi), mesh=True)
}}
sage出力：
&ref(basic-graph-13.png);
*** 等高線グラフ [#tecd33c4]
前のグラフと同じものを等高線グラフで表示すると以下のよう...
sage入力：
#pre{{
var('x y')
contour_plot(sin(x*y), [x, -pi, pi], [y, -pi, pi], aspect...
}}
sage出力：
&ref(basic-graph-14.png);
*** 3次元の媒介変数表示 [#a916339f]
3次元の媒介変数表示の例です。
$$
\begin{eqnarray}
f_x &=& u v \\
f_y &=& u \\
f_z &=& v^2
\end{eqnarray}
$$
を媒介変数表示したのが、以下のグラフです。
sage入力：
#pre{{
x, y = var('x y')
u, v = var('u v')
fx = u*v
fy = u
fz = v^2
parametric_plot3d([fx, fy, fz], (u, -1, 1), (v, -1, 1), f...
}}
sage出力：
&ref(basic-graph-15.png);
** 方程式 [#vad58931]
*** 方程式の解法 [#w92b4c70]
数式処理で助かる機能は、方程式の解法でしょう。
方程式の解法は、solve関数を使って以下のように行います。
solve(解きたい方程式また方程式のリスト, 解を得る変数)
例として、\(f(x)=ax+b\)の一次方程式を解いてみます。結果か...
sage入力：
#pre{{
# 方程式の解法
var('x a b')
f = a*x + b
sln = solve(f, x); sln
}}
sage出力：
#pre{{
[x == -b/a]
}}
値を得るためには、
sage入力：
#pre{{
x = -b/a
x(a=2, b=1)
}}
sage出力：
#pre{{
-1/2
}}
としなくてはなりません。
これだと、解から式をもう一度入力しなくてはなりません。そ...
sage入力：
#pre{{
var('x a b')
g = lambda x, a, b : a*x + b
solve(g(x, a, b), x)
}}
sage出力：
#pre{{
[x == -b/a]
}}
sage入力：
#pre{{
solve(g(x, a=2, b=1), x)
}}
sage出力：
#pre{{
[x == (-1/2)]
}}
解の値がほしいときには、find_root関数で数値解析で値を計算...
#pre{{
find_root(関数, 計算する変数の開始値、終了値)
}}
先ほどの関数g(x):a=2,b=1の解をx=−1,1 の範囲で数値解析しま...
sage入力：
#pre{{
find_root(g(x, a=2, b=1), -1, 1)
}}
sage出力：
#pre{{
-0.5
}}
*** 多項式の解 [#d89f8a47]
一次方程式ではあまりに簡単なので、多項式の解を求めてみま...
$$
f(x) = x^3 -2 x + 4
$$
のプロット図を以下に表示します。x軸と交差しているのはx=−2...
sage入力：
#pre{{
var('x')
f = x^3-2*x+4
plot(f, -2.5, 2.5)
}}
sage出力：
&ref(basic-eq-1.png);
solveの結果、\(x=−2, x=1−i, x=1+i\) を得ました。
sage入力：
#pre{{
solve(f, x)
}}
sage出力：
#pre{{
[x == (-I + 1), x == (I + 1), x == -2]
}}
これは、因数分解をした結果とも合致します。
sage入力：
#pre{{
factor(f)
}}
sage出力：
#pre{{
(x + 2)*(x^2 - 2*x + 2)
}}
plotで可視化することでfind_rootでの計算範囲を大まかに把握...
図から-3から3の間に解があることが見て取れたので、この範囲...
sage入力：
#pre{{
find_root(f, -3, 3)
}}
sage出力：
#pre{{
-1.9999999999999949
}}
*** 連立方程式の解 [#o75c9554]
同様に関連式のリストをsolveに与えることで、連立方程式の解...
次の連立方程式を使って、solveで連立方程式を解いてみましょ...
$$
\begin{eqnarray}
x^2 - 2y &=& 2 \\
x - y &=& 1
\end{eqnarray}
$$
sage入力：
#pre{{
var('x y')
f = [y == 1/2*x^2 - 1, y == x - 1]
solve(f, x, y)
}}
sage出力：
#pre{{
[[x == 2, y == 1], [x == 0, y == -1]]
}}
上記関数のプロット結果は、以下のとおりです。
sage入力：
#pre{{
p1 = plot(lambda x : 1/2*x^2 - 1, (x, -1, 3))
p2 = plot(lambda x : x - 1, (x, -1, 3))
(p1+p2).show()
}}
sage出力：
&ref(basic-eq-2.png);
もう一つ、図化によって簡単に解の見つかる例を見てみましょ...
$$
\begin{eqnarray}
x^2 + y^2 &=& 1 \\
y &=& x - 1
\end{eqnarray}
$$
sage入力：
#pre{{
var('x y')
f = [x^2 + y^2 == 1, y == x - 1]
solve(f, x, y)
}}
sage出力：
#pre{{
[[x == 1, y == 0], [x == 0, y == -1]]
}}
この関数をimplicit_plotを使って関係式の形をそのまま使って...
このようにsolveの解と図化を使った解のチェックを簡単に行え...
sage入力：
#pre{{
var('x y')
p1 = implicit_plot(x^2 + y^2 == 1, (x, -1, 1), (y, -1, 1))
p2 = plot(lambda x : x - 1, (x, -1, 1))
(p1+p2).show(aspect_ratio=1)
}}
sage出力：
&ref(basic-eq-3.png);
** 微分・積分 [#mba88d14]
*** 導関数 [#k366b624]
高校で習った、導関数をsageを使って導いてみましょう。
関数f(x)を\( f(x)=\frac{1}{2} x^3 \)とします。
平均変化率は、g(x)=hf(x+h)−f(x) を以下のように定義します。
sage入力：
#pre{{
# 導関数
def f(x):
return (x^3)/2
var('x, h')
g=(f(x + h) - f(x))/h; g
}}
sage出力：
#pre{{
1/2*((h + x)^3 - x^3)/h
}}
g(x)を展開し、整理すると
sage入力：
#pre{{
g1= g.rational_simplify(); g1
}}
sage出力：
#pre{{
1/2*h^2 + 3/2*h*x + 3/2*x^2
}}
となります。
ここで、h→0 の極値を取ったものが求める導関数\(f′(x)\)です。
sage入力：
#pre{{
limit(g1, h=0)
}}
sage出力：
#pre{{
3/2*x^2
}}
結果は、f(x)をxで微分したものと同じになります。
sage入力：
#pre{{
diff(f(x), x)
}}
sage出力：
#pre{{
3/2*x^2
}}
*** 微分 [#g74eae6d]
高校で習う微分の公式は、
- \(\frac{d}{dx} c = 0\)
- \(\frac{d}{dx} x^n = n x^{x-1}\)
です。
sageの微分diff関数は、以下のように使います。
#pre{{
diff(関数, 微分する変数)
}}
上記公式をsageで計算させると、以下のようになります。
sage入力：
#pre{{
# 微分
x, c, n = var('x c n')
f = function('f', x)
g = function('g', x)
print diff(c, x)
print diff(x^n, x)
}}
sage出力：
#pre{{
0
n*x^(n - 1)
}}
*** いろんな関数の微分 [#v3657bbb]
以下の関数を微分した結果を示します。
- \( (x^2+1)^3 \)
- \( sin(x)^3 \)
- \( sin^{-1}(x) \)
- \( e^x \)
- \( log(x) \)
sage入力：
#pre{{
# いろんな関数の微分
print diff((x^2+1)^3).factor()
print diff(sin(x)^3, x)
print diff(arcsin(x), x)
print diff(exp(x), x)
print diff(log(x), x)
}}
sage出力：
#pre{{
6*(x^2 + 1)^2*x
3*sin(x)^2*cos(x)
1/sqrt(-x^2 + 1)
e^x
1/x
}}
*** 微分の応用 [#f7474b10]
微分の応用として、接線の方程式とその法線を計算してみまし...
$$
y−y_1=f′(x_1)(x−x_1)
$$
となります。 また、法線は、
$$
y−y_1=\frac{1}{f′(x1)}(x−x_1)
$$
となります。
微分した式をgとし、x=1での接線と法線を計算し、プロットし...
sage入力：
#pre{{
f = exp(-(x^2)); f
}}
sage出力：
#pre{{
e^(-x^2)
}}
sage入力：
#pre{{
g = diff(f, x); g
}}
sage出力：
#pre{{
-2*x*e^(-x^2)
}}
mにx=1でのgの値をセット
sage入力：
#pre{{
m = g(x=1); m
}}
sage出力：
#pre{{
-2*e^(-1)
}}
関数fと接線、法線をプロットすると、
sage入力：
#pre{{
p1 = plot(f, (x, -2, 2))
p2 = plot(m*(x-1)+f(1), (x, -2, 2))
p3 = plot(-(x-1)/m+f(1), (x, -2, 2))
pt = point([1, f(1)])
(p1+p2+p3+pt).show(aspect_ratio=1, ymin=-1, ymax=2)
}}
sage出力：
&ref(basic-diff-1.png);
*** 積分 [#xcd5a3b2]
sageでは積分は、integrate関数で計算します。
#pre{{
integrate(被積分関数, 積分変数)
また
integrate(被積分関数, 積分変数, 積分範囲)
とすると定積分を計算します。
}}
関数f = xを積分し、その後に求めた積分関数を微分してみます。
sage入力：
#pre{{
f = integrate(x); f
}}
sage出力：
#pre{{
1/2*x^2
}}
sage入力：
#pre{{
diff(f)
}}
sage出力：
#pre{{
1
}}
*** いろんな関数の積分 [#mcba0e3e]
以下の関数をsageを使って積分した結果を示します。
- \( \int cdx \)定数c
- \( \int \frac{1}{x}dx \)
- \( \int e^x dx \)
- \( \int log(x)dx \)
- \( \int sin(x)dx \)
- \( \int \frac{f'(x)}{f(x)}dx \)
sage入力：
#pre{{
# いろんな関数の積分
x = var('x')
c = var('c')
f = function('f', x)
g = function('g', x)
print integrate(c, x)
print integrate(1/x, x)
print integrate(exp(x), x)
print integrate(log(x), x)
print integrate(sin(x), x)
print integrate(diff(f, x)/f(x), x)
}}
sage出力：
#pre{{
c*x
log(x)
e^x
x*log(x) - x
-cos(x)
log(f(x))
}}
*** 定積分 [#nc1bef84]
次に定積分\( \int_0^{\pi/2}sin(x)dx \)の計算を示します。
sage入力：
#pre{{
integrate(sin(x), x, 0, pi/2)
}}
sage出力：
#pre{{
1
}}
もう一つ
$$
\int_0^3 x^2 + 2 sin(x) dx
$$
もやってみましょう。
sage入力：
#pre{{
integrate(x^2+2*sin(x), x, 0, 3)
}}
sage出力：
#pre{{
-2*cos(3) + 11
}}
上記結果を数値とすると以下のようになります。関数Nは数値に...
sage入力：
#pre{{
N(_)
}}
sage出力：
#pre{{
12.9799849932009
}}
*** 数値積分 [#v480ee4b]
数値積分は、numerical_integral関数を使って計算します。
numerical_integral(被関分館数, 積分範囲)
で戻り値は、積分結果と誤差を返します。
以下は、\( sin(x) \)を0から\( \pi \)で数値積分した結果で...
sage入力：
#pre{{
sigma, error = numerical_integral(sin(x), 0, pi); (sigma...
}}
sage出力：
#pre{{
(1.9999999999999998, 2.2204460492503128e-14)
}}
*** 積分の応用 [#kb365ad6]
積分の応用例として、サイクロイドの計算をしてみます。
サイクロイドは、tをパラメータとして、以下のような式で表さ...
$$
\begin{eqnarray}
x &=& 2(t-sin(t)) \\
y &=& 2(1-cos(t))
\end{eqnarray}
$$
この曲線の長さを計算すると、
$$
\int_{0}^{2\pi}\sqrt{\left(\frac{dx}{dt}\right)^2+\left(\...
$$
となります。解析解では、この結果は16となります。
sage入力：
#pre{{
t = var('t')
x = 2*(t-sin(t))
y = 2*(1-cos(t))
parametric_plot([x, y], (t, 0, 2*pi))
}}
sage出力：
&ref(basic-int-1.png);
残念ながら、定積分は結果がでませんでした。
sage入力：
#pre{{
cycloid = sqrt(diff(x,t)^2+diff(y,t)^2)
integrate(cycloid, t, 0, 2*pi)
}}
sage出力：
#pre{{
integrate(sqrt(4*(cos(t) - 1)^2 + 4*sin(t)^2), t, 0, 2*pi)
}}
それで、数値積分で計算し、期待した結果となりました。
sage入力：
#pre{{
numerical_integral(cycloid, 0, 2*pi)
}}
sage出力：
#pre{{
(15.999999999999998, 1.7763568394002502e-13)
}}
** ベクトルと行列 [#e5768075]
* 応用編 [#pc4d1af6]
** インタラクティブな計算 [#h1486b4e]
** 微分方程式を解く [#sed60b4d]
** 経済モデルを解く [#iebce591]
この例題は、以下のURLでノートブックを参照することができま...
http://www.sagenb.org/pub/1182/
*** 計量経済モデルとは [#g31c7c4e]
計量経済学では、経済の構造をモデルで表現し、モデルの要素...
ここでは、「MathematicaによるOR」4章で紹介されているモデ...
以下のモデルでは、矩形を要素としY:所得、C:消費、I:投資、 ...
&ref(eco-0.png);
このモデルを連立方程式で表すと、以下のようになります。
$$
\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 + Y_t + G_t \\
\end{array}
$$
*** データの読み込み [#pb6c3cf1]
sageでは、統計情報の分析に有益なRパッケージとのインタフェ...
経済データは、Econometrics.txtファイルに入っています。 最...
sage入力：
#pre{{
# 計量経済学(Econometrics)
# Rのread.table関数を使ってデータを読み込む
fileName = DATA + 'Econometrics.txt'
d = r("read.table('%s', header=T)" %fileName) ; d
}}
sage出力：
#pre{{
Year Ct Yt It Gt
1 1965 57176.9 89662.8 10325.1 26408.3
2 1966 63042.7 100033.0 12874.9 29014.0
3 1967 69239.4 110974.0 16394.3 32237.9
4 1968 75771.6 125997.0 19831.2 36869.6
5 1969 83092.1 141402.0 25595.8 39998.9
6 1970 88847.4 153915.0 28909.3 45352.6
7 1971 94196.4 161688.0 27715.0 47546.6
8 1972 103848.0 176628.0 29409.4 53396.2
9 1973 110290.0 184569.0 33396.1 56225.9
10 1974 111694.0 183798.0 30345.5 52550.3
11 1975 115765.0 190875.0 29288.2 53985.8
12 1976 120366.0 199630.0 29478.3 56413.9
13 1977 125394.0 210235.0 29722.5 60068.7
14 1978 133154.0 221243.0 32393.5 64567.0
15 1979 140184.0 232878.0 35569.1 65560.4
16 1980 141398.0 242131.0 38293.4 63846.2
17 1981 144257.0 250159.0 39917.0 64394.4
18 1982 150360.0 258241.0 40698.0 64297.5
19 1983 154910.0 267700.0 42711.7 63274.3
20 1984 158910.0 281399.0 47631.2 64695.8
21 1985 163395.0 293982.0 53899.5 64404.0
22 1986 169016.0 301846.0 56236.0 68598.6
23 1987 176556.0 318109.0 61897.9 74350.9
24 1988 185318.0 334969.0 72617.1 76545.4
25 1989 191233.0 351877.0 84584.1 77764.1
}}
読み込んだデータdをsageのデータ形式に変換するには、_sage_...
sage入力：
#pre{{
# 読み込んだデータdをsageのデータに変換
ds = d._sage_()
# Ct, Yt, It, Gtを取り出す
CtData = ds['DATA']['Ct']
YtData = ds['DATA']['Yt']
ItData = ds['DATA']['It']
GtData = ds['DATA']['Gt']
}}
*** 係数の決定 [#j9af8fcf]
連立方程式のCt, Itの係数a_0, a_1, a_2, b_0, b_1の値を回帰...
CtList, YtList, GtListには、1966年から1984年のデータをセ...
sage入力：
#pre{{
# Ctのフィッティングは、前年度のCtであるCt-1を使うため196...
CtList = CtData[1:20]
Ct1List= CtData[0:19]
YtList = YtData[1:20]
GtList = GtData[1:20]
ItList = ItData[1:20]
}}
*** 回帰分析関数find_fitの使い方 [#e690ffb0]
回帰分析には関数find_fitを使います。
#pre{{
find_fit(data, model, options)
dataは、変数1の値, 変数2の値, ... , 変数nの値, 関数の値を...
modelは、model(変数1の値, 変数2の値, ... , 変数nの値)の引...
optionsで、solution_dict=Trueを指定すると係数名をキーとす...
}}
zip関数を使って各年度のYt, Ct1, Ctをタプルにまとめたリス...
sage入力：
#pre{{
# 各年度のYt, Ct1, Ctをタプルにまとめたリストを作成
data = zip(YtList, Ct1List, CtList); data
}}
sage出力：
#pre{{
[(100033, 57176.900000000001, 63042.699999999997), (110974,
63042.699999999997, 69239.399999999994), (125997,
69239.399999999994, 75771.600000000006), (141402,
75771.600000000006, 83092.100000000006), (153915,
83092.100000000006, 88847.399999999994), (161688,
88847.399999999994, 94196.399999999994), (176628,
94196.399999999994, 103848), (184569, 103848, 110290), (1...
110290, 111694), (190875, 111694, 115765), (199630, 11576...
(210235, 120366, 125394), (221243, 125394, 133154), (2328...
140184), (242131, 140184, 141398), (250159, 141398, 14425...
(258241, 144257, 150360), (267700, 150360, 154910), (2813...
158910)]
}}
*** 消費の係数を求める [#ieeaf898]
変数Yt, Ct1, a0, a1, a2とモデルCtModelを定義します。CtMod...
\( C_t=a0+a_1 Y_t+a_2 C_{t−1} \)の部分です。
- a0: 7518.3096037375299
- a1: 0.21857976971146531
- a2: 0.59268191617908206
が求まります。
sage入力：
#pre{{
# モデル関数と変数を定義
var('Yt 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{{
{a0: 7518.3096037375299, a1: 0.21857976971146531, a2:
0.59268191617908206}
}}
*** 投資の係数を求める [#h8f4edf2]
同様に、変数Yt, b0, b1とモデルItModelを定義します。CtMode...
- b0: -11736.272856134668
- b1: 0.22718320402830394
が求まります。
sage入力：
#pre{{
# 同様にItに対するb0, b1を求める
data = zip(YtData, ItData)
var('Yt b0 b1')
ItModel(Yt) = b0 + b1*Yt
ItFit = find_fit(data, ItModel, solution_dict=True); ItFit
}}
sage出力：
#pre{{
{b0: -11736.272856134668, b1: 0.22718320402830394}
}}
*** 連立方程式を解く [#z5451f63]
係数a0, a1, b2, b0, b1が決まったので、solve関数を使って連...
係数を入力して式を定義する代わりに、
#pre{{
CtEq = (Ct == CtModel(Yt, Ct1)).subs(CtFit)
}}
のようにsubs関数を使ってCt == a0 + a1*Yt + a2*Ct1のa0, a1...
sage入力：
#pre{{
# 連立方程式をCt, It, Ytに対して解く
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.592681916179082*Ct1 + 0.218579769711465*Yt +
7518.30960373753, It == 0.227183204028304*Yt - 11736.2728...
== Ct + Gt + It]
}}
solveにsolution_dict=Trueのオプションを指定し、辞書形式で...
sage入力：
#pre{{
sol = solve(eqn, [Ct, It, Yt], solution_dict=True); sol
}}
sage出力
#pre{{
[{Yt: 721680419112/674867923015*Ct1 + 506300280/280610363...
164437808222620/21606997951, Ct: 557726748944/67486792301...
110666997/280610363*Gt + 379515962234357/64820993853, It:
163953670168/674867923015*Ct1 + 115022920/280610363*Gt -
872829386902217/64820993853}]
}}
*** 計算結果と実データの比較 [#b6ac6d60]
前年度のzip関数を使って前年度消費と政府支出をタプルとする...
sage入力：
#pre{{
# 入力データをセット
data = zip(Ct1List, GtList); data
}}
sage出力：
#pre{{
[(57176.900000000001, 29014), (63042.699999999997,
32237.900000000001), (69239.399999999994, 36869.599999999...
(75771.600000000006, 39998.900000000001), (83092.10000000...
45352.599999999999), (88847.399999999994, 47546.599999999...
(94196.399999999994, 53396.199999999997), (103848,
56225.900000000001), (110290, 52550.300000000003), (111694,
53985.800000000003), (115765, 56413.900000000001), (120366,
60068.699999999997), (125394, 64567), (133154, 65560.3999...
(140184, 63846.199999999997), (141398, 64394.400000000001...
64297.5), (150360, 63274.300000000003), (154910,
64695.800000000003)]
}}
*** 消費データの比較 [#s25a6962]
連立方程式からCt式を取り出し、前年度消費Ct1と政府支出Gtを...
sage入力：
#pre{{
ctFunc(Ct1, Gt) = sol[0][Ct]; ctFunc
}}
sage出力：
#pre{{
(Ct1, Gt) |--> 557726748944/674867923015*Ct1 +
110666997/280610363*Gt + 379515962234357/64820993853
}}
*** 消費データ計算 [#zc1d9aad]
dataの各タプルに対し、リスト内for文を使って計算結果をcalC...
sage入力：
#pre{{
calCtList = [ctFunc(Ct1, Gt).n() for (Ct1, Gt) in data]; ...
}}
sage出力：
#pre{{
[64549.6971260734, 70668.7727491564, 77616.5194992103,
84249.0154719521, 92410.2390710246, 98031.8233620791,
104759.325753524, 113851.611245606, 117725.850212143,
119452.280771693, 123774.244126530, 129017.997396518,
134947.292700448, 141752.116171432, 146885.828296538,
148105.305379392, 150429.834945473, 155069.968723451,
159390.806476880]
}}
*** 消費データのプロット [#qc97165a]
list_plotを使って計算結果（青い線）と実データ（赤い点）を...
かなりよい精度で計算できていることが分かります。
sage入力：
#pre{{
# 計算値（青）と実測値（赤）
calCtPlot = list_plot(calCtList, plotjoined=true)
ctPlot = list_plot(CtList, rgbcolor='red')
(calCtPlot + ctPlot).show(ymin=0)
}}
sage出力：
&ref(eco-1.png);
*** 投資データのプロット [#q97d81a7]
同様に投資データも計算し、プロットします。
こちらはかなりずれていますが、傾向はつかんでいる感じです...
sage入力：
#pre{{
# 同様にItに対する計算値と実測値をプロット
itFunc(Ct1, Gt) = sol[0][It]
calItList = [itFunc(Ct1, Gt) for (Ct1, Gt) in data]
calItPlot = list_plot(calItList, plotjoined=true)
itPlot = list_plot(ItList, rgbcolor='red')
(calItPlot + itPlot).show(ymin=0)
}}
sage出力：
&ref(eco-2.png);
*** 所得データのプロット [#ed5af4f0]
同様に所得データも計算し、プロットします。
よく一致しています。
sage入力：
#pre{{
# 同様にYtに対する計算値と実測値をプロット
ytFunc(Ct1, Gt) = sol[0][Yt]
calYtList = [ytFunc(Ct1, Gt) for (Ct1, Gt) in data]
calYtPlot = list_plot(calYtList, plotjoined=true)
ytPlot = list_plot(YtList, rgbcolor='red')
(calYtPlot + ytPlot).show(ymin=0)
}}
sage出力：
&ref(eco-3.png);
*** 未来の予測 [#w162823b]
計量経済モデルの目標である、未来の予測をしてみましょう。
1985年から1989年までの政府支出Gtをプロットすると、直線上...
sage入力：
#pre{{
# 予測：Gtを予測し、それに対するIt, Yt, Ctを計算する
gtPlot = list_plot(GtList, rgbcolor='red'); gtPlot.show()
}}
sage出力：
&ref(eco-4.png);
sage入力：
#pre{{
# 1985年から1989年のデータを予測
var('x c0 c1')
data = zip(range(20), GtList[0:20])
GtModel(x) = c0 + c1*x
GtFit = find_fit(data, GtModel, solution_dict=True); GtFit
}}
sage出力：
#pre{{
{c1: 1960.3166656122035, c0: 35741.149903570993}
}}
*** 政府支出の予測結果 [#tc7a9709]
直線回帰で求めた政府支出と実データをプロットすると以下の...
sage入力：
#pre{{
gtFunc(x) = GtModel.subs(GtFit)
calGtPlot = plot(gtFunc, [x, 0, 25])
(calGtPlot + gtPlot).show()
}}
sage出力：
&ref(eco-5.png);
*** 1985年から1989年までの予測 [#s5a1b5e6]
政府支出モデルgtFuncを使って、1985年から1989年まで政府支...
*** 消費の予測 [#ge027c50]
予測した消費を1985年から1989年までの実データと重ねてプロ...
sage入力：
#pre{{
# 予測
predCt1List = Ct1List
predCtList = calCtList
predItList = calItList
predYtList = calYtList
for i in range(19, 24):
predCt1List.append(ctFunc(predCt1List[i-1], gtFunc(i)))
predCtList.append(ctFunc(predCt1List[i], gtFunc(i)))
predItList.append(itFunc(predCt1List[i], gtFunc(i)))
predYtList.append(ytFunc(predCt1List[i], gtFunc(i)))
# 予測値（青）と実測値（赤）
predCtPlot = list_plot(predCtList, plotjoined=true)
ctPlot = list_plot(CtData[1:25], rgbcolor='red')
(predCtPlot + ctPlot).show(ymin=0)
}}
sage出力：
&ref(eco-6.png);
*** 投資の予測 [#bac8062b]
同様に投資の予測と実データを重ねてプロットします。 計算の...
sage入力：
#pre{{
# 同様にItに対する予測値と実測値をプロット
predItPlot = list_plot(predItList, plotjoined=true)
itPlot = list_plot(ItData[1:25], rgbcolor='red')
(predItPlot + itPlot).show(ymin=0)
}}
sage出力：
&ref(eco-7.png);
*** 所得の予測 [#wdd6ddb3]
最後に所得の予測と実データを重ねてプロットします。 こちら...
このようにsageを使うと簡単に計量経済モデルの計算をするこ...
sage入力：
#pre{{
# 同様にYtに対する予測値と実測値をプロット
predYtPlot = list_plot(predYtList, plotjoined=true)
ytPlot = list_plot(YtData[1:25], rgbcolor='red')
(predYtPlot + ytPlot).show(ymin=0)
}}
sage出力：
&ref(eco-8.png);
* インストール方法 [#ef739fbc]
** ダウンロード [#e6577808]
sageのホームページ http://www.sagemath.org/ で、Download...
現在、sageがサポートしているOSは、
- MacOSX (10.4 or 10.5)
- Linux (Debian, Fedora)
- Solaris
- Windows (VM PlayerとVMイメージを使ってサポート)
です。
** MacOSX(10.5)へのインストール [#g80700a5]
MacOSXのインストールは、
- sage-4.1.1-OSX10.5-Intel-32bit-i386-Darwin.dmgのディス...
をダウンロードし、をマウント（ダブルクリック）するとsage...
ここでは、ホームディレクトリのlocalフォルダ以下にコピーし...
次にコピーしたsageフォルダ内のsageアイコンをダブルクリッ...
#pre{{
sage:
}}
ここで、notebook()と入力するとブラウザにSageのNotebook画...
- New Worksheetをクリックして、
plot(cos, -5, 5)
と入力し、シフトキーとリターンキーを同時に押すと最初にご...
** Windowsへのインストール [#f734f8bc]
*** VM Playerのインストール [#z66259f3]
Windowsの場合、最初にVM Playerをインストールしてください...
http://www.vmware.com/products/player/
*** VMイメージのダウンロード [#ga53c28f]
次にVMイメージファイルをダウンロードします
- VMWare Image of Sageをクリックし、イメージファイル（778...
をダウンロードし、解凍します。
*** VMイメージファイルの起動 [#cd20681f]
VM Playerで
- sage_vmx.vmxファイルを起動
すると、以下のような画面がでます。
&ref(vm_snap.jpg);
sage login: プロンプトにnotebookと入力すると、
&ref(vm_login.jpg);
が表示されますので、
Open Firefox to the address http://172.16.137.131
の部分に記載されたURL（例ではhttp://172.16.137.131）をブ...
&ref(vm_notebook.jpg);
の画面がでたら、成功です。
*** sageの終了 [#l818ad3c]
VM Playerの画面で、Ctrl-Cを入力し、sage login: 画面になっ...
終了行:
[[FrontPage]]
#contents
2009/12/02からのアクセス回数 &counter;
ここで紹介したSageワークシートは、以下のURLからダウンロー...
http://sage.math.canterbury.ac.nz/home/pub/89/
また、Sageのサーバを公開しているサイト（http://sage.math....
アップロードし、実行したり、変更していろいろ動きを試すこ...
* sageとは [#ra3b3502]
sageは、mathematicaのような数式処理を行うオープンソースの...
&ref(notebook.jpg);
sageの目標は、商用のMagma, Maple, MathematicaそしてMatlab...
** sageの誕生 [#lb35db09]
sageは、ウィリアム・スタイン (William Stein)によって2005...
最新のバージョンは、4.2.1（2009-11-14）です。
** sageに含まれるオープンソースのコンポーネント [#q5ca65e7]
sageに含まれるオープンソースのコンポーネントは、以下のURL...
http://www.sagemath.org/links-components.html
主なコンポーネントは、以下の通りです。((William Steinの20...
|計算処理|GMP, MPFR|
|可換環論|Singular (libcf, libfactor y)|
|暗号処理|OpenSSL, PyOpenSSL, PyCrypto|
|組み合わせ|GAP|
|グラフ理論|NetworkX|
|数論|PARI, NTL|
|数値計算|GSL, Numpy|
|計算機代数|Maxima|
|統計処理|R|
|特殊数値計算|多くの C/C++個別プログラム|
|コマンドライン処理|IPython|
|グラフィックインタフェース|Notebook, jsmath, Moin wiki|
|プロット処理|Matplotlib, Tachyon, libgd|
|ネットワーク処理|Twisted|
|データベース|ZODB, Python Pickles|
|プログラミング言語|Python, SageX (compiled python)|
* 基礎編 [#xe37bb90]
** インストールなしで使えるsage [#acf91637]
Sageのおもしろいところは、sageをインストールしなくてもオ...
Sageのホームページ(http://www.sagemath.org/)のTry Sage On...
ログインが完了すると以下のようなNotebook画面になります。
&ref(online.jpg);
** ノートブックの使い方 [#eb354044]
** 表記方法 [#e706e24d]
*** 変数の宣言と値の代入 [#l6a544f0]
sageの式で変数として認識させるには変数をvar関数で宣言しな...
変数の宣言は、
var('変数名')
複数の変数を宣言する場合には、スペースを空けて指定します...
x, y = var('x y')
のようにpython変数x, yに宣言したsage変数を代入します。
変数の値を指定して、関数の返す値を調べるには、
f1(x=1.5)
のように式を表す変数にカッコを付けて、x=1.5と値を指定する...
sage入力：
#pre{{
x = var('x')
f1 = (x - 1)*(x^2 - 1)
f1(x=1.5)
}}
sage出力：
#pre{{
0.625000000000000
}}
変数をx, yと複数存在する場合の例をです。f5にf5(x, y)=x^2+...
sage入力：
#pre{{
var('x y w')
f5=x^2 + y
f5(x=2, y=3)
}}
sage出力：
#pre{{
7
}}
*** 規則の代入 [#kcfd0935]
規則の代入には、subs_exprを使います。
sage入力：
#pre{{
f6=x^2 + 1; f6
}}
sage出力：
#pre{{
x^2 + 1
}}
sage入力：
#pre{{
f6.subs_expr(x^2 == w)
}}
sage出力：
#pre{{
w + 1
}}
*** 関数の定義 [#p480044a]
sageで関数を定義する方法が２つあります。
- pythonの関数を定義する
- lambda式を使って関数を定義する
です。
*** python関数の定義 [#lbff504c]
手続き的な関数を定義する場合には、python関数を使うと便利...
pythonでは
#pre{{
def 関数名(引数のリスト):
関数の本体
}}
のように定義します。
以下にx^3を返す関数cubeをpythonの関数を使って定義します。
sage入力：
#pre{{
def cube(x):
return x^3
cube(2)
}}
sage出力：
#pre{{
8
}}
sage入力：
#pre{{
x, y = var('x y')
cube(x+y)
}}
sage出力：
#pre{{
(x + y)^3
}}
*** lambda関数の定義 [#od48df9e]
つぎにlambda式を使って関数を定義する方法を示します。 x^3...
#pre{{
lambda 変数のリスト: 式
のように定義します。 先ほどと同じx3を返す関数の例を以下に...
}}
sage入力：
#pre{{
f = lambda x: x^3
f(2)
}}
sage出力：
#pre{{
f(x+y)
}}
** 基本的な計算 [#i2fc0a51]
** グラフ表示 [#w338af36]
*** グラフの使い方 [#ibafe2d2]
sageでのグラフの使い方について、説明します。
レファレンスマニュアル (( http://sagemath.com/doc/referen...
*** 基本図形 [#b11b16c5]
計算結果の表示の他に、補足説明などのために基本図形を表示...
- 円: circle
- 文字列: text
- 線: line
- 点: point
- ポリゴン: poligon
を示します。
*** 円 [#l0500bd1]
円は以下のように表示します。 circle((座標), 半径) circle...
sage入力：
#pre{{
circle((0,0), 1)
}}
sage出力：
&ref(basic-graph-1.png);
*** 文字列 [#c6a68586]
次に文字列textです。表示する文字列には、$で囲んでlatexの...
textは以下の形式で使用します。文字列の中心が指定した座標...
text(文字列, (座標))
textの例を以下に示します。原点近くで日本語文字が化けてい...
sage入力：
#pre{{
text('test', (1, 1)) + text('$f(x) = x^2 + 1$', (0.5,0.5)...
}}
sage出力：
&ref(basic-graph-2.png);
*** 線 [#baa25291]
線（line）は、指定された座標のリストを線で結びます。 line...
line([(開始座標), (終了座標)])
*** 点 [#cccbe72f]
これまで、説明しないで使ってきた点（point）です。
point((座標), オプション属性（pointsize, rgccolor, facet...
以下に例を示します。日本語が使えないため、タイトルをhtml...
sage入力：
#pre{{
html('<center>テスト</center>')
c = circle((0.5,0.5), 1)
l = line([(0,0), (1, 1)])
pt = point((0.5, 0.5), rgbcolor='white', pointsize=30, fa...
(c + l + pt).show(xmin=-1, xmax=2)
}}
sage出力：
&ref(basic-graph-3.png);
*** ポリゴンの塗りつぶし [#ga802d84]
polygon関数を使うとリストで指定した座標の図形を塗りつぶし...
sage入力：
#pre{{
polygon([(0,0), (1,1), (0,1)])
}}
sage出力：
&ref(basic-graph-4.png);
*** 2次元グラフ [#u9239084]
sageでは、簡単に2次元グラフを表示することができます。
2次元グラフには、plot関数を使用します。 plot関数のもっと...
plot(関数, 最小値, 最大値)
です。
例として、\( y=cos(x) \) のグラフを\(−2 \pi \) から\(2 \p...
sage入力：
#pre{{
plot(cos, -2*pi, 2*pi)
}}
sage出力：
&ref(basic-graph-5.png);
*** グラフの基本 [#zc8a0692]
もう少しplotの使い方を調べてみましょう。
plotの呼び出しは、以下の形式で覚えると便利です。
plot(関数, [変数名, 最小, 最大], オプション)
オプションは、省略可能です。plotのオプションは、plot.opti...
- 描画範囲指定のxmin, xmax, ymin, ymax
- グラフの比率を指定するaspect_ratio
- 線の色指定のrgbcolor
などです。
プロットのオプションは、以下のように知ることができます。
sage入力：
#pre{{
# plotのオプションを知る
plot.options
}}
sage出力：
#pre{{
{'fillalpha': 0.5, 'detect_poles': False, 'plot_points': ...
'adaptive_tolerance': 0.01, 'fillcolor': 'automatic','alp...
'rgbcolor': (0, 0, 1), 'fill':None}
}}
sage入力：
#pre{{
# Graphicsのオプションを知る
Graphics.show?
}}
sage出力：
&ref(basic-graph-6.png);
*** 重ね書き [#w8b71a83]
重ね書きの例として、直線に○をプロットしてみます。
関数 \( f(x)=\frac{(x2−1)}{(x−1)} \)を-1から3の範囲で表...
描画範囲は、(0, 0)から(3, 4)とします。
sage入力：
#pre{{
# 重ね書きの例
p = plot((x^2 - 1)/(x - 1), x, -1, 3)
pt = point((1, 2), rgbcolor='white', pointsize=30, facete...
g = p+pt
g.show(xmin=0, ymin=0, xmax=3, ymax=4)
}}
sage出力：
&ref(basic-graph-7.png);
*** plotはlimitを計算する [#n6f1a188]
グラフをみて、おやっと思われた方もいると思いますが、\(f(x...
sage入力：
#pre{{
limit((x^2 - 1)/(x - 1), x=1)
}}
sage出力：
#pre{{
2
}}
もう一つ0で不連続な関数\( \frac{sin(x)}{x} \)をグラフにし...
sage入力：
#pre{{
plot(sin(x)/x, x, -100, 100)
}}
sage出力：
&ref(basic-graph-8.png);
*** 場合分けのグラフ [#v47bd963]
以下のような場合分けのグラフを表示する例を示します。
$$
f(x)=\left\{
\begin{array}{l l}
x^2, & 0 \le x \le 1 \\
2 -x , & 1 \lt x \le 2 \\
x^2-3x+2& 2 \lt x \le 3 \\
\end{array}
\right.
$$
sage入力：
#pre{{
p1 = plot(x^2, x, 0, 1)
p2 = plot(-x+2, x, 1, 2)
p3 = plot(x^2-3*x+2, x, 2, 3)
pt1 = point((0, 0), rgbcolor='black', pointsize=30)
pt2 = point((3, 2), rgbcolor='black', pointsize=30)
(p1+p2+p3+pt1+pt2).show(xmin=0, xmax=3, ymin=0, ymax=2)
}}
sage出力：
&ref(basic-graph-9.png);
*** パラメトリック方程式の可視化 [#j7f615f0]
高校で習う運動力学で円運動がありますが、これを時間tをパラ...
$$
x=cos(t), y=sin(t)
$$
となります。xyを時間tをパラメータでので、媒介変数表示と呼...
媒介変数表示を行う関数が、parametric_plot関数です。
parametric_plot([座標を示す関数リスト], (パラメータ名, ...
の形式で使います。 例）単位円上の円運動を媒介変数表示で表...
sage入力：
#pre{{
var('t')
parametric_plot([cos(t),sin(t)],[t,0,2*pi], aspect_ratio=1)
}}
sage出力：
&ref(basic-graph-10.png);
*** リストプロット [#g5dabb05]
リスト (( []で括られた値のリスト )) をプロットするlist_pl...
list_plotの使い方は、以下の通りです。
list_plot(プロットするリスト)
sage入力：
#pre{{
list_plot([1, 2, 4, 3, 6])
}}
sage出力：
&ref(basic-graph-11.png);
*** 関係式の値をプロット [#tb4cc34f]
これまでのプロットは、プロットする値を明示的に示すもので...
sage入力：
#pre{{
implicit_plot(x^2 + y^2 == 1, [x, -1, 1], [y, -1, 1], asp...
}}
sage出力：
&ref(basic-graph-12.png);
*** 3次元グラフ [#wf265d38]
2次元グラフと同様に3次元のグラフを表示することができます...
sage入力：
#pre{{
x, y = var('x y')
plot3d(sin(x*y),(x,-pi,pi),(y,-pi,pi), mesh=True)
}}
sage出力：
&ref(basic-graph-13.png);
*** 等高線グラフ [#tecd33c4]
前のグラフと同じものを等高線グラフで表示すると以下のよう...
sage入力：
#pre{{
var('x y')
contour_plot(sin(x*y), [x, -pi, pi], [y, -pi, pi], aspect...
}}
sage出力：
&ref(basic-graph-14.png);
*** 3次元の媒介変数表示 [#a916339f]
3次元の媒介変数表示の例です。
$$
\begin{eqnarray}
f_x &=& u v \\
f_y &=& u \\
f_z &=& v^2
\end{eqnarray}
$$
を媒介変数表示したのが、以下のグラフです。
sage入力：
#pre{{
x, y = var('x y')
u, v = var('u v')
fx = u*v
fy = u
fz = v^2
parametric_plot3d([fx, fy, fz], (u, -1, 1), (v, -1, 1), f...
}}
sage出力：
&ref(basic-graph-15.png);
** 方程式 [#vad58931]
*** 方程式の解法 [#w92b4c70]
数式処理で助かる機能は、方程式の解法でしょう。
方程式の解法は、solve関数を使って以下のように行います。
solve(解きたい方程式また方程式のリスト, 解を得る変数)
例として、\(f(x)=ax+b\)の一次方程式を解いてみます。結果か...
sage入力：
#pre{{
# 方程式の解法
var('x a b')
f = a*x + b
sln = solve(f, x); sln
}}
sage出力：
#pre{{
[x == -b/a]
}}
値を得るためには、
sage入力：
#pre{{
x = -b/a
x(a=2, b=1)
}}
sage出力：
#pre{{
-1/2
}}
としなくてはなりません。
これだと、解から式をもう一度入力しなくてはなりません。そ...
sage入力：
#pre{{
var('x a b')
g = lambda x, a, b : a*x + b
solve(g(x, a, b), x)
}}
sage出力：
#pre{{
[x == -b/a]
}}
sage入力：
#pre{{
solve(g(x, a=2, b=1), x)
}}
sage出力：
#pre{{
[x == (-1/2)]
}}
解の値がほしいときには、find_root関数で数値解析で値を計算...
#pre{{
find_root(関数, 計算する変数の開始値、終了値)
}}
先ほどの関数g(x):a=2,b=1の解をx=−1,1 の範囲で数値解析しま...
sage入力：
#pre{{
find_root(g(x, a=2, b=1), -1, 1)
}}
sage出力：
#pre{{
-0.5
}}
*** 多項式の解 [#d89f8a47]
一次方程式ではあまりに簡単なので、多項式の解を求めてみま...
$$
f(x) = x^3 -2 x + 4
$$
のプロット図を以下に表示します。x軸と交差しているのはx=−2...
sage入力：
#pre{{
var('x')
f = x^3-2*x+4
plot(f, -2.5, 2.5)
}}
sage出力：
&ref(basic-eq-1.png);
solveの結果、\(x=−2, x=1−i, x=1+i\) を得ました。
sage入力：
#pre{{
solve(f, x)
}}
sage出力：
#pre{{
[x == (-I + 1), x == (I + 1), x == -2]
}}
これは、因数分解をした結果とも合致します。
sage入力：
#pre{{
factor(f)
}}
sage出力：
#pre{{
(x + 2)*(x^2 - 2*x + 2)
}}
plotで可視化することでfind_rootでの計算範囲を大まかに把握...
図から-3から3の間に解があることが見て取れたので、この範囲...
sage入力：
#pre{{
find_root(f, -3, 3)
}}
sage出力：
#pre{{
-1.9999999999999949
}}
*** 連立方程式の解 [#o75c9554]
同様に関連式のリストをsolveに与えることで、連立方程式の解...
次の連立方程式を使って、solveで連立方程式を解いてみましょ...
$$
\begin{eqnarray}
x^2 - 2y &=& 2 \\
x - y &=& 1
\end{eqnarray}
$$
sage入力：
#pre{{
var('x y')
f = [y == 1/2*x^2 - 1, y == x - 1]
solve(f, x, y)
}}
sage出力：
#pre{{
[[x == 2, y == 1], [x == 0, y == -1]]
}}
上記関数のプロット結果は、以下のとおりです。
sage入力：
#pre{{
p1 = plot(lambda x : 1/2*x^2 - 1, (x, -1, 3))
p2 = plot(lambda x : x - 1, (x, -1, 3))
(p1+p2).show()
}}
sage出力：
&ref(basic-eq-2.png);
もう一つ、図化によって簡単に解の見つかる例を見てみましょ...
$$
\begin{eqnarray}
x^2 + y^2 &=& 1 \\
y &=& x - 1
\end{eqnarray}
$$
sage入力：
#pre{{
var('x y')
f = [x^2 + y^2 == 1, y == x - 1]
solve(f, x, y)
}}
sage出力：
#pre{{
[[x == 1, y == 0], [x == 0, y == -1]]
}}
この関数をimplicit_plotを使って関係式の形をそのまま使って...
このようにsolveの解と図化を使った解のチェックを簡単に行え...
sage入力：
#pre{{
var('x y')
p1 = implicit_plot(x^2 + y^2 == 1, (x, -1, 1), (y, -1, 1))
p2 = plot(lambda x : x - 1, (x, -1, 1))
(p1+p2).show(aspect_ratio=1)
}}
sage出力：
&ref(basic-eq-3.png);
** 微分・積分 [#mba88d14]
*** 導関数 [#k366b624]
高校で習った、導関数をsageを使って導いてみましょう。
関数f(x)を\( f(x)=\frac{1}{2} x^3 \)とします。
平均変化率は、g(x)=hf(x+h)−f(x) を以下のように定義します。
sage入力：
#pre{{
# 導関数
def f(x):
return (x^3)/2
var('x, h')
g=(f(x + h) - f(x))/h; g
}}
sage出力：
#pre{{
1/2*((h + x)^3 - x^3)/h
}}
g(x)を展開し、整理すると
sage入力：
#pre{{
g1= g.rational_simplify(); g1
}}
sage出力：
#pre{{
1/2*h^2 + 3/2*h*x + 3/2*x^2
}}
となります。
ここで、h→0 の極値を取ったものが求める導関数\(f′(x)\)です。
sage入力：
#pre{{
limit(g1, h=0)
}}
sage出力：
#pre{{
3/2*x^2
}}
結果は、f(x)をxで微分したものと同じになります。
sage入力：
#pre{{
diff(f(x), x)
}}
sage出力：
#pre{{
3/2*x^2
}}
*** 微分 [#g74eae6d]
高校で習う微分の公式は、
- \(\frac{d}{dx} c = 0\)
- \(\frac{d}{dx} x^n = n x^{x-1}\)
です。
sageの微分diff関数は、以下のように使います。
#pre{{
diff(関数, 微分する変数)
}}
上記公式をsageで計算させると、以下のようになります。
sage入力：
#pre{{
# 微分
x, c, n = var('x c n')
f = function('f', x)
g = function('g', x)
print diff(c, x)
print diff(x^n, x)
}}
sage出力：
#pre{{
0
n*x^(n - 1)
}}
*** いろんな関数の微分 [#v3657bbb]
以下の関数を微分した結果を示します。
- \( (x^2+1)^3 \)
- \( sin(x)^3 \)
- \( sin^{-1}(x) \)
- \( e^x \)
- \( log(x) \)
sage入力：
#pre{{
# いろんな関数の微分
print diff((x^2+1)^3).factor()
print diff(sin(x)^3, x)
print diff(arcsin(x), x)
print diff(exp(x), x)
print diff(log(x), x)
}}
sage出力：
#pre{{
6*(x^2 + 1)^2*x
3*sin(x)^2*cos(x)
1/sqrt(-x^2 + 1)
e^x
1/x
}}
*** 微分の応用 [#f7474b10]
微分の応用として、接線の方程式とその法線を計算してみまし...
$$
y−y_1=f′(x_1)(x−x_1)
$$
となります。 また、法線は、
$$
y−y_1=\frac{1}{f′(x1)}(x−x_1)
$$
となります。
微分した式をgとし、x=1での接線と法線を計算し、プロットし...
sage入力：
#pre{{
f = exp(-(x^2)); f
}}
sage出力：
#pre{{
e^(-x^2)
}}
sage入力：
#pre{{
g = diff(f, x); g
}}
sage出力：
#pre{{
-2*x*e^(-x^2)
}}
mにx=1でのgの値をセット
sage入力：
#pre{{
m = g(x=1); m
}}
sage出力：
#pre{{
-2*e^(-1)
}}
関数fと接線、法線をプロットすると、
sage入力：
#pre{{
p1 = plot(f, (x, -2, 2))
p2 = plot(m*(x-1)+f(1), (x, -2, 2))
p3 = plot(-(x-1)/m+f(1), (x, -2, 2))
pt = point([1, f(1)])
(p1+p2+p3+pt).show(aspect_ratio=1, ymin=-1, ymax=2)
}}
sage出力：
&ref(basic-diff-1.png);
*** 積分 [#xcd5a3b2]
sageでは積分は、integrate関数で計算します。
#pre{{
integrate(被積分関数, 積分変数)
また
integrate(被積分関数, 積分変数, 積分範囲)
とすると定積分を計算します。
}}
関数f = xを積分し、その後に求めた積分関数を微分してみます。
sage入力：
#pre{{
f = integrate(x); f
}}
sage出力：
#pre{{
1/2*x^2
}}
sage入力：
#pre{{
diff(f)
}}
sage出力：
#pre{{
1
}}
*** いろんな関数の積分 [#mcba0e3e]
以下の関数をsageを使って積分した結果を示します。
- \( \int cdx \)定数c
- \( \int \frac{1}{x}dx \)
- \( \int e^x dx \)
- \( \int log(x)dx \)
- \( \int sin(x)dx \)
- \( \int \frac{f'(x)}{f(x)}dx \)
sage入力：
#pre{{
# いろんな関数の積分
x = var('x')
c = var('c')
f = function('f', x)
g = function('g', x)
print integrate(c, x)
print integrate(1/x, x)
print integrate(exp(x), x)
print integrate(log(x), x)
print integrate(sin(x), x)
print integrate(diff(f, x)/f(x), x)
}}
sage出力：
#pre{{
c*x
log(x)
e^x
x*log(x) - x
-cos(x)
log(f(x))
}}
*** 定積分 [#nc1bef84]
次に定積分\( \int_0^{\pi/2}sin(x)dx \)の計算を示します。
sage入力：
#pre{{
integrate(sin(x), x, 0, pi/2)
}}
sage出力：
#pre{{
1
}}
もう一つ
$$
\int_0^3 x^2 + 2 sin(x) dx
$$
もやってみましょう。
sage入力：
#pre{{
integrate(x^2+2*sin(x), x, 0, 3)
}}
sage出力：
#pre{{
-2*cos(3) + 11
}}
上記結果を数値とすると以下のようになります。関数Nは数値に...
sage入力：
#pre{{
N(_)
}}
sage出力：
#pre{{
12.9799849932009
}}
*** 数値積分 [#v480ee4b]
数値積分は、numerical_integral関数を使って計算します。
numerical_integral(被関分館数, 積分範囲)
で戻り値は、積分結果と誤差を返します。
以下は、\( sin(x) \)を0から\( \pi \)で数値積分した結果で...
sage入力：
#pre{{
sigma, error = numerical_integral(sin(x), 0, pi); (sigma...
}}
sage出力：
#pre{{
(1.9999999999999998, 2.2204460492503128e-14)
}}
*** 積分の応用 [#kb365ad6]
積分の応用例として、サイクロイドの計算をしてみます。
サイクロイドは、tをパラメータとして、以下のような式で表さ...
$$
\begin{eqnarray}
x &=& 2(t-sin(t)) \\
y &=& 2(1-cos(t))
\end{eqnarray}
$$
この曲線の長さを計算すると、
$$
\int_{0}^{2\pi}\sqrt{\left(\frac{dx}{dt}\right)^2+\left(\...
$$
となります。解析解では、この結果は16となります。
sage入力：
#pre{{
t = var('t')
x = 2*(t-sin(t))
y = 2*(1-cos(t))
parametric_plot([x, y], (t, 0, 2*pi))
}}
sage出力：
&ref(basic-int-1.png);
残念ながら、定積分は結果がでませんでした。
sage入力：
#pre{{
cycloid = sqrt(diff(x,t)^2+diff(y,t)^2)
integrate(cycloid, t, 0, 2*pi)
}}
sage出力：
#pre{{
integrate(sqrt(4*(cos(t) - 1)^2 + 4*sin(t)^2), t, 0, 2*pi)
}}
それで、数値積分で計算し、期待した結果となりました。
sage入力：
#pre{{
numerical_integral(cycloid, 0, 2*pi)
}}
sage出力：
#pre{{
(15.999999999999998, 1.7763568394002502e-13)
}}
** ベクトルと行列 [#e5768075]
* 応用編 [#pc4d1af6]
** インタラクティブな計算 [#h1486b4e]
** 微分方程式を解く [#sed60b4d]
** 経済モデルを解く [#iebce591]
この例題は、以下のURLでノートブックを参照することができま...
http://www.sagenb.org/pub/1182/
*** 計量経済モデルとは [#g31c7c4e]
計量経済学では、経済の構造をモデルで表現し、モデルの要素...
ここでは、「MathematicaによるOR」4章で紹介されているモデ...
以下のモデルでは、矩形を要素としY:所得、C:消費、I:投資、 ...
&ref(eco-0.png);
このモデルを連立方程式で表すと、以下のようになります。
$$
\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 + Y_t + G_t \\
\end{array}
$$
*** データの読み込み [#pb6c3cf1]
sageでは、統計情報の分析に有益なRパッケージとのインタフェ...
経済データは、Econometrics.txtファイルに入っています。 最...
sage入力：
#pre{{
# 計量経済学(Econometrics)
# Rのread.table関数を使ってデータを読み込む
fileName = DATA + 'Econometrics.txt'
d = r("read.table('%s', header=T)" %fileName) ; d
}}
sage出力：
#pre{{
Year Ct Yt It Gt
1 1965 57176.9 89662.8 10325.1 26408.3
2 1966 63042.7 100033.0 12874.9 29014.0
3 1967 69239.4 110974.0 16394.3 32237.9
4 1968 75771.6 125997.0 19831.2 36869.6
5 1969 83092.1 141402.0 25595.8 39998.9
6 1970 88847.4 153915.0 28909.3 45352.6
7 1971 94196.4 161688.0 27715.0 47546.6
8 1972 103848.0 176628.0 29409.4 53396.2
9 1973 110290.0 184569.0 33396.1 56225.9
10 1974 111694.0 183798.0 30345.5 52550.3
11 1975 115765.0 190875.0 29288.2 53985.8
12 1976 120366.0 199630.0 29478.3 56413.9
13 1977 125394.0 210235.0 29722.5 60068.7
14 1978 133154.0 221243.0 32393.5 64567.0
15 1979 140184.0 232878.0 35569.1 65560.4
16 1980 141398.0 242131.0 38293.4 63846.2
17 1981 144257.0 250159.0 39917.0 64394.4
18 1982 150360.0 258241.0 40698.0 64297.5
19 1983 154910.0 267700.0 42711.7 63274.3
20 1984 158910.0 281399.0 47631.2 64695.8
21 1985 163395.0 293982.0 53899.5 64404.0
22 1986 169016.0 301846.0 56236.0 68598.6
23 1987 176556.0 318109.0 61897.9 74350.9
24 1988 185318.0 334969.0 72617.1 76545.4
25 1989 191233.0 351877.0 84584.1 77764.1
}}
読み込んだデータdをsageのデータ形式に変換するには、_sage_...
sage入力：
#pre{{
# 読み込んだデータdをsageのデータに変換
ds = d._sage_()
# Ct, Yt, It, Gtを取り出す
CtData = ds['DATA']['Ct']
YtData = ds['DATA']['Yt']
ItData = ds['DATA']['It']
GtData = ds['DATA']['Gt']
}}
*** 係数の決定 [#j9af8fcf]
連立方程式のCt, Itの係数a_0, a_1, a_2, b_0, b_1の値を回帰...
CtList, YtList, GtListには、1966年から1984年のデータをセ...
sage入力：
#pre{{
# Ctのフィッティングは、前年度のCtであるCt-1を使うため196...
CtList = CtData[1:20]
Ct1List= CtData[0:19]
YtList = YtData[1:20]
GtList = GtData[1:20]
ItList = ItData[1:20]
}}
*** 回帰分析関数find_fitの使い方 [#e690ffb0]
回帰分析には関数find_fitを使います。
#pre{{
find_fit(data, model, options)
dataは、変数1の値, 変数2の値, ... , 変数nの値, 関数の値を...
modelは、model(変数1の値, 変数2の値, ... , 変数nの値)の引...
optionsで、solution_dict=Trueを指定すると係数名をキーとす...
}}
zip関数を使って各年度のYt, Ct1, Ctをタプルにまとめたリス...
sage入力：
#pre{{
# 各年度のYt, Ct1, Ctをタプルにまとめたリストを作成
data = zip(YtList, Ct1List, CtList); data
}}
sage出力：
#pre{{
[(100033, 57176.900000000001, 63042.699999999997), (110974,
63042.699999999997, 69239.399999999994), (125997,
69239.399999999994, 75771.600000000006), (141402,
75771.600000000006, 83092.100000000006), (153915,
83092.100000000006, 88847.399999999994), (161688,
88847.399999999994, 94196.399999999994), (176628,
94196.399999999994, 103848), (184569, 103848, 110290), (1...
110290, 111694), (190875, 111694, 115765), (199630, 11576...
(210235, 120366, 125394), (221243, 125394, 133154), (2328...
140184), (242131, 140184, 141398), (250159, 141398, 14425...
(258241, 144257, 150360), (267700, 150360, 154910), (2813...
158910)]
}}
*** 消費の係数を求める [#ieeaf898]
変数Yt, Ct1, a0, a1, a2とモデルCtModelを定義します。CtMod...
\( C_t=a0+a_1 Y_t+a_2 C_{t−1} \)の部分です。
- a0: 7518.3096037375299
- a1: 0.21857976971146531
- a2: 0.59268191617908206
が求まります。
sage入力：
#pre{{
# モデル関数と変数を定義
var('Yt 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{{
{a0: 7518.3096037375299, a1: 0.21857976971146531, a2:
0.59268191617908206}
}}
*** 投資の係数を求める [#h8f4edf2]
同様に、変数Yt, b0, b1とモデルItModelを定義します。CtMode...
- b0: -11736.272856134668
- b1: 0.22718320402830394
が求まります。
sage入力：
#pre{{
# 同様にItに対するb0, b1を求める
data = zip(YtData, ItData)
var('Yt b0 b1')
ItModel(Yt) = b0 + b1*Yt
ItFit = find_fit(data, ItModel, solution_dict=True); ItFit
}}
sage出力：
#pre{{
{b0: -11736.272856134668, b1: 0.22718320402830394}
}}
*** 連立方程式を解く [#z5451f63]
係数a0, a1, b2, b0, b1が決まったので、solve関数を使って連...
係数を入力して式を定義する代わりに、
#pre{{
CtEq = (Ct == CtModel(Yt, Ct1)).subs(CtFit)
}}
のようにsubs関数を使ってCt == a0 + a1*Yt + a2*Ct1のa0, a1...
sage入力：
#pre{{
# 連立方程式をCt, It, Ytに対して解く
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.592681916179082*Ct1 + 0.218579769711465*Yt +
7518.30960373753, It == 0.227183204028304*Yt - 11736.2728...
== Ct + Gt + It]
}}
solveにsolution_dict=Trueのオプションを指定し、辞書形式で...
sage入力：
#pre{{
sol = solve(eqn, [Ct, It, Yt], solution_dict=True); sol
}}
sage出力
#pre{{
[{Yt: 721680419112/674867923015*Ct1 + 506300280/280610363...
164437808222620/21606997951, Ct: 557726748944/67486792301...
110666997/280610363*Gt + 379515962234357/64820993853, It:
163953670168/674867923015*Ct1 + 115022920/280610363*Gt -
872829386902217/64820993853}]
}}
*** 計算結果と実データの比較 [#b6ac6d60]
前年度のzip関数を使って前年度消費と政府支出をタプルとする...
sage入力：
#pre{{
# 入力データをセット
data = zip(Ct1List, GtList); data
}}
sage出力：
#pre{{
[(57176.900000000001, 29014), (63042.699999999997,
32237.900000000001), (69239.399999999994, 36869.599999999...
(75771.600000000006, 39998.900000000001), (83092.10000000...
45352.599999999999), (88847.399999999994, 47546.599999999...
(94196.399999999994, 53396.199999999997), (103848,
56225.900000000001), (110290, 52550.300000000003), (111694,
53985.800000000003), (115765, 56413.900000000001), (120366,
60068.699999999997), (125394, 64567), (133154, 65560.3999...
(140184, 63846.199999999997), (141398, 64394.400000000001...
64297.5), (150360, 63274.300000000003), (154910,
64695.800000000003)]
}}
*** 消費データの比較 [#s25a6962]
連立方程式からCt式を取り出し、前年度消費Ct1と政府支出Gtを...
sage入力：
#pre{{
ctFunc(Ct1, Gt) = sol[0][Ct]; ctFunc
}}
sage出力：
#pre{{
(Ct1, Gt) |--> 557726748944/674867923015*Ct1 +
110666997/280610363*Gt + 379515962234357/64820993853
}}
*** 消費データ計算 [#zc1d9aad]
dataの各タプルに対し、リスト内for文を使って計算結果をcalC...
sage入力：
#pre{{
calCtList = [ctFunc(Ct1, Gt).n() for (Ct1, Gt) in data]; ...
}}
sage出力：
#pre{{
[64549.6971260734, 70668.7727491564, 77616.5194992103,
84249.0154719521, 92410.2390710246, 98031.8233620791,
104759.325753524, 113851.611245606, 117725.850212143,
119452.280771693, 123774.244126530, 129017.997396518,
134947.292700448, 141752.116171432, 146885.828296538,
148105.305379392, 150429.834945473, 155069.968723451,
159390.806476880]
}}
*** 消費データのプロット [#qc97165a]
list_plotを使って計算結果（青い線）と実データ（赤い点）を...
かなりよい精度で計算できていることが分かります。
sage入力：
#pre{{
# 計算値（青）と実測値（赤）
calCtPlot = list_plot(calCtList, plotjoined=true)
ctPlot = list_plot(CtList, rgbcolor='red')
(calCtPlot + ctPlot).show(ymin=0)
}}
sage出力：
&ref(eco-1.png);
*** 投資データのプロット [#q97d81a7]
同様に投資データも計算し、プロットします。
こちらはかなりずれていますが、傾向はつかんでいる感じです...
sage入力：
#pre{{
# 同様にItに対する計算値と実測値をプロット
itFunc(Ct1, Gt) = sol[0][It]
calItList = [itFunc(Ct1, Gt) for (Ct1, Gt) in data]
calItPlot = list_plot(calItList, plotjoined=true)
itPlot = list_plot(ItList, rgbcolor='red')
(calItPlot + itPlot).show(ymin=0)
}}
sage出力：
&ref(eco-2.png);
*** 所得データのプロット [#ed5af4f0]
同様に所得データも計算し、プロットします。
よく一致しています。
sage入力：
#pre{{
# 同様にYtに対する計算値と実測値をプロット
ytFunc(Ct1, Gt) = sol[0][Yt]
calYtList = [ytFunc(Ct1, Gt) for (Ct1, Gt) in data]
calYtPlot = list_plot(calYtList, plotjoined=true)
ytPlot = list_plot(YtList, rgbcolor='red')
(calYtPlot + ytPlot).show(ymin=0)
}}
sage出力：
&ref(eco-3.png);
*** 未来の予測 [#w162823b]
計量経済モデルの目標である、未来の予測をしてみましょう。
1985年から1989年までの政府支出Gtをプロットすると、直線上...
sage入力：
#pre{{
# 予測：Gtを予測し、それに対するIt, Yt, Ctを計算する
gtPlot = list_plot(GtList, rgbcolor='red'); gtPlot.show()
}}
sage出力：
&ref(eco-4.png);
sage入力：
#pre{{
# 1985年から1989年のデータを予測
var('x c0 c1')
data = zip(range(20), GtList[0:20])
GtModel(x) = c0 + c1*x
GtFit = find_fit(data, GtModel, solution_dict=True); GtFit
}}
sage出力：
#pre{{
{c1: 1960.3166656122035, c0: 35741.149903570993}
}}
*** 政府支出の予測結果 [#tc7a9709]
直線回帰で求めた政府支出と実データをプロットすると以下の...
sage入力：
#pre{{
gtFunc(x) = GtModel.subs(GtFit)
calGtPlot = plot(gtFunc, [x, 0, 25])
(calGtPlot + gtPlot).show()
}}
sage出力：
&ref(eco-5.png);
*** 1985年から1989年までの予測 [#s5a1b5e6]
政府支出モデルgtFuncを使って、1985年から1989年まで政府支...
*** 消費の予測 [#ge027c50]
予測した消費を1985年から1989年までの実データと重ねてプロ...
sage入力：
#pre{{
# 予測
predCt1List = Ct1List
predCtList = calCtList
predItList = calItList
predYtList = calYtList
for i in range(19, 24):
predCt1List.append(ctFunc(predCt1List[i-1], gtFunc(i)))
predCtList.append(ctFunc(predCt1List[i], gtFunc(i)))
predItList.append(itFunc(predCt1List[i], gtFunc(i)))
predYtList.append(ytFunc(predCt1List[i], gtFunc(i)))
# 予測値（青）と実測値（赤）
predCtPlot = list_plot(predCtList, plotjoined=true)
ctPlot = list_plot(CtData[1:25], rgbcolor='red')
(predCtPlot + ctPlot).show(ymin=0)
}}
sage出力：
&ref(eco-6.png);
*** 投資の予測 [#bac8062b]
同様に投資の予測と実データを重ねてプロットします。 計算の...
sage入力：
#pre{{
# 同様にItに対する予測値と実測値をプロット
predItPlot = list_plot(predItList, plotjoined=true)
itPlot = list_plot(ItData[1:25], rgbcolor='red')
(predItPlot + itPlot).show(ymin=0)
}}
sage出力：
&ref(eco-7.png);
*** 所得の予測 [#wdd6ddb3]
最後に所得の予測と実データを重ねてプロットします。 こちら...
このようにsageを使うと簡単に計量経済モデルの計算をするこ...
sage入力：
#pre{{
# 同様にYtに対する予測値と実測値をプロット
predYtPlot = list_plot(predYtList, plotjoined=true)
ytPlot = list_plot(YtData[1:25], rgbcolor='red')
(predYtPlot + ytPlot).show(ymin=0)
}}
sage出力：
&ref(eco-8.png);
* インストール方法 [#ef739fbc]
** ダウンロード [#e6577808]
sageのホームページ http://www.sagemath.org/ で、Download...
現在、sageがサポートしているOSは、
- MacOSX (10.4 or 10.5)
- Linux (Debian, Fedora)
- Solaris
- Windows (VM PlayerとVMイメージを使ってサポート)
です。
** MacOSX(10.5)へのインストール [#g80700a5]
MacOSXのインストールは、
- sage-4.1.1-OSX10.5-Intel-32bit-i386-Darwin.dmgのディス...
をダウンロードし、をマウント（ダブルクリック）するとsage...
ここでは、ホームディレクトリのlocalフォルダ以下にコピーし...
次にコピーしたsageフォルダ内のsageアイコンをダブルクリッ...
#pre{{
sage:
}}
ここで、notebook()と入力するとブラウザにSageのNotebook画...
- New Worksheetをクリックして、
plot(cos, -5, 5)
と入力し、シフトキーとリターンキーを同時に押すと最初にご...
** Windowsへのインストール [#f734f8bc]
*** VM Playerのインストール [#z66259f3]
Windowsの場合、最初にVM Playerをインストールしてください...
http://www.vmware.com/products/player/
*** VMイメージのダウンロード [#ga53c28f]
次にVMイメージファイルをダウンロードします
- VMWare Image of Sageをクリックし、イメージファイル（778...
をダウンロードし、解凍します。
*** VMイメージファイルの起動 [#cd20681f]
VM Playerで
- sage_vmx.vmxファイルを起動
すると、以下のような画面がでます。
&ref(vm_snap.jpg);
sage login: プロンプトにnotebookと入力すると、
&ref(vm_login.jpg);
が表示されますので、
Open Firefox to the address http://172.16.137.131
の部分に記載されたURL（例ではhttp://172.16.137.131）をブ...
&ref(vm_notebook.jpg);
の画面がでたら、成功です。
*** sageの終了 [#l818ad3c]
VM Playerの画面で、Ctrl-Cを入力し、sage login: 画面になっ...
ページ名:
SmartDoc
![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
