Hiroshi TAKEMOTO (take@pwv.co.jp)

はじめに

エンジニアによって憧れのツールの一つに数式処理システムがあります。しかし、数式処理システムは、高価でかつ多機能なため簡単には使うことができませんでした。

ここでは、フリーの数式処理システムSageを例に、基礎編、応用編に分けて数式処理システムの使い方を紹介します。

基礎編：Sageの基本操作とよく使われる関数の使い方を紹介
応用編：「計量経済モデル」を例にSageを実務に使う方法を紹介

sageとは

sageは、Mathematicaのような数式処理を行うオープンソースのソフトウェアです。 Maxima, R等の既存のオープンソースパッケージをPythonをベースとしたインタフェースで結合し、ブラウザー（Firefox等）から簡単に操作できるノートブックを提供しています。(図参照)

sageの目標は、商用のMagma, Maple, MathematicaそしてMatlabの代わりとなるフリーかつオープンソースのシステムを開発することです。

sageの誕生

sageは、2005年2月にウィリアム・スタイン氏(William Stein)によって開発がスタートし、 2006年2月のUCSD SAGE Days 1でSage 1.0が公開されました注1。

sageに含まれるオープンソースのコンポーネント

sageに含まれる主なコンポーネントは、以下の通りです注2。

計算処理	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)

基礎編

インストールなしで使えるsage

Sageのおもしろいところは、sageをインストールしなくてもオンラインでsageを使うことができるところです。

