NumPy
:高速化#
# 警告メッセージを非表示
import warnings
warnings.filterwarnings("ignore")
NumPy
は,数値計算をする上で重要な役割を果たすパッケージであり,特に,行列計算に威力を発揮する。NumPy
は「ナンパイ」と読む。慣例としてnp
としてインポートする。
import numpy as np
array
#
基本となる関数がnp.array()
であり,次のコードでは1次元配列を作る。
arr = np.array([10, 20, 30, 40, 50])
arr
array([10, 20, 30, 40, 50])
次に2次元配列を作成しよう。
mat = np.array([[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12], [13, 14, 15, 16]])
mat
array([[ 1, 2, 3, 4],
[ 5, 6, 7, 8],
[ 9, 10, 11, 12],
[13, 14, 15, 16]])
mat
は行列と解釈できる。データ型を確認してみよう。
print(type(arr))
print(type(mat))
<class 'numpy.ndarray'>
<class 'numpy.ndarray'>
arr
とmat
はNumPy
のndarray
(n
次元array
)というデータ型(クラス)である。[]
に囲まれた数字が並んでいてリストと似ているように感じるかも知れないが、処理速度が格段に早く、より簡単なコードで同じ処理ができるなど使い勝手が大きく向上することがNumPy
の大きな特徴となる。科学計算やデータ処理では欠かせないパッケージとなっている。
次のコードが示すように、リストからndarray
を作成することができる
l = [1,2,3,4,5]
np.array(l)
array([1, 2, 3, 4, 5])
np.array([l,l])
array([[1, 2, 3, 4, 5],
[1, 2, 3, 4, 5]])
後述するが、ndarray
をリストに変換するメソッドも用意されている。
array
についてこのサイトでは図を使って説明しているので興味があるひとは参考にしてはどうだろうか。
1次元配列#
要素の抽出#
リストと同じようにarr
の要素を抽出するには要素のインデックスを使う。
arr[1]
20
複数の要素を抽出するにはインデックスをリストとして使う。
arr[[0,1,3]]
array([10, 20, 40])
要素を連続で抽出するスライシング(slicing)も使うことができる。
arr[1:3+1]
array([20, 30, 40])
リストの場合と同じ様に,:
の左側が選択する最初の要素であり,:
の右側が選択する最後の次の番号である(即ち,:
の右側の番号の要素は含まれない。
要素の代入と削除#
要素の入れ替え方法を説明するために,上で作成したリストarr
を考えてみよう。
arr
array([10, 20, 30, 40, 50])
30
を3000
に入れ替えてみよう。リストや辞書と同じように=
を使い値を代入する。
arr[2] = 3000
arr
array([ 10, 20, 3000, 40, 50])
要素を削除するには関数np.delete()
を使う。第一引数にarray
を,第二引数に削除したい要素のインデックスを指定する。例えば,3000
を削除したいとしよう。
np.delete(arr,2)
array([10, 20, 40, 50])
ここでもう一度arr
を表示してみよう。
arr
array([ 10, 20, 3000, 40, 50])
3000
は残ったままである。np.delete()
を実行すると,もし3000
を削除するとどうなるかを示している。arr
自体を変更したい場合は=
を使って再割り当てをする必要がある。
arr = np.delete(arr,2)
arr
array([10, 20, 40, 50])
これはnp.delete()
特有なものではなく,他の関数(例えば,np.log()
やnp.sqrt()
)も同じである。
2次元配列:要素の抽出#
第\(i\)行目の第\(j\)列目要素にアクセスするためにはmat[i,j]
と書く。[,]
の,
を挟んで左が行を表し,右が列を示す。
[行のインデックス,列のインデックス]
コードはキーストロークが少なく,簡単なものが良いと考えられている。一方で,Python
は初心者に比較的簡単な言語だと言われるが,それでも関数・メソッドの数は膨大で,そのオプションとの組み合わせを考えるとまさしく「無数」にあると言っても過言ではない。そのため初心者にとって関数の使い方やオプションの書き方を間違う可能性は小さくない。さらに,自分が書いたコードを数週間・数ヶ月後に読み直すと,何をしようとしているのか分からないという状況が生じることもある。従って,以下の点が非常に重要になる。
間違いにくいコードの書き方を覚える。
高い可読性を意識したコードの書き方を覚える。
このスタンスに基づいて以下のルールに従って説明する。
[,]
内の,
を省略可能な場合であっても省略しない。[,]
内の,
の左または右を省略できる場合であっても省略しない。
行または列を連続選択する(slicing)場合を考えよう。上でも説明したように:
を使うが
start:end
となる。ここでstart
とは選択する要素の最初インデックスであり,end
は選択する最後の要素の次のインデックスである(リストの場合と同じ)。このルールに従うと,
:
の左を省略すると「最初から」という意味になる。:
の右を省略すると「最後まで」という意味になる。,
の左または右が:
のみの場合,「全て」という意味になる。
これを読むだけでは分かりにくいと思うので,以下の例に目を通してもう一度この箇所の説明を読んむことを推奨する。
mat[0,1]
2
第\(i\)行の抽出はmat[i,:]
で可能である。:
は「列の全て」という意味になる。
r0 = mat[0,:]
r0
array([1, 2, 3, 4])
抽出した行は上で説明した1次元配列なのでインデックスやスライシングを使って要素にアクセスすることができる。
複数の行の抽出は次の方法で可能である。
mat[1:3,:]
array([[ 5, 6, 7, 8],
[ 9, 10, 11, 12]])
第3行目は含まれないことに注意しよう。リストの要素の取り出し方と同じように:
の右のインデックスの行は含まれない。
(注意)
,:
を省略してmat[1:3]
と書いてもエラーは発生せず同じ結果が返されるが,,:
があることにより,行を抽出しており列は全て選択されていることが明示的になる。
第\(i\) 列目を抽出したい場合は次のコードになる。
mat[:, 1]
array([ 2, 6, 10, 14])
複数列の抽出は以下のようにする。
mat[:, 1:3]
array([[ 2, 3],
[ 6, 7],
[10, 11],
[14, 15]])
array
を使った計算#
まず2つの行列(2次元配列array
)を作成する。
m1 = np.array([[1, 1], [1, 1]])
m1
array([[1, 1],
[1, 1]])
m2 = np.array([[2, 2], [2, 2]])
m2
array([[2, 2],
[2, 2]])
これらは数学で習った行列のように見えるが,少し異なる性質(特に積)を持っているので注意しよう。
加算#
要素ごとの和となる。
m2 + m1
array([[3, 3],
[3, 3]])
減算#
要素ごとの差となる。
m2 - m1
array([[1, 1],
[1, 1]])
スカラー乗算#
与えられた定数とそれぞれの要素の積となる。
10 * m1
array([[10, 10],
[10, 10]])
m1/2 # 1/2をかけるのと同じ
array([[0.5, 0.5],
[0.5, 0.5]])
乗算:バージョン1#
*
を使うと要素どうしの積となる。数学で学んだ行列の計算とは異なるので注意すること。
m1*m2
array([[2, 2],
[2, 2]])
乗算:バージョン2#
@
を使うと数学で学んだ行列の積となる。
m1@m2
array([[4, 4],
[4, 4]])
numpy
の関数dot()
は@
と同じとなる。
np.dot(m1,m2)
array([[4, 4],
[4, 4]])
転置#
数学で学ぶ転置行列と同じ。m3
を使って説明する。
m3 = np.array([[1,2,3],[4,5,6]])
m3
array([[1, 2, 3],
[4, 5, 6]])
m3
のメソッドtranspose()
を使う。
m3.transpose()
array([[1, 4],
[2, 5],
[3, 6]])
.transpose()
の省略形として.T
を使うこともできる。
m3.T
array([[1, 4],
[2, 5],
[3, 6]])
linalg
サブパッケージ#
2次元配列のarray
は行列として使える。その特徴を活かしてNumPy
のサブパッケージであるlinalg
(linear algebraの略であり「線形代数」という意味)を使うと連立方程式を解いたり固有値を簡単に計算することができる。
連立一次方程式の解#
連立一次方程式は行列表記できる。例えば,
ここで
解は逆行列を使い計算できる。
solve
#
z
を求めるためにNumPy
のサブパッケージlinalg
の中にあるsolve
関数を使うが,次のような書き方となる。
np.linalg.solve(a,b)
このコードは.
を助詞「の」として,()
を「実行する」と読むと,次のように日本語として読むことができる。「NumPy
のlinalg
のsolve
関数を引数a
とb
を使い実行する」と読める。では実際にsolve
を使って解を求めよう。
a = np.array([[1,2],
[3,4]])
b = np.array([20,10])
sol = np.linalg.solve(a,b)
ここでsol
は解からなるarray
である。
sol
array([-30., 25.])
分かりやすく次のように表示しよう。
print(f'y = {sol[0]}\nx = {sol[1]}')
y = -30.0
x = 25.0
inv
#
array
の逆行列は,サブパッケージlinalg
の中にあるinv
で計算できる。a
を例にとると,書き方は次のようになる
np.linalg.inv(a)
solve
と同じように日本語として読むことができる。
ではa
の逆行列を計算して解を求めてみよう。
a = np.array([[1,2],[3,4]])
a_inv = np.linalg.inv(a)
b = np.array([20,10])
a_inv@b
array([-30., 25.])
当たり前だが,同じ解である。
固有値:eigvals
#
行列の固有値を計算するeigvals()
関数を紹介する。a
を例に使うと次の書き方となる。
np.linalg.eigvals(a)
この関数は動学モデルの安定性(収束、発散)を調べる場合に役に立つ。より具体的な説明は後程おこなうことにする。
a = np.array([[1,2],[3,4]])
np.linalg.eigvals(a)
array([-0.37228132, 5.37228132])
NumPy使用時によく使う属性とメソッド#
dir()
を使ってm3
の属性とメソッドを調べてみる。
dir(m3)
['T',
'__abs__',
'__add__',
'__and__',
'__array__',
'__array_finalize__',
'__array_function__',
'__array_interface__',
'__array_prepare__',
'__array_priority__',
'__array_struct__',
'__array_ufunc__',
'__array_wrap__',
'__bool__',
'__class__',
'__class_getitem__',
'__complex__',
'__contains__',
'__copy__',
'__deepcopy__',
'__delattr__',
'__delitem__',
'__dir__',
'__divmod__',
'__dlpack__',
'__dlpack_device__',
'__doc__',
'__eq__',
'__float__',
'__floordiv__',
'__format__',
'__ge__',
'__getattribute__',
'__getitem__',
'__gt__',
'__hash__',
'__iadd__',
'__iand__',
'__ifloordiv__',
'__ilshift__',
'__imatmul__',
'__imod__',
'__imul__',
'__index__',
'__init__',
'__init_subclass__',
'__int__',
'__invert__',
'__ior__',
'__ipow__',
'__irshift__',
'__isub__',
'__iter__',
'__itruediv__',
'__ixor__',
'__le__',
'__len__',
'__lshift__',
'__lt__',
'__matmul__',
'__mod__',
'__mul__',
'__ne__',
'__neg__',
'__new__',
'__or__',
'__pos__',
'__pow__',
'__radd__',
'__rand__',
'__rdivmod__',
'__reduce__',
'__reduce_ex__',
'__repr__',
'__rfloordiv__',
'__rlshift__',
'__rmatmul__',
'__rmod__',
'__rmul__',
'__ror__',
'__rpow__',
'__rrshift__',
'__rshift__',
'__rsub__',
'__rtruediv__',
'__rxor__',
'__setattr__',
'__setitem__',
'__setstate__',
'__sizeof__',
'__str__',
'__sub__',
'__subclasshook__',
'__truediv__',
'__xor__',
'all',
'any',
'argmax',
'argmin',
'argpartition',
'argsort',
'astype',
'base',
'byteswap',
'choose',
'clip',
'compress',
'conj',
'conjugate',
'copy',
'ctypes',
'cumprod',
'cumsum',
'data',
'diagonal',
'dot',
'dtype',
'dump',
'dumps',
'fill',
'flags',
'flat',
'flatten',
'getfield',
'imag',
'item',
'itemset',
'itemsize',
'max',
'mean',
'min',
'nbytes',
'ndim',
'newbyteorder',
'nonzero',
'partition',
'prod',
'ptp',
'put',
'ravel',
'real',
'repeat',
'reshape',
'resize',
'round',
'searchsorted',
'setfield',
'setflags',
'shape',
'size',
'sort',
'squeeze',
'std',
'strides',
'sum',
'swapaxes',
'take',
'tobytes',
'tofile',
'tolist',
'tostring',
'trace',
'transpose',
'var',
'view']
この中にtranspose()
やT
があるのが確認できる。この中でよく使うのがshape
であり,行と列の数を確認する場合に有用なデータ属性である。
m3.shape
(2, 3)
2x3のarray
であることが確認できる。
次に,抽出した行または列をリストに変換するメソッドについて説明する。
r0 = mat[0,:]
r0
array([1, 2, 3, 4])
これをリストに変換するためにtolist()
というメソッドを使う。
r0_list = r0.tolist()
r0_list
[1, 2, 3, 4]
type(r0_list)
list
よく使うNumPy関数#
浮動小数点の比較#
次のコードを考えよう。
0.3 == 0.1 + 0.2
False
これはバグ(bug)ではない。浮動小数点は近似値として扱われるために発生する(詳細は「浮動小数点型について」を参照)。従って,浮動小数点が等しいかを確認するために==
を使うのは避けるべきである。その代わりにNumPy
の関数.isclose()
を使うようにしよう。使い方は比べる値を引数に指定するだけで良い。
np.isclose(0.3, 0.1+0.2)
True
様々なNumPy
関数#
ルート \(\left(\sqrt x\right)\)
np.sqrt(4)
2.0
底が\(e\)の指数関数(\(e^x\))
np.exp(10)
22026.465794806718
自然対数(\(\log_ex\)または\(\ln x\))
np.log(10)
2.302585092994046
0
が\(N\)個のarray
を作る。
np.zeros(N)
np.zeros(10)
array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0.])
1.0
が\(N\)個のarray
を作る。
np.ones(N)
np.ones(10)
array([1., 1., 1., 1., 1., 1., 1., 1., 1., 1.])
\(a\)から\(b\)までの区間を等間隔に割った\(N\)個の数字を返す。
np.linspace(a,b,N)
np.linspace(0,1,5)
array([0. , 0.25, 0.5 , 0.75, 1. ])
\(a\)から\(b\)までの区間で\(m\)ステップずつ増加し等間隔に割った数字を返す(\(b\)は含まない)。
np.arange(a,b,m)
m = 1
の場合,組み込み関数のrange(a,b)
と同じ数字を生成するが,返り値がarray
であることが異なる。
np.arange(5,10,0.5)
array([5. , 5.5, 6. , 6.5, 7. , 7.5, 8. , 8.5, 9. , 9.5])
標本平均
x
が数字のarray
やリストの場合
np.mean(x)
\(=\bar{x}=\frac{1}{n}\sum_{i=1}^{n}x_i\)
xx = [1,2,3,4,5,6]
np.mean(xx)
3.5
標本中央値
x
が数字のarray
やリストの場合
np.median(x)
np.median(xx)
3.5
標本分散
x
が数字のarray
やリストの場合
np.var(x, ddof=0)
\(=s_x^2=\dfrac{1}{1-\text{ddof}}\sum_{i=1}^n\left(x_i-\bar{x}\right)^2\)(ddof=0
がデフォルト)
(注意)計量経済学で習う分散の不偏推定量はddof=1
が必要!
np.var(xx,ddof=1)
3.5
標本標準偏差
x
が数字のarray
やリストの場合
np.std(x, ddof=0)
\(=s_x=\sqrt{s_x^2}\) (ddof=0
がデフォルト)
(注意)標本標準偏差の場合,必ずしもddof=1
により不偏推定量とはならないが,通常ddof=1
を用いる。
np.std(xx,ddof=1)
1.8708286933869707
標本共分散
2次元以上のarray
やリストの場合
np.cov(xy, ddof=0)
\(=c_{xy}=\dfrac{1}{1-\text{ddof}}\sum_{i=1}^n(x_i-\bar{x})(y_i-\bar{y})\)(ddof=0
がデフォルト)
(注意)計量経済学で習う分散の不偏推定量はddof=1
が必要!
下の計算結果
\(c_{xy}=-0.6\)
\(s_x^2=3.5\) ([1,2,3,4,5,6]の分散)
\(s_y^2=4.4\) ([1,6,2,5,3,1]の分散)
xy = [[1,2,3,4,5,6],[1,6,2,5,3,1]]
np.cov(xy,ddof=1)
array([[ 3.5, -0.6],
[-0.6, 4.4]])
標本相関係数
2次元以上のarray
やリストの場合
np.corrcoef(xy)
\(=r_{xy}=\dfrac{c_{xy}}{s_x\cdot s_y}\)
(注意)ddof
の影響はない。
下の計算結果
\(r_{xy}=-0.152...\)
\(r_{xx}=r_{yy}=1\)
np.corrcoef(xy)
array([[ 1. , -0.15289416],
[-0.15289416, 1. ]])
効用関数#
関数にif
文を組み込む例を考える。相対的リスク回避度一定効用関数は次のように定義される。
def
を使って表すと次のようになる。
def utility(x,s):
if s != 1:
return (x**(1-s)-1) / (1-s)
else:
return np.log(x)
utility(10, 1)
2.302585092994046
utility(10, 0.2)
6.636966806002416
array
vs list
#
ここではlist
とarray
の重要な違いについて説明する。
次のリストのそれぞれの要素に10
を足したいとしよう。
list0 = [1.0, 2.0, 3.0, 4.0, 5.0]
for
ループを使うと次のようになる。
list1 = []
for i in list0:
list1.append(i + 10)
list1
[11.0, 12.0, 13.0, 14.0, 15.0]
複雑さが残る。また次のコードでは10
を最後に追加し,リスト同士を結合するだけである。
list0 + [10]
[1.0, 2.0, 3.0, 4.0, 5.0, 10]
より簡単なコードで実行できれば良いが,以下のコードではエラーが発生する。
list0 + 10
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
Cell In[65], line 1
----> 1 list0 + 10
TypeError: can only concatenate list (not "int") to list
これを実現するのがNumPy
のarray
である。
まずarray
を作成する。
arr0 = np.array(list0)
arr0
array([1., 2., 3., 4., 5.])
arr0 + 10
array([11., 12., 13., 14., 15.])
この機能はブロードキャスティング(Broadcasting)もしくはベクトル演算(Vectorization)と呼ばれ、ループを使わずに個々の要素に直接働きかけ計算している。この機能により、より高速な計算が可能となるばかりか、より短いコードでそれを実現できる。+
, -
, *
, **
や他の関数にも同様に使うことができる。以下で例を挙げる。
arr0 - 5
array([-4., -3., -2., -1., 0.])
arr0 * 10
array([10., 20., 30., 40., 50.])
arr0 ** 2
array([ 1., 4., 9., 16., 25.])
np.sqrt(arr0)
array([1. , 1.41421356, 1.73205081, 2. , 2.23606798])
np.log(arr0)
array([0. , 0.69314718, 1.09861229, 1.38629436, 1.60943791])
次の計算も可能である。
y = arr0 * 2 + np.sqrt(arr0) + 10
y
array([13. , 15.41421356, 17.73205081, 20. , 22.23606798])
この機能はNumPy
の行列でも有効である。
mat0 = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
mat0
array([[1, 2, 3],
[4, 5, 6],
[7, 8, 9]])
mat0 * 10
array([[10, 20, 30],
[40, 50, 60],
[70, 80, 90]])
np.sqrt(mat0)
array([[1. , 1.41421356, 1.73205081],
[2. , 2.23606798, 2.44948974],
[2.64575131, 2.82842712, 3. ]])
np.log(mat0)
array([[0. , 0.69314718, 1.09861229],
[1.38629436, 1.60943791, 1.79175947],
[1.94591015, 2.07944154, 2.19722458]])
マクロ経済学を考える#
ケインズの45度線モデル#
所得:\(Y\)
計画支出:\(E=C+I+G\)
消費関数:\(C=d+cY\)
自発的消費:\(d=10\)
限界消費性向:\(c=0.7\)
投資:\(I=15\)
政府支出:\(G=10\)
均衡:\(E=Y\)
2つの式で表す。
i = 15 # I
g = 10 # G
d = 10 # d
c = 0.7 # c
a = np.array([[1,-c],[1,-1]])
b = np.array([i+g-d, 0])
sol = np.linalg.solve(a,b)
print(f'E = {sol[0]:.1f}\nY = {sol[1]:.1f}')
E = 50.0
Y = 50.0
IS曲線の導出:2つの方程式を使う#
for
ループを使う:方法1#
上の45度線モデルの投資を次のように変更する。
投資関数:\(I=f-hr\)
\(r\):実質利子率
\(f=20\):\(r=0\)の場合の投資
\(h=1.0\):実質利子率が1%上昇した場合,減少する投資量
実質利子率は次の値を取るとする。
r_list = range(0, 10)
list(r_list)
[0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
それぞれの実質利子率の値での均衡での所得と投資を計算する。まず2つの式で表す。
f = 20.0
h = 1.0
y_list = []
inv_list = []
for r in r_list:
# 異なる実質利子率の値で解を計算する
a = np.array([[1-c, -1],[0,1]])
b = np.array([g-d,f-h*r])
sol = np.linalg.solve(a,b)
y_list.append(sol[0].round(1)) # round(x)は小数点第x位まで四捨五入するメソッド
inv_list.append(sol[1])
print('所得:', y_list)
print('投資:', inv_list)
所得: [66.7, 63.3, 60.0, 56.7, 53.3, 50.0, 46.7, 43.3, 40.0, 36.7]
投資: [20.0, 19.0, 18.0, 17.0, 16.0, 15.0, 14.0, 13.0, 12.0, 11.0]
実質利子率が上昇すると,所得と投資は減少することがわかる。即ち,IS曲線上では所得と実質利子率は負の関係にあることが確認できる。
for
ループを使う:方法2#
n = len(r_list) # r_listの要素数
y_arr = np.zeros(n) # n個のゼロarray
inv_arr = np.zeros(n) # n個のゼロarray
for i in range(n):
a = np.array([[1-c, -1],[0,1]])
b = np.array([g-d,f-h*r_list[i]]) # i番目のr_listを使う
sol = np.linalg.solve(a,b)
y_arr[i] = sol[0].round(1) # i番目のy_arrに割り当てる
inv_arr[i] = sol[1] # i番目のinv_arrに割り当てる
print('所得:', y_arr)
print('投資:', inv_arr)
所得: [66.7 63.3 60. 56.7 53.3 50. 46.7 43.3 40. 36.7]
投資: [20. 19. 18. 17. 16. 15. 14. 13. 12. 11.]
IS曲線の導出:1つの方程式を使う#
上の2つの式を1つにする。
この式を使ってfor
ループで計算する。
y_list = []
for r in r_list:
y = (f-h*r+g-d)/(1-c)
y = round(y,1) # round(x)は浮動小数点を四捨五入する関数
y_list.append(y)
print('所得:', y_list)
所得: [66.7, 63.3, 60.0, 56.7, 53.3, 50.0, 46.7, 43.3, 40.0, 36.7]
次の方法でも良い。
n = len(r_list) # r_listの要素数
y_arr = np.zeros(n) # n個のゼロarray
for i in range(n):
y_arr[i] = (f-h*r_list[i]+g-d)/(1-c)
print('所得:', y_arr.round(1)) # round()は四捨五入するメソッド
所得: [66.7 63.3 60. 56.7 53.3 50. 46.7 43.3 40. 36.7]
IS曲線の導出:ベクトル演算を使う#
for
ループを使わずに,array
のベクトル演算を使うともっと簡単なコードで同じ結果を得ることができる。上の1つにまとめた式に基づいて次の関数を定義する。
def y_equil(r):
y = (f-h*r+g-d)/(1-c)
return y.round(1) # round()は四捨五入のメソッド
r_list
をarray
に変換する。
r_arr = np.array(r_list)
r_arr
array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
r_arr
を使って直接評価する。
print('所得:', y_equil(r_arr))
所得: [66.7 63.3 60. 56.7 53.3 50. 46.7 43.3 40. 36.7]