NumPy：高速化#
import py4macro
# 警告メッセージを非表示
import warnings
warnings.filterwarnings("ignore")
リストの用途は広く使いやすい。一方で、以下で説明する高速化のためのベクトル演算が使えない。それにより、コードが長くなり、比較的に実行速度が落ちることになる。それを補うのがNumPyである。NumPyは数値計算をする上で重要な役割を果たすパッケージであり，特に，行列計算に威力を発揮する。NumPyは「ナンパイ」と読む。慣例としてnpとしてインポートする。
import numpy as np
array#
基本となる関数がnp.array()であり，次のコードでは１次元配列を作る。
arr = np.array([10, 20, 30, 40, 50])
arr
array([10, 20, 30, 40, 50])
次に２次元配列を作成しよう。
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])
リストと見てくれはあまり変わらないが、乗用車をレースカーにチューンアップしたのと等しい。2次元arrayもリストから作ることができる。まず、入れ子になるリストを作成する。
lst = [l, l, l, l, l]
lst
[[1, 2, 3, 4, 5],
[1, 2, 3, 4, 5],
[1, 2, 3, 4, 5],
[1, 2, 3, 4, 5],
[1, 2, 3, 4, 5]]
np.array()関数の引数にlstを使うとarrayを作成できる。
np.array([lst])
array([[[1, 2, 3, 4, 5],
[1, 2, 3, 4, 5],
[1, 2, 3, 4, 5],
[1, 2, 3, 4, 5],
[1, 2, 3, 4, 5]]])
後述するが、ndarrayをリストに変換するメソッドも用意されている。
arrayについてこのサイトでは図を使って説明しているので興味があるひとは参考にしてはどうだろうか。
１次元配列#
要素の抽出#
１次元配列の場合、リストと同じ方法で要素を抽出できる。即ち、[]と要素のインデックスを使う。1番目の要素を抽出してみよう。
arr[1]
20
また、要素を連続で抽出するスライシング（slicing）も使うことができる。
arr[1:3+1]
array([20, 30, 40])
リストの場合と同じ様に，:の左側が選択する最初の要素であり，:の右側が選択する最後の次の番号である（即ち，:の右側の番号の要素は含まれない。
一方で、リストとは異なり、インデックスをリストとして指定することにより、複数の要素を簡単に抽出するできる。
arr[[0,1,3]]
array([10, 20, 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()）も同じである。
２次元配列：要素の抽出#
第\(i\)行目の第\(j\)列目の要素にアクセスするためにはmat[i,j]と書く。[,]の,を挟んで左が行を表し，右が列を示す。
[行のインデックス,列のインデックス]
コードはキーストロークが少なく，簡単なものが良いと考えられている。一方で，Pythonは初心者に比較的簡単な言語だと言われるが，それでも関数・メソッドの数は膨大である。そのため初心者にとって関数の使い方やオプションの書き方を間違う可能性は小さくない。さらに，自分が書いたコードを数週間・数ヶ月後に読み直すと，何をしようとしているのか分からないという状況が生じることもある。従って，以下の点が非常に重要になる。
間違いにくいコードの書き方を覚える。
高い可読性を意識したコードの書き方を覚える。
このスタンスに基づいて以下のルールに従って説明する。
[,]内の,を省略可能な場合であっても省略しない。[,]内の,の左または右を省略できる場合であっても省略しない。
行または列を連続選択する（slicing）場合を考えよう。上でも説明したように:を使うが
start:end
となる。ここでstartとは選択する要素の最初インデックスであり，endは選択する最後の要素の次のインデックスである（リストの場合と同じ）。このルールに従うと，
:の左を省略すると「最初から」という意味になる。:の右を省略すると「最後まで」という意味になる。,の左または右が:のみの場合，「全て」という意味になる。
これを読むだけでは分かりにくいと思うので，以下の例に目を通してもう一度この箇所の説明を読んむことを推奨する。
0行目、1列目の要素を抽出するとしよう
mat[0,1]
2
第\(i\)行の抽出はmat[i,:]で可能である。:は「列の全て」という意味になる。
r0 = mat[0,:]
r0
array([1, 2, 3, 4])
抽出した行は上で説明した１次元配列なのでインデックスやスライシングを使って要素にアクセスすることができる。
複数の行の抽出は次の方法で可能である。
mat[1:3,:]
array([[ 5, 6, 7, 8],
[ 9, 10, 11, 12]])
第３行目は含まれないことに注意しよう。リストの要素の取り出し方と同じように:の右のインデックスの行は含まれない。
（注意）
,:を省略してmat[1:3]と書いてもエラーは発生せず同じ結果が返されるが，,:があることにより，行を抽出しており列は全て選択されていることが明示的になる。
第\(i\) 列目を抽出したい場合は次のコードになる。
mat[:, 1]
array([ 2, 6, 10, 14])
複数列の抽出は以下のようにする。
mat[:, 1:3]
array([[ 2, 3],
[ 6, 7],
[10, 11],
[14, 15]])
arrayを使った計算#
まず２つ２次元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]])
乗算：バージョン１#
*を使うと要素どうしの積となる。数学で学んだ行列の計算とは異なるので注意すること。
m1*m2
array([[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サブパッケージ#
２次元配列の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使用時によく使う属性とメソッド#
py4macro.see()関数を使ってm3の属性とメソッドを調べてみる。
py4macro.see(m3)
.T .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[66], 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.])
この機能はベクトル演算（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の2次元配列でも有効である。
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\)
２つの式で表す。
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曲線の導出：２つの方程式を使う#
forループを使う：方法１#
上の４５度線モデルの投資を次のように変更する。
投資関数：\(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]
それぞれの実質利子率の値での均衡での所得と投資を計算する。まず２つの式で表す。
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ループを使う：方法２#
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曲線の導出：１つの方程式を使う#
上の２つの式を１つにする。
この式を使って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のベクトル演算を使うともっと簡単なコードで同じ結果を得ることができる。上の１つにまとめた式に基づいて次の関数を定義する。
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]