Sageのホームページ(http://www.sagemath.org/)のTry Sage Onlineをクリックして、 Sign up for a new Sage Notebook accountでアカウントを作ったら早速使ってみよう！

ログインが完了すると以下のようなNotebook画面になります。

ノートブックの使い方

Sageを使うときには、クイックレファレンス注3を印刷して傍らに置いておくと便利です。

ワークシートの操作

ワークシートの基本操作は、以下のように行います。

新規作成：Notebook画面で、New Worksheet リンクをクリックすると新しいワークシートが作成されます。
オープン：Notebook画面で開きたいワークシートをクリックすると新しいウィンドウにワークシートが表示されます。
削除：Notebook画面で、削除したいワークシートの左側のチェックボックスにマークを付け、上部のDeleteボタンをクリックするとワークシートが削除されます。

セルの操作

ワークシートで数式を評価するには、セルと呼ばれるテキストエリアを利用します。
セルの操作には、以下の種類があります。

セルの評価：セルに記述した数式を評価するには、シフトキーとエンターキーを同時に押します（shit-enterと記す）。またはevaluateリンクをクリックします。
セルの追加（分割）：カーソル位置で、コントロールキーと；（セミコロン）を同時に押すとその位置でセルが分割されます。またはセルの上部にマウスを移動したときに表示されるブルーラインをクリックするとセルが挿入されます
セルの削除：空のセルでバックスペースキーを押すとセルが削除されます。

ヘルプ

sageの関数は、関数名の後に?を付けてshit-enterでヘルプが表示されます。?を2個付けるとさらに詳しいヘルプが表示されます。

abs?

補完機能

関数名やpythonの変数の属性は、タブキーで補完することができます。

以下の例では、plを入力した後にタブキーを入力すると、以下のような補完の候補が表示されます。マウスまたはカーソルで選択してください。

表記方法

sageを実行する際に必要な表記に関する約束事のうち、代表的なものを以下に示します。

python文法

sageは、pythonのインタプリタを使っているため、pythonの構文がそのまま使えます。

pythonでのコメント（実行に影響しない注釈）は、#記号から行末がコメントとして扱われます。

# 1 + 2 この行はコメントです。 1 + 2 + 3

複数の式を1行にまとめる

複数の式を1行にまとめるには、式を;(カンマ)で区切ります。

1 + 2 + 3; 4 + 5 # 複数の式を1行にまとめるには、式を;(カンマ)で区切ります。 a = 1; a # 代入の結果は、値を返さないので;で区切って変数名を書いて代入結果を表示します。

python変数への代入

計算結果をpython変数に代入するときには、= (イコール)を使います。 python変数と断っているのは、sageの数式で使われる変数と区別するためです。

変数aに$\frac{1}{2} + \frac{1}{3}$の結果を代入します。結果は、$\frac{5}{6}$と分数で表示されます。 aの値を数値として出力するには、Nまたはn関数を使います。

a = 1/2 + 1/3; print a N(a)

5/6
0.833333333333333

5/6
0.833333333333333

タプル

タプルは、カンマで区切られた値からなるシーケンスデータ型です。タプルの出力は()カッコで括られたカンマ区切りの形式で表示されます。

タプルの変数への代入をタプル・パッキングと呼び、逆にタプルから変数への代入をシーケンス・アンパッキングと呼びます。

t = 1, 2, 3; print t # 変数ｔへのタプル(1, 2, 3)の代入（タプル・パッキング） x, y, z = t; print x, y, z # タプルtから変数x, y, zへの代入（シーケンス・アンパッキング）

(1, 2, 3)
1 2 3

(1, 2, 3)
1 2 3

リストの生成：要素をカンマ区切りで並べたリストをカギ括弧で括って表します。
空のリスト： []で空のリストを生成します。
要素の取得: リスト[要素のインデックス]
要素の追加： append関数を使用します。
指定範囲のサブリスト: リスト[開始インデックス:終了インデックス]
指定以降のサブリスト: リスト[開始インデックス:]
先頭から指定までのサブリスト: リスト[:終了インデックス]
リストの結合：＋でリストを連結します。

L1 = [1, 2, 3]; print L1 # リストの生成 L2 = [9, 8, 7] L1.append(4); print L1 # 要素の追加 print L1[2] # 要素の取得 print L1[1:3] # 部分リストの取得 print L1[2:] print L1[:2] print L1 + L2 # リストの結合

[1, 2, 3]
[1, 2, 3, 4]
3
[2, 3]
[3, 4]
[1, 2]
[1, 2, 3, 4, 9, 8, 7]

[1, 2, 3]
[1, 2, 3, 4]
3
[2, 3]
[3, 4]
[1, 2]
[1, 2, 3, 4, 9, 8, 7]

range関数

範囲生成関数rangeを使うと簡単に指定した範囲のインデックスのリストを生成することができます。

range(終了): ０から終了の前までの整数のリストを生成します。
range(開始, 終了): 開始インデックスから終了インデックスの前まで整数のリストを生成します。
range(開始, 終了, 増分): 開始から終了の前までの値を増分ごとに追加したリストを生成します。

print range(3) print range(1, 4) print range(-1, 7, 2)

[0, 1, 2]
[1, 2, 3]
[-1, 1, 3, 5]

[0, 1, 2]
[1, 2, 3]
[-1, 1, 3, 5]

その他の便利な関数

zip関数: 複数のリストの各要素をタプルとして結合したい場合に、zip関数を使います。
sum関数: リストの和を求めるときに、sum関数を使います。
prod関数: リストの積を求めるときに、prod関数を使います。

L1 = [1,2,3,4] t1 = ['a', 'b', 'c', 'd']; print L1, t1 print zip(L1, t1) # リストL1とリストt1の各要素をタプルとしたリストを作成します print sum(L1) # リストL1の値の和を求めます print prod(L1) # リストL1の値の積を求めます

[1, 2, 3, 4] ['a', 'b', 'c', 'd']
[(1, 'a'), (2, 'b'), (3, 'c'), (4, 'd')]
10
24

[1, 2, 3, 4] ['a', 'b', 'c', 'd']
[(1, 'a'), (2, 'b'), (3, 'c'), (4, 'd')]
10
24

リスト内包

リスト内包は、式とそれにつづくfor節からなり、各要素に対し、式を評価した結果をリストとして返します。最初は分かりにくいですがリストをコンパクトに生成できるので重宝します。

s1 = [sqrt(s) for s in L1]; print s1 # リストL1の各要素をsqrtを取った値のリストを返します print [n(s, digits=5) for s in s1] # リストs1の要素を5桁の数値としたリストを返します view(s1) # リストを整形して出力します

[1, sqrt(2), sqrt(3), 2]
[1.0000, 1.4142, 1.7320, 2.0000]
\newcommand{\Bold}[1]{\mathbf{#1}}\left[1, \sqrt{2}, \sqrt{3}, 2\right]

[1, sqrt(2), sqrt(3), 2]
[1.0000, 1.4142, 1.7320, 2.0000]
\newcommand{\Bold}[1]{\mathbf{#1}}\left[1, \sqrt{2}, \sqrt{3}, 2\right]

ディクショナリ

ディクショナリは、別名「連想配列」と言われ、リストのインデックスの代わりにキーで値を取得することができます。

ディクショナリの生成：キーと値のペアを:(コロン)ではさみ,(カンマ)で連結し｛｝カッコで括ります。
要素の取得；[]カッコにキーを指定して要素を取得します。

ages = {'John':24, 'Sarah':28, 'Mike':31}; # ディクショナリの生成 print ages # ディクショナリの出力 print ages['Sarah'] # 要素の取得

{'Sarah': 28, 'Mike': 31, 'John': 24}
28

{'Sarah': 28, 'Mike': 31, 'John': 24}
28

変数の宣言

sageの式で変数として認識させるには変数をvar関数で宣言しなくてなりません。変数の宣言は、

var('変数名')

複数の変数を宣言する場合には、スペースを空けて指定します。宣言される変数を参照したりする場合には、

x, y = var('x y')

のようにpython変数x, yに宣言したsage変数を代入します。

関数の定義

sageの関数には、以下の3種類の指定方法があります。必要に応じて使い分けて下さい。

シンボリック関数: f(x) = 式の形式で定義します。
python関数: pythonのdefを使って関数を定義します。
lambda関数: lambda 変数のリスト: 式の形式で定義します。

例では、$x^3$を返す関数をシンボリック関数、python関数、lambda関数をそれぞれ symCube, pyCube, lamCube として定義し、その使い方を説明します。当然すべて、同じ結果になります。

# シンボリック関数 symCube(x) = x^3 # python関数 def pyCube(x): return x^3 # lambda関数 lamCube = lambda x: x^3 # 関数の呼び出し print symCube(2), pyCube(2), lamCube(2) # 数式が引数の場合 x, y = var('x y') print symCube(x+y), pyCube(x+y), lamCube(x+y)

8 8 8
(x + y)^3 (x + y)^3 (x + y)^3

8 8 8
(x + y)^3 (x + y)^3 (x + y)^3

値の代入

変数の値を指定して関数の値を取得するには、以下の3つの方法があります。

値を指定したい変数を引数とするシンボリック関数を定義する
関数の引数にx=2のように変数に値を代入して参照する
subs関数を使って値を代入する

# 関数f、変数x, a, bを定義 var('x a b') f(x) = a*x + b g(x, a, b) = f print g(x, 2, 1) print f(a=2, b=1) print f.subs(a=2, b=1)

2*x + 1
2*x + 1
x |--> 2*x + 1

2*x + 1
2*x + 1
x |--> 2*x + 1

規則の代入

規則の代入には、subs_expr関数を使います。例として、関数hの$x^2$をwに置き換えてみます。

var('w') h = x^2 + 1; print h # 関数hの定義 print h.subs_expr(x^2 == w) # x^2をwに置き換える

x^2 + 1
w + 1

x^2 + 1
w + 1

前の出力結果の参照

直前の実行結果を _ で参照することができます。

[1, 2, 3]

[1, 2, 3]

[1, 2, 3]

sum(_) # 直前の実行結果を_で参照し、リストの和を求めます。

数の計算・基本的な関数

よく使われる定数

よく使われる定数を以下に示します注4。

hdr=["円周率", "自然体数の底", "虚数単位", "無限大"] sts=[pi, e, I, oo] html.table([(h, st, "%s"%st) for (h, st) in zip(hdr, sts)])





円周率
\pi
pi


自然体数の底
e
e


虚数単位
I
I


無限大
\infty
oo





円周率
\pi
pi


自然体数の底
e
e


虚数単位
I
I


無限大
\infty
oo

基本的な計算

sageで使われる基本的な計算の表現方法を以下に示します。

# よく使われる表現 var('a b x n') hdr=["積", "商", "累乗", "平方根", "n乗根", "絶対値", "自然対数", "階乗"] sts=[a*b, a/b, a^b, sqrt(x), x^(1/n), abs(x), log(x), factorial(n)] html.table([(h, st, "%s"%st) for (h, st) in zip(hdr, sts)])





積
a b
a*b


商
\frac{a}{b}
a/b


累乗
a^{b}
a^b


平方根
\sqrt{x}
sqrt(x)


n乗根
x^{\frac{1}{n}}
x^(1/n)


絶対値
{\left| x \right|}
abs(x)


自然対数
\ln\left(x\right)
log(x)


階乗
n!
factorial(n)





積
a b
a*b


商
\frac{a}{b}
a/b


累乗
a^{b}
a^b


平方根
\sqrt{x}
sqrt(x)


n乗根
x^{\frac{1}{n}}
x^(1/n)


絶対値
{\left| x \right|}
abs(x)


自然対数
\ln\left(x\right)
log(x)


階乗
n!
factorial(n)

式の計算・簡単化

数式処理システムの恩恵は、中学や高校の数学で悩まされた計算から解放？されることです。まずは、

展開: expand
因数分解: factor
簡単化: simplify

を試してみます。

f1 = (x - 1)*(x^2 - 1); print f1 # f1を定義 f2 = expand(f1); print f2 # f1を展開し、f2に代入 print factor(f2) # f2を因数分解する f3 = 1/(x+2)+1/(x-2); print factor(f3) # 分数式f3の因数分解 print simplify(I + x -x) # 簡単化で不要な項を消去

(x - 1)*(x^2 - 1)
x^3 - x^2 - x + 1
(x - 1)^2*(x + 1)
2*x/((x - 2)*(x + 2))
I

(x - 1)*(x^2 - 1)
x^3 - x^2 - x + 1
(x - 1)^2*(x + 1)
2*x/((x - 2)*(x + 2))
I

三角関数の簡単化

三角関数や指数関数を含む式の簡単化には simplify_full関数を使います注5。

f4 = cos(x)^2-sin(x)^2; print f4 # cos(x)^2 + sin(x)^2 = 1を使って f4.simplify_full() # 簡単化する

-sin(x)^2 + cos(x)^2
2*cos(x)^2 - 1

-sin(x)^2 + cos(x)^2
2*cos(x)^2 - 1

グラフ表示

数式処理システムのもう一つの便利な機能にグラフ表示があります。関数の形を可視化することによってその性質を把握することができます。

plot関数

plot関数の呼び出しは、以下のように指定します。

plot(関数, [変数名, 最小, 最大], オプション)

オプションは、省略可能です。plotのオプションは、plot.optionsで知ることができます。それ以外にもGraphicsのオプションも使えます。よく使うオプションは、

描画範囲指定のxmin, xmax, ymin, ymax
グラフの比率を指定するaspect_ratio
線の色指定のrgbcolor

があります。

plotの例として、$f(x) = \frac{sin(x)}{x}$を-100から100の範囲で表示してみます。

plot(sin(x)/x, [x, -100, 100])

plotはlimitを計算する

グラフをみて、おやっと思われた方もいると思いますが、$f(x) = \frac{x}{sin(x)}$ は、$x = 0$で分母が０になりエラーとなるはずですが、plotは、limitを計算して表示しています。これが、数式処理システムのグラフ化の強みです。

limit(sin(x)/x, x=0)

基本図形と重ね書き

計算結果の表示の他に、補足説明などのために基本図形を表示したい場合があります。以下によく使う基本図形（グラフオブジェクト）を示します。

円: circle
文字列: text注6
線: line
点: point
ポリゴン: polygon

重ね書きは、各グラフオブジェクトを+で結合し、show関数で表示します。

# 点: point((座標), オプション属性) 形式で指定 # 文字列: text(文字列, (座標)) 形式で指定（latex表現は$で囲む） # ポリゴン: polygon(座標のリスト) 形式で指定 c = circle((0.5,0.5), 1) # 円: circle((座標), 半径) 形式で指定 l = line([(0,0), (1, 1)]) # 線: line([(開始座標), (終了座標)]) 形式で指定 txt = text('$f(x) = x^2 + 1$', (1,1)) + text('テスト', (0,0)) # テストは文字化けします pt = point((0.5, 0.5), rgbcolor='white', pointsize=30, faceted=True) pol = polygon([(0,0), (-1,1), (0,1)], rgbcolor='yellow') # 重ね書き (c + l + pt + txt + pol).show(xmin=-1, xmax=2)

リストプロット

リストの値をプロットするには、list_plot関数を使用します。list_plotの使い方は以下の通りです。

list_plot(プロットするリスト, [オプションにはxmax, ymax, plotjoined等を指定])

# 0から24までの値の2乗の値をプロット list_plot([i^2 for i in range(25)], xmax=25)

3次元プロット

2次元グラフと同様に3次元のグラフを表示することができます。表示された図形はマウスで自由に回転することができます。（驚きました）

x, y = var('x y') plot3d(sin(x*y),(x,-pi,pi),(y,-pi,pi), mesh=True)

Sleeping...
If no image appears re-execute the cell. 3-D viewer has been updated.

その他の表示方法

その他の表示方法としては、

媒介変数表示: parametric_plot関数
関係式の表示: implicit_plot関数
等高線グラフ: contour_plot

があります。詳しくはヘルプとマニュアルを参照してください。

方程式

方程式の解法

数式処理システムのありがたい機能の一つに、方程式の解法があります。

方程式の解法は、solve関数を使って以下のように行います。

solve(解きたい方程式また方程式のリスト, 解を得る変数)

例として、$f(x)=ax+b$ の一次方程式を解いてみます。解として$x = \frac{−b}{a}$ が得られたので、$a=2, b=1$を代入して解の値を取得します。

# 方程式の解法 var('x a b') f = a*x + b sln = solve(f, x); view(sln) # a=2, b=1の場合の値を得るには print sln[0].subs(a=2, b=1)

\newcommand{\Bold}[1]{\mathbf{#1}}[
x == -b/a
]
x == (-1/2)

\newcommand{\Bold}[1]{\mathbf{#1}}[
x == -b/a
]
x == (-1/2)

多項式と数値解

多項式の解法を例に、解法、グラフ化、数値解の求め方を説明します。

多項式 $$ f(x) = x^3 - 2 x + 4 $$ をプロットし、方程式の解と数値解を比べてみます。数値解には、find_root関数を使用します。

var('x') f = x^3-2*x+4 # 多項式のグラフ化 plot(f, -2.5, 2.5)

# 解を求める sln = solve(f, x); view(sln) # グラフから-3から３の範囲で数値解を求める find_root(f, -3, 3)

\newcommand{\Bold}[1]{\mathbf{#1}}[
x == (-I + 1),
x == (I + 1),
x == -2
]
-1.9999999999999949

\newcommand{\Bold}[1]{\mathbf{#1}}[
x == (-I + 1),
x == (I + 1),
x == -2
]
-1.9999999999999949

連立方程式

方程式のリストをsolve関数に与えることで、連立方程式の解を得ることができます。

solve関数を使って、次の連立方程式を解いてみましょう。多項式と同様にグラフを使った解法も試してみます。 $$ \begin{eqnarray} x^2 - 2y & = & 2 \\ x - y & = & 1 \end{eqnarray} $$

var('x y') eq = [y == 1/2*x^2 - 1, y == x - 1] # 解を求める view(solve(eq, [x, y])) # グラフから解を求める p1 = plot(1/2*x^2 - 1, [x, -1, 3]) p2 = plot(x - 1, [x, -1, 3]) (p1+p2).show()

\newcommand{\Bold}[1]{\mathbf{#1}}[
[x == 2, y == 1],
[x == 0, y == -1]
]

\newcommand{\Bold}[1]{\mathbf{#1}}[
[x == 2, y == 1],
[x == 0, y == -1]
]

微分・積分

微分

微分は、diff関数を使って以下のように行います。

diff(関数, 微分する変数)

いろんな関数の微分を以下に示します。

var('x c n') sts=[c, x^n, sin(x), cos(x), exp(x), log(x)] html.table([("関数", "微分結果")] + [(st, diff(st, x)) for st in sts], header=true)





関数
微分結果


c
0


x^{n}
n x^{{(n - 1)}}


\sin\left(x\right)
\cos\left(x\right)


\cos\left(x\right)
-\sin\left(x\right)


e^{x}
e^{x}


\ln\left(x\right)
\frac{1}{x}





関数
微分結果


c
0


x^{n}
n x^{{(n - 1)}}


\sin\left(x\right)
\cos\left(x\right)


\cos\left(x\right)
-\sin\left(x\right)


e^{x}
e^{x}


\ln\left(x\right)
\frac{1}{x}

積分

積分には、integrate関数を使用します。sageでのintegrate関数の使い方は、以下の通りです。

integrate(被積分関数, 積分変数)
また
integrate(被積分関数, 積分変数, 積分範囲)
とすると定積分を計算します。

いろんな関数の積分を以下に示します注7。

var('x c') f = function('f', x) sts=[c, 1/x, exp(x), log(x), sin(x), diff(f, x)/f(x)] html.table([("関数", "積分結果")] + [(st, integrate(st, x)) for st in sts], header=true)





関数
積分結果


c
c x


\frac{1}{x}
\ln\left(x\right)


e^{x}
e^{x}


\ln\left(x\right)
x \ln\left(x\right) - x


\sin\left(x\right)
-\cos\left(x\right)


\frac{f'\left(x\right)}{f\left(x\right)}
\ln\left(f\left(x\right)\right)





関数
積分結果


c
c x


\frac{1}{x}
\ln\left(x\right)


e^{x}
e^{x}


\ln\left(x\right)
x \ln\left(x\right) - x


\sin\left(x\right)
-\cos\left(x\right)


\frac{f'\left(x\right)}{f\left(x\right)}
\ln\left(f\left(x\right)\right)

定積分

定積分の例として、$\int_0^{\pi/2}sin(x)dx$を計算してみます。

integrate(sin(x), x, 0, pi/2)

数値積分

数値積分には、numerical_integral関数を使用します。

numerical_integral(関数, 積分範囲)

数値積分の例として、サイクロイドを計算してみます。サイクロイドは、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(\frac{dy}{dt}\right)^2}dt $$ 解析解では、この結果は16となります。

t = var('t') x = 2*(t-sin(t)) y = 2*(1-cos(t)) # サイクロイドのプロット（パラメトリックプロットの例） parametric_plot([x, y], (t, 0, 2*pi))

# 数値積分の結果 cycloid = sqrt(diff(x,t)^2+diff(y,t)^2) numerical_integral(cycloid, 0, 2*pi)

(15.999999999999998, 1.7763568394002502e-13)

(15.999999999999998, 1.7763568394002502e-13)

ベクトルと行列

ベクトルの生成

ベクトルは、vector関数で生成します。

vector(値のリスト)
または、
vector(環, 値のリスト)

よく使われる環を以下に示します。

整数: ZZ
実数: R
有理数: QQ
複素数: CC
倍精度小数(real double): RDF

ベクトルの生成例として、$\mathbf{v}=(2,1,3),\mathbf{w}=(1,1,−4)$

v = vector([2,1,3]); print v w = vector([1,1,-4]); print w

(2, 1, 3)
(1, 1, -4)

(2, 1, 3)
(1, 1, -4)

ベクトルの内積

ベクトルの内積$$\mathbf{v}\cdot\mathbf{w}$ は、dot_product関数で計算します。

また、ベクトルの内積は、 $$ \mathbf{v}\cdot\mathbf{w} = |\mathbf{v}||\mathbf{w}|cos(\theta) $$ と表現され、2つのベクトルの角度$cos(\theta)$を知るために利用されます。

# ベクトルの内積 print v.dot_product(w) # 同様の計算は、ベクトルの積を使っても計算できます print v*transpose(w)

-9
(-9)

-9
(-9)

ベクトルの外積

ベクトルの外積$\mathbf{v}\times\mathbf{w}$は、cross_product関数で計算します。

ベクトルの外積の大きさは、 $$ |\mathbf{v}\times\mathbf{w}| = |\mathbf{v}||\mathbf{w}|sin(\theta) $$ であり、その方向は、ベクトルvからベクトルwにねじを回して進む方向となります。

# ベクトルの外積 print v.cross_product(w)

(-7, 11, 1)

(-7, 11, 1)

行列の生成

行列は、matrix関数で生成します。

matrix(行列の要素のリスト)
または
matrix(環, 行列の要素のリスト)

# 行列の生成 A = matrix([[1,2,3],[3,2,1],[1,2,1]]); A

[1 2 3]
[3 2 1]
[1 2 1]

[1 2 3]
[3 2 1]
[1 2 1]

特殊な行列（単位行列、対角行列）

特殊な行列として、

単位行列: identity_matrixで生成
対角行列: diagonal_matrixで生成

があります。

identity_matrixの呼び出し形式は以下の通りです。

identity_matrix(対角要素の数)

diagonal_matrixの呼び出し形式は以下の通りです。

diagonal_matrix(対角要素のリスト)

# 単位行列の生成 print "I="; print identity_matrix(3) # 対角行列 D = diagonal_matrix([1, 2, 3]); print "D="; print D

I=
[1 0 0]
[0 1 0]
[0 0 1]
D=
[1 0 0]
[0 2 0]
[0 0 3]

I=
[1 0 0]
[0 1 0]
[0 0 1]
D=
[1 0 0]
[0 2 0]
[0 0 3]

ベクトル・行列の基本演算

ベクトル、行列も基本演算子+, -, *を使って和、差、積を計算することができます。

# ベクトル、行列の積の例 print w*A print v+w

(0, -4, 0)
(3, 2, -1)

(0, -4, 0)
(3, 2, -1)

転置行列

転置行列は、要素の行と列を入れ替えた行列です。転置行列は、transpose関数で取得できます。

転置行列の性質には、

$ (\mathbf{A}^T)^T = \mathbf{A} $
$ (\mathbf{A+B})^T = \mathbf{A}^{T} + \mathbf{B}^{T}$
$ (\mathbf{A}\mathbf{B})^T = \mathbf{B}^T\mathbf{A}^T $

があり、最後の転置行列の積の公式には、OpenGLの座標変換処理でお世話になりました。

# 転置行列の例 print "Aの転置行列="; print A.transpose() # ベクトルの転置の例 print "vの転置行列="; print v.transpose()

Aの転置行列=
[1 3 1]
[2 2 2]
[3 1 1]
vの転置行列=
[2]
[1]
[3]

Aの転置行列=
[1 3 1]
[2 2 2]
[3 1 1]
vの転置行列=
[2]
[1]
[3]

逆行列

逆行列は、行列$\mathbf{A}$とその逆行列$\mathbf{A}^{-1}$を掛けると単位行列になる行列です。逆行列は、inverse関数で取得できます。

逆行列の性質には、

$ (\mathbf{A}^{-1})^{-1} = \mathbf{A} $
$ (\mathbf{A}^T)^{-1} = (\mathbf{A}^{-1})^T $
$ (\mathbf{A}\mathbf{B})^{-1} = \mathbf{B}^{-1}\mathbf{A}^{-1} $

があります。

# 逆行列の例 print A.inverse()

[   0  1/2 -1/2]
[-1/4 -1/4    1]
[ 1/2    0 -1/2]

[   0  1/2 -1/2]
[-1/4 -1/4    1]
[ 1/2    0 -1/2]

行列式

行列式は、det関数で取得できます。行列式は逆行列の計算に使われるので、逆行列が存在するか否かの判別に使用されます。

# 行列式の例 print A.det()

行列の解

solve_right関数を使って行列の解を得ることができます。 $$ \mathbf{Y} = \mathbf{A}\mathbf{X} $$ となるようなXを求めるのが、solve_right関数です。

# 行列の解 Y = vector([0,-4,-1]) X = A.solve_right(Y); print X # AXを計算して、Yと同じことを確認する print A*X

(-3/2, 0, 1/2)
(0, -4, -1)

(-3/2, 0, 1/2)
(0, -4, -1)

応用編

応用編では、

Rを使って入力データを処理（読み込む）する
find_fit関数を使って回帰分析を行う
回帰分析で求めた係数を使ってモデルの連立方程式を解く
求めた解から未来を予想する

方法を紹介します。

計量経済モデル

計量経済学では、経済の構造をモデルで表現し、モデルの要素間の相互関係を連立方程式で表し、過去の経済統計データからモデルの係数を求め、未来の予測を行います。

ここでは、「MathematicaによるOR」4章で紹介されているモデルと経済データを使って、sageの活用例を紹介します。

以下のモデルでは、矩形を要素としY:所得、C:消費、I:投資、 G:政府支出とし、矢印の方向は影響与える向きを表します。消費は自分自身へ矢印があり、これは前年度の消費が今年度の消費に影響を与えていることを表しています。

このモデルを連立方程式で表すと、以下のようになります。 $$ \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} $$

データの読み込み

sageでは、統計情報の分析に有益なRパッケージとのインタフェースを持っています。ここでは、Rのデータ取り込みを関数であるread.table関数を使ってデータを読み込み、取り込んだデータをsageの形式に変換する方法を紹介します。

Rとのインタフェースは、以下のようになっています。

r("Rのコマンド")

Rでデータを読み込むコマンドは、read.table関数でオプションにheader=Tを付けて最初の行がヘッダであることを知らせます。

read.table('入力ファイル名' [,オプション] )

ファイル名に'%s'とし、%filenameをコマンド文字列の後に付けることで %sが変数filenameの値に置き換えれます。

経済データは、Econometrics.txtファイルに入っています注8。最初の1行がYear, Ct, Yt, It, Gtのヘッダで、その後に 1965年から1989年までの各年度の消費、所得、投資、政府支出の値が続きます。

# 計量経済学(Econometrics) # Rのread.table関数を使ってデータを読み込む fileName = DATA + 'Econometrics.txt' d = r("read.table('%s', header=T)" %fileName) ; d

   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

   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_() 関数を呼び出します。 Rインタフェースの_sage_関数の仕様で、変換されたデータdsの 'DATA'要素にヘッダをキーとした列データのディクショナリがセットされます。 'DATA'要素からCt, Yt, It, Gtをキーとして消費、所得、投資、政府支出の列データを取り出します。

# 読み込んだデータ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']

係数の決定

連立方程式のCt, Itの係数$a_0, a_1, a_2, b_0, b_1$の値を回帰分析を使って求めます。

CtList, YtList, GtListには、1966年から1984年のデータをセットし、 Ct1Listには、1年ずらしたCtの値をセットします。

# Ctの計算には前年度のCtであるCt-1を使うため1966年から1984年のデータを作成 CtList = CtData[1:20] Ct1List= CtData[0:19] YtList = YtData[1:20] GtList = GtData[1:20] ItList = ItData[1:20]

回帰分析関数find_fitの使い方

回帰分析には関数find_fitを使います。 find_fit関数の使い方は、以下の通りです。

find_fit(data, model, options)
data: 変数1の値, 変数2の値, ... , 変数nの値, 関数の値をタプルとするリストをセット
model: model(変数1の値, 変数2の値, ... , 変数nの値)の引数を持つ関数をセット
options: solution_dict=Trueを指定すると係数名をキーとする辞書が戻されます

zip関数を使って各年度のYt, Ct1, Ctをタプルにまとめたリストを作成します。

# 各年度のYt, Ct1, Ctをタプルにまとめたリストを作成 data = zip(YtList, Ct1List, CtList)

消費の係数を求める

変数Yt, Ct1, a0, a1, a2とモデルCtModelを定義します。CtModelの式は、連立方程式の$C_t = a_0 + a_1 Y_t + a_2 C_{t-1}$の部分です。

a0: 7518.3096037375299
a1: 0.21857976971146531
a2: 0.59268191617908206

が求まりました。

# モデル関数と変数を定義 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

{a0: 7518.3096037375299, a1: 0.21857976971146531, a2:
0.59268191617908206}

{a0: 7518.3096037375299, a1: 0.21857976971146531, a2: 0.59268191617908206}

投資の係数を求める

消費と同様に、変数Yt, b0, b1とモデルItModelを定義します。 CtModelの式は、連立方程式の$I_t = b_0 + b_1 Y_t$の部分です。

b0: -11736.272856134668
b1: 0.22718320402830394

が求まりました。

# 同様に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

{b0: -11736.272856134668, b1: 0.22718320402830394}

{b0: -11736.272856134668, b1: 0.22718320402830394}

連立方程式を解く

係数a0, a1, b2, b0, b1が決まったので、solve関数を使って連立方程式を解きます。

係数を入力して式を定義する代わりに、

CtEq = (Ct == CtModel(Yt, Ct1)).subs(CtFit)

のようにsubs関数を使ってCt == a0 + a1*Yt + a2*Ct1のa0, a1, a2を代入しています。これで、モデルを変更しても再計算が楽にできます。解の結果、Ct, It, YtはGtとCt1からなる式に変換されています。

# 連立方程式を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

[Ct == 0.592681916179082*Ct1 + 0.218579769711465*Yt + 7518.30960373753,
It == 0.227183204028304*Yt - 11736.2728561347, Yt == Ct + Gt + It]

[Ct == 0.592681916179082*Ct1 + 0.218579769711465*Yt + 7518.30960373753, It == 0.227183204028304*Yt - 11736.2728561347, Yt == Ct + Gt + It]

sol = solve(eqn, [Ct, It, Yt], solution_dict=True); sol

[{It: 163953670168/674867923015*Ct1 + 115022920/280610363*Gt -
872829386902217/64820993853, Yt: 721680419112/674867923015*Ct1 +
506300280/280610363*Gt - 164437808222620/21606997951, Ct:
557726748944/674867923015*Ct1 + 110666997/280610363*Gt +
379515962234357/64820993853}]

[{It: 163953670168/674867923015*Ct1 + 115022920/280610363*Gt - 872829386902217/64820993853, Yt: 721680419112/674867923015*Ct1 + 506300280/280610363*Gt - 164437808222620/21606997951, Ct: 557726748944/674867923015*Ct1 + 110666997/280610363*Gt + 379515962234357/64820993853}]

計算結果と実データの比較

zip関数を使って前年度消費と政府支出をタプルとするリストを作成します。

# 入力データをセット data = zip(Ct1List, GtList)

消費データの比較

連立方程式からCt式を取り出し、前年度消費Ct1と政府支出Gtを引数とするctFuncを定義します。

ctFunc(Ct1, Gt) = sol[0][Ct]; ctFunc

(Ct1, Gt) |--> 557726748944/674867923015*Ct1 + 110666997/280610363*Gt
+ 379515962234357/64820993853

(Ct1, Gt) |--> 557726748944/674867923015*Ct1 + 110666997/280610363*Gt + 379515962234357/64820993853

消費データ計算

dataの各タプルに対し、リスト内包を使ってctFunc関数から消費の計算結果のリストを作成し、calCtListにセットします。

calCtList = [ctFunc(Ct1, Gt).n() for (Ct1, Gt) in data]

消費データのプロット

list_plotを使って計算結果（青い線）と実データ（赤い点）をプロットします。単純なモデルの割には計算値と実データが合っています。

# 計算値（青）と実測値（赤） calCtPlot = list_plot(calCtList, plotjoined=true) ctPlot = list_plot(CtList, rgbcolor='red') (calCtPlot + ctPlot).show(ymin=0)

投資データと所得データの計算

同様に投資データ、所得データを計算します。

# 同様にItに対する計算値と実測値をプロット itFunc(Ct1, Gt) = sol[0][It] calItList = [itFunc(Ct1, Gt) for (Ct1, Gt) in data] # 同様にYtに対する計算値と実測値をプロット ytFunc(Ct1, Gt) = sol[0][Yt] calYtList = [ytFunc(Ct1, Gt) for (Ct1, Gt) in data]

未来の予測

計量経済モデルの目標である、未来の予測をしてみましょう。

問題を簡単にするために、政府支出Gtを$G_t(x) = c_0 + c_1 x$として回帰分析をします。

# 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

{c0: 35741.149903570993, c1: 1960.3166656122035}

{c0: 35741.149903570993, c1: 1960.3166656122035}

政府支出の予測結果

直線回帰で求めた政府支出と実データをプロットすると以下のようになります。

gtPlot = list_plot(GtList, rgbcolor='red') gtFunc(x) = GtModel.subs(GtFit) calGtPlot = plot(gtFunc, [x, 0, 25]) (calGtPlot + gtPlot).show()

1985年から1989年までの予測

政府支出関数gtFuncを使って、1985年から1989年まで政府支出を求め、前年度の消費Ct1List の値とgtFuncの値から消費、投資、所得を順番に計算します。

消費の予測

予測した消費を1985年から1989年までの実データと重ねてプロットします。ラフな政府支出のモデルを使った割にはうまく予測できています。

# 予測 predCt1List = Ct1List predCtList = [calCtList[18]] predItList = [calItList[18]] predYtList = [calYtList[18]] for i in range(19, 24): predCt1List.append(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))) # 計算値（青）、予測値（緑）と実測値（赤） calcCtPlot = list_plot(calCtList, plotjoined=true) predCtPlot = list_plot(zip(range(18,24), predCtList), plotjoined=true, rgbcolor ='green', linestyle='dashed') ctPlot = list_plot(CtData[1:25], rgbcolor='red') (calcCtPlot + predCtPlot + ctPlot).show(ymin=0)

投資と所得の予測結果

同様に投資と所得の予測と実データを重ねてプロットします。

# 同様にItに対する予測値と実測値をプロット calcItPlot = list_plot(calItList, plotjoined=true) predItPlot = list_plot(zip(range(18,24), predItList), plotjoined=true, rgbcolor ='green', linestyle='dashed') itPlot = list_plot(ItData[1:25], rgbcolor='red') (calcItPlot + predItPlot + itPlot).show(ymin=0) # 同様にYtに対する予測値と実測値をプロット calcYtPlot = list_plot(calYtList, plotjoined=true) predYtPlot = list_plot(zip(range(18,24), predYtList), plotjoined=true, rgbcolor ='green', linestyle='dashed') ytPlot = list_plot(YtData[1:25], rgbcolor='red') (calcYtPlot + predYtPlot + ytPlot).show(ymin=0)

まとめ

このようにsageを使うと簡単に計量経済モデルを計算をすることができます。モデルやデータの変更も容易なので、モデルの改良も短時間にできます。

また、sageのノートブックは計算だけではなく、手順などの説明も付け加えることができるので、動く解説書のような役割もします。

本稿を執筆するにあたり、永山純一さん、戸張一夫さん、鈴木由宇さんには原稿に対し貴重なコメントをいただきました。この場を借りて感謝申し上げます。

脚注

William Steinの2007年UW CSE コロキアムの資料(http://boxen.math.washington.edu/talks/2007-01-15-sage-cse/sage.pdf)から引用。

William Steinの2007年UW CSE コロキアムの資料にRを追加、sageに含まれるオープンソースのコンポーネントは、 http://www.sagemath.org/links-components.html で参照できます。

http://wiki.sagemath.org/quickref?action=AttachFile&do=view&target=quickref.pdf で参照できます。

無限大の出力は、+∞と+Infinityと出力されるため、∞とooに修正しました。

変数の後に.(ピリオド)を付けて関数（メソッド）を呼び出すことができます。

残念ながら日本語の文字は化けてしまいします。

diff(f, x)の出力がD[0](f)(x)となるため、f'に修正しました。

ワークシートのトップにあるFile...メニューからUpload worksheet from fileを選択し、Econometrics.txtをアップロードしておきます。

Sage_introduction

4831 days ago by takepwave

はじめに

sageとは

sageの誕生

sageに含まれるオープンソースのコンポーネント

基礎編

インストールなしで使えるsage

ノートブックの使い方

ワークシートの操作

セルの操作

ヘルプ

補完機能

表記方法

python文法

コメント

複数の式を1行にまとめる

python変数への代入

タプル

リスト

range関数

その他の便利な関数

リスト内包

ディクショナリ

変数の宣言

関数の定義

値の代入

規則の代入

前の出力結果の参照

数の計算・基本的な関数

よく使われる定数

基本的な計算

式の計算・簡単化

三角関数の簡単化

グラフ表示

plot関数

plotはlimitを計算する

基本図形と重ね書き

リストプロット

3次元プロット

その他の表示方法

方程式

方程式の解法

多項式と数値解

連立方程式

微分・積分

微分

積分

定積分

数値積分

ベクトルと行列

ベクトルの生成

ベクトルの内積

ベクトルの外積

行列の生成

特殊な行列（単位行列、対角行列）

ベクトル・行列の基本演算

転置行列

逆行列

行列式

行列の解

応用編

計量経済モデル

データの読み込み

係数の決定

回帰分析関数find_fitの使い方

消費の係数を求める

投資の係数を求める

連立方程式を解く

計算結果と実データの比較

消費データの比較

消費データ計算

消費データのプロット

投資データと所得データの計算

未来の予測

政府支出の予測結果

1985年から1989年までの予測

消費の予測

投資と所得の予測結果

まとめ