NumPyモジュールの基本的な使用法を説明する.

基本的なn次元配列の生成 (zeros, ones, arange関数)

NumPyは科学技術計算のための基本モジュールであり,リストに似た配列を扱う.

リストとの違いは:

  • サイズが変更できない.

  • 同じ種類のものしか保管できない(既定値は浮動小数点数).

zerosは$0$が入ったベクトルや行列を生成するための関数である. 引数は形状(shape: 配列の大きさ)で返値はNumPyのndarray($n$次元配列)オブジェクトである.

import numpy as np  # まずはインポート;別名は np が標準
z = np.zeros(5)  # 長さ5の0ベクトルを生成
z
array([0., 0., 0., 0., 0.])

同様に,ones関数ですべての要素が $1$ のベクトルや行列を生成できる.

e = np.ones(10)  # 長さ10の単位ベクトル
e
array([1., 1., 1., 1., 1., 1., 1., 1., 1., 1.])

$0,1,2,\cdots,n-1$ の数列はarange関数で生成できる.これはPython標準のrange関数に相当する.

$n$ 行 $m$ 列の行列を生成するためには,引数をタプル $(n,m)$ とすれば良い.

seq = np.arange(10)  # 0,1,2,...9 の数列
seq
array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
Z = np.zeros((2, 2))  # 2行2列の0行列を生成;引数に行数,列数のタプルを入力
Z
array([[0., 0.],
       [0., 0.]])

行列 $A$ の $i$ 行 $j$ 列の要素は,$A[i,j]$ でアクセスでき,代入も可能である.

(Pythonではインデックスは 0 から始まることに注意.)

たとえば,$3 \times 3$ の単位行列は最初に$3 \times 3$ の $0$ 行列を作って,後から対角成分に $1$ を代入することにより生成することができる.

I = np.zeros( (3,3) ) 
I[0,0] =1. 
I[1,1] =1. 
I[2,2] =1.
>>> 結果
array([[ 1.,  0.,  0.],
       [ 0.,  1.,  0.],
       [ 0.,  0.,  1.]])

実は単位行列はeye関数を使えば1行で生成できる.

I = np.eye(3)

もしくは,以下のようにfor文による反復を用いた方が,より一般的である.

I = np.zeros((3, 3))
for i in range(3):
    I[i, i] = 1.0
I
array([[1., 0., 0.],
       [0., 1., 0.],
       [0., 0., 1.]])

問題

  1. 長さ $10$ の $1$ が並んだベクトルを生成せよ.

  2. $3$ 行 $4$ 列の $1$ だけが入った行列を生成せよ.

  3. $(0,1)$, $(1,2)$, $(2,0)$ の成分だけが $1$ で,他の要素が $0$ の $3 \times 3$ 行列を生成せよ.

  4. $i$ 行,$i+1$ 列だけ $1$ で,他の要素が $0$ の $10 \times 10$ 行列を生成せよ. (ただし最後の行はすべて $0$ である.)

ベクトルと行列の生成(array関数の使用法)と形状(shape属性)

数字を自分で入れたベクトルや行列はNumPyのarray関数で生成できる.この関数の返値は多次元配列(ndarray)のオブジェクトである.

引数にリストを入れて,

v = np.array( [4,5,6] )

とすればベクトル $v = (4,5,6)$ が生成される.

行列の場合には,行ごとのリストをもう1つのリストで入れ子にすることによって生成される.

たとえば,

A  = np.array( [ [1,0,0],
               [-1, 1, 0],
               [0, -1, 1] ] )

とすれば,$3$ 行 $3$ 列の行列 $A$ が生成される.

また,引数dtypeでデータのタイプを指定できる.既定値は浮動小数点数であるが,上のようにすべて整数を入れると整数の要素をもつベクトルや配列が生成される.浮動小数点数を指定したい場合には,dtype=floatとする.

v = np.array( [4,5,6], dtype=float) # array([4.,5.,6.])

NumPyの多次元配列は形状(shape)属性をもつ.たとえば,上で生成した行列 $A$ の形状は $A.shape = (3,3)$ である. これは$3$行$3$列の行列を意味する.

NumPyの多次元配列の形状(shape)はshape属性にタプルを代入することによって変更できる. ただし要素数(size)の変更はできない.たとえば,行列 $A$ を 1行9列の行列に変更するには,

A.shape=(1,9) # array([[ 1.,  0.,  0., -1.,  1.,  0.,  0., -1.,  1.]])

とすれば良い.

上のベクトル $v$ の形状 $v.shape$ はタプル $(3,)$ である.これは長さ $3$ のベクトル(列がない行列)を表す. $1$ 行 $3$ 列の行列とは異なることに注意されたい.$1$ 行 $3$ 列に変更するには,

v.shape=(1,3) # array( [[4, 5, 6]] )

とする.よく見ると外側に $[]$ が1つ増えている(元は $array( [4, 5, 6] )$ だった).

五十音によるNumPyの例

x = np.array(["あ", "い", "う", "え"])
print(x.shape)
x = np.array([["あ", "い", "う", "え"]])
print(x.shape)
x = np.array([["あ"], ["い"], ["う"], ["え"]])
print(x.shape)
x = np.array([["あ", "い"], ["う", "え"]])
print(x.shape)
x = np.array([[["あ"]], [["い"]], [["う"]], [["え"]]])
print(x.shape)
(4,)
(1, 4)
(4, 1)
(2, 2)
(4, 1, 1)

問題

  1. 以下の行列 $A$ を $9$ 行 $1$ 列の形状に変更せよ.
  2. 以下のベクトル $v$ を$3$ 行 $1$ 列の形状に変更せよ.カッコがどのように変更されたか観察せよ.
  3. 以下の行列 $A$ を $4$ 行 $2$ 列の形状に変更することを試みよ.何が発生するか観察せよ.
A = np.array([[1, 0, 0], [-1, 1, 0], [0, -1, 1]])
v = np.array([4, 5, 6])

添え字(インデックス)とスライシングと演算

添え字とスライス表記はリストと同じ.

L=[1,2,3,4] 
L[2]     # 3
L[1:3]   # [2,3]
x= np.array([1,2,3,4]) 
x[2]     # 3
x[1:3]   # array( [2,3] )

加算の結果は異なり,リストは結合,arrayの場合は要素ごとの和になる.(他の演算子に対しても同じである.)

[1,2,3] + [4,5,6]                         # [1,2,3,4,5,6]
np.array( [1,2,3] ) + np.array( [4,5,6] ) # [5,7,9]

NumPyのarray(ベクトル)とスカラーの和は,(後述するブロードキャストが行われ)スカラーをベクトルにコピーした後に和をとり,ベクトルとなる. (他の演算子に対しても同じである.)

np.array( [1,2,3] ) + 1    # [2,3,4] 
np.array( [1,2,3] ) * 10   # [10,20,30]

五十音によるNumPyの例

import numpy as np

L = np.array(["あ", "い", "う", "え", "お", "か"])

print("L[1] =", L[1])
print("L[-2] =", L[-2])
print("L[1:4] =", L[1:4])
print("L[:-2:2] =", L[:-2:2])
L[1] = い
L[-2] = お
L[1:4] = ['い' 'う' 'え']
L[:-2:2] = ['あ' 'う']

問題

  1. 以下のベクトル $x,y$ の和を計算せよ.
  2. $x$ の最初の2つの要素と $y$ の最後の2つの要素の和を計算せよ.
  3. ベクトル $x$ とスカラー $10$ との和を計算せよ.
  4. ベクトル $x$ をスカラー $3$ で乗じた値を計算せよ.
  5. $x$ と $y$ の積を計算せよ.(内積ではなく,要素ごとの積をとれ.)
x = np.array([1, 2, 3])
y = np.array([5, 6, 7])

行列に対する添え字(インデックス)とスライシングと演算

2次元の配列(行列)の添え字は配列[行番号,列場号]でアクセスできる.

X = np.array( [[1, 2, 3],
               [4, 5, 6],
               [7, 8, 9]])
X[0,1]   # 2

スライシングはリストのときと同様に,開始番号:終了番号:ステップサイズ で開始番号から終了番号 $-1$ まで刻み幅($=$ステップサイズ)で切り出しを行う.

行と列に対して別々に切り出しを行うことができる.たとえば,すべての行(:)と1列目を指定すると,1列目だけを切り出すことができる.

X[ : , 1]  #array([2, 5, 8])

同様に1行目だけを切り出すには,

X[ 1, : ]  #X[1] としても同じ!  array([4, 5, 6]) を返す.

とする.

演算もベクトルと同様に行われる.スカラーに対する演算は,後述するブロードキャスト(同じ大きさになるようにコピーしてから演算)が行われる.

Image("../figure/numpy_slice.PNG", width=800)
import numpy as np

X = np.array(
    [["あ", "い", "う", "え", "お"], ["か", "き", "く", "け", "こ"], ["さ", "し", "す", "せ", "そ"]]
)

print("X[0,1] =", X[0, 1])
print("X[:, 1] =", X[:, 1])
print("X[1, 2:4] =", X[1, 2:4])
print("X[1:, 1:] =", X[1:, 1:])
print(X[::2, 1::2])
X[0,1] = い
X[:, 1] = ['い' 'き' 'し']
X[1, 2:4] = ['く' 'け']
X[1:, 1:] = [['き' 'く' 'け' 'こ']
 ['し' 'す' 'せ' 'そ']]
[['い' 'え']
 ['し' 'せ']]

問題

  1. 行列 $X$ の2行1列目の要素を出力せよ.
  2. 行列 $X$ の2行目だけを切り出せ.
  3. 行列 $X$ の0から1行と1から2列から成る部分行列 $[[2,3], [5,6]]$ を切り出せ.
  4. 行列 $X$ の0行目と2行目,0列目と2列目から成る部分行列 $[[1,3][7,9]]$を切り出せ.
  5. $X$ と $Y$ の和を求めよ.
  6. $X$ と $Y$ の積を求めよ.また $X$ とスカラー $1$ の積,ベクトル $(1,1,1)$ との積と較べてみよ.
X = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
Y = np.ones((3, 3))

問題

  1. 行列 $A$ を以下のように生成せよ.

    A = np.arange(25)
    A.shape = (5,5)
  2. 行列 $A$ から

    array([[ 5,  7,  9],
        [20, 22, 24]])

    を.1行で切り出せ.

  3. 行列 $A$ から

    array([[ 0,  2],
        [10, 12],
        [20, 22]])

    を.1行で切り出せ.

array([[ 0,  2],
       [10, 12],
       [20, 22]])

ユニバーサル関数

NumPyの以下の関数はユニバーサル関数である.

  • Arithmetic Operators: + - * / // % **
  • Bitwise Operators: & | ~ ^ >> <<
  • Comparison Oper’s: < > <= >= == !=
  • Trig Family: np.sin, np.cos, np.tan ...
  • Exponential Family: np.exp, np.log, np.log10 ...

これらを用いることによって多次元配列の各要素に対する演算を高速に行うことができる.

たとえば,0から10000-1の整数の配列に対して,すべての要素に $5$ を加えるには,以下のようにする.

x = np.arange(10000) #0から10000-1の整数の配列を生成
x + 5 # スカラーを加えると後述のブロードキャストによって全ての要素にスカラーが足される

同じ配列に対して平方根をとるには,

np.sqrt(x)

とする.NumPyモジュールのsqrt関数はユニバーサル関数なので,すべての要素に対して平方根を計算して返す. (数学モジュールmathの平方根だとエラーする.)

x = np.arange(10000)
np.sqrt(x)
array([ 0.        ,  1.        ,  1.41421356, ..., 99.98499887,
       99.9899995 , 99.99499987])

問題

  1. 0から99の整数が入った配列を生成し,正弦(sin)を計算せよ.
  2. 1から999の整数が入った配列の要素を2倍し,それに対数(log)をとったものに100を加えた配列を生成せよ.

擬似乱数発生サブモジュール random

通常のrandomモジュールと類似の擬似乱数がNumPyでも生成できる.

ただし引数sizeでランダムな要素をもつ配列を一度に生成できる.

たとえば,5つの[0,1)の一様乱数の生成は以下のようになる.

np.random.random(size=5) # array([ 0.19672059,  0.91704669,  0.05184376,  0.13490049,  0.13663051])

整数の乱数を生成するにはrandintを用いる.引数は通常の下限(low),上限(high)の他に,sizeを指定できる.(標準のrandomモジュールと異なり上限は含まれないことに注意.) sizeにタプルを入れると多次元の配列が生成される.

np.random.randint(1, 6, size=(3,3))
>>> 結果
array([[4, 5, 5],
       [1, 4, 1],
       [1, 4, 1]])

正規分布にしたがうランダムな値を得るにはnormalが使える.

np.random.normal(100,10,(2,2)) #平均100,標準偏差10にしたがう正規乱数を2行2列の行列に入れる
>>> 結果
array([[ 92.24228997,  85.74292368],
       [ 83.08779098,  99.39865564]])

他にも様々な分布にしたがうランダムな値を生成できる.

問題

  1. さいころを $10$回振ったシミュレーション結果を表す配列を生成せよ.
  2. 平均 $10$,標準偏差 $10$ の正規乱数を $10 \times 10$ の行列に入れよ.

問題

平均 $100$, 標準偏差 $10$ の需要を $10$個入れた配列 demand を生成せよ.

集約関数

NumPyの関数で多次元配列内の値を集約して計算してくれるものがある.

min, max, sum, mean, average, argmin, argmax, std などなど
x = np.random.randint(low=1, high=100, size=10000000)
x
array([60, 32, 83, ..., 42, 97, 99])
np.max(x)  # 最小値を求める
99
np.argmin(x)  # 最小値になるインデックス
117
np.sum(x)  # 合計を求める
499809263
np.std(x)  # 標準偏差を求める
28.579390628457567
A = np.random.randint(low=1, high=100, size=(3, 5))  # ランダムな要素を含んだ行列 A を生成
A
array([[90, 44, 30, 39, 68],
       [24, 42, 94, 63, 65],
       [39, 34, 61, 80, 17]])
np.sum(A, axis=0)  # 行方向(軸番号=0)で合計
array([153, 120, 185, 182, 150])
np.sum(A, axis=1)  # 列方向(軸番号=1)で合計
array([271, 288, 231])

問題

  1. 上の行列Aの行ごとの最小値を求めよ.
  2. 同じ行列の列ごとの最小値を求めよ.
  3. さいころを $10$回振ったときの結果の配列を作成し,試行の中の最小値,最大値,平均値を出力せよ.
  4. 平均 $10$,標準偏差 $10$ の正規乱数を $10 \times 10$ の行列に入れ,全体の平均,行ごとの平均を出力せよ.

ブロードキャスト

多次元配列 = 多次元配列 + スカラー

になるのは,同じ形状になるように変形してから演算を行うためである.

x= np.arange(5)  # [0 1 2 3 4]  形状は (5,)
x + 10           # 結果は array([10, 11, 12, 13, 14])

スカラー10をxと同じ形状(5,)になるようにコピーしてから加算をしている.

x+ np.array([10,10,10,10,10]) #スカラー10を1次元配列にしてから,同じ長さになるまでコピー(ブロードキャスト)

同様に3行3列の行列 A と長さ3のベクトルxの加算を試してみる.

A = np.ones((3, 3), dtype=int)
x = np.array([10, 20, 30])
print("A=", A)
print("x=", x)
print(A.shape, x.shape)
A= [[1 1 1]
 [1 1 1]
 [1 1 1]]
x= [10 20 30]
(3, 3) (3,)
A + x  # ブロードキャストを用いた加算
array([[11, 21, 31],
       [11, 21, 31],
       [11, 21, 31]])
A + np.array([[10, 20, 30], [10, 20, 30], [10, 20, 30]])  # [10,20,30] のコピーを作成してから加算
array([[11, 21, 31],
       [11, 21, 31],
       [11, 21, 31]])
x = np.array([[1, 2, 3]]).T  # 3行1列の行列(縦ベクトル)
y = np.array([10, 20, 30])  # 長さ3の横ベクトル
print("x=", x)
print("y=", y)
print(x.shape, y.shape)
x= [[1]
 [2]
 [3]]
y= [10 20 30]
(3, 1) (3,)
x * y
array([[10, 20, 30],
       [20, 40, 60],
       [30, 60, 90]])
x + y  # ブロードキャストによる加算
array([[11, 21, 31],
       [12, 22, 32],
       [13, 23, 33]])
np.array([[1, 1, 1], [2, 2, 2], [3, 3, 3]]) + np.array(
    [[10, 20, 30], [10, 20, 30], [10, 20, 30]]
)  # コピーを作成してから加算
array([[11, 21, 31],
       [12, 22, 32],
       [13, 23, 33]])

インデックス配列

NumPyでは,添え字を別の配列で指定できる.この機能はpandasでよく使うので,慣れておく必要がある.

添え字の配列(もしくはリスト)を用いて配列 $x$ から添え字に対応する要素を取り出す.

x = np.array( [1,2,3,4,5,6] )
ind = np.array( [1,3,5] )  # リスト [1,3,5] でも同じ
x[ind]

>>> 結果
array([2, 4, 6])

五十音によるNumPyの例

import numpy as np

L = np.array(["あ", "い", "う", "え", "お", "か"])
ind = np.array([0, 2, 5])
L[ind]
array(['あ', 'う', 'か'], dtype='<U1')

Boolインデックス配列の例

NumPyの配列に対して要素がTrueかFalseになる条件式を書くと,要素がTrueかFalseの配列が生成される. これをBoolインデックス配列とよぶ.

x = np.array( [5,1,6,9,3,4] )

x >= 5 
>>> 
array([ True, False,  True,  True, False, False])

x[ x>=5 ]
>>>
array([5, 6, 9])

配列から奇数だけを抜き出すことも、同様にできる。

x = np.array( [1,2,3,4,5,6] )
x%2 == 1 #2で割ったときの余りが1か?

>>> 結果
array([ True, False,  True, False,  True, False])

この配列の要素がTrueの部分だけを切り出すことができる.

x = np.array( [1,2,3,4,5,6] )
bool_index = x%2 == 1 
x[ bool_index ]

>>> 結果
array([1, 3, 5])

これを1行で行うと

x[ x%2==1 ]

と書ける.

複数の論理条件を入れる場合には,集合(set)に対する演算子を用いる.

かつ(and)は &,または(or)は | である.論理条件には必ず( )をつけて区分を明確にする. たとえば,偶数で $4$以上のものを切り出すには,

x[ (x%2==0) & (x>=4) ] # array([4, 6])

とする.

例として,ランダムに生成した長さ10の配列から,条件を満たすものを切り出してみる.

x = np.random.randint(low=1, high=100, size=10)  # [1,99]の一様乱数による配列
x
array([47, 70, 26, 98,  6, 16, 19, 62,  8, 77])
bool_index = x < 30  # 条件を満たすなら True の配列(Boolインデックス配列)
bool_index
array([False, False,  True, False,  True,  True,  True, False,  True,
       False])
x[bool_index]  # 配列 bool_index がTrueのインデックスから成る配列
array([26,  6, 16, 19,  8])
x[x < 30]  # 上の操作を1行で書く
array([26,  6, 16, 19,  8])
x[(x <= 30) & (x >= 10)]  # 30以下でかつ10以上のもの(and は & を使う;論理条件には( )をつける)
array([26, 16, 19])
x[(x < 30) | (x % 2 == 0)]  # 30以下かもしくは2で割り切れるもの
array([70, 26, 98,  6, 16, 19, 62,  8])

問題

NumPyの配列 $x$ を以下のように生成する.

x = np.array( [1,2,3,4,5,6] )
  1. 添え字が $0$ 番目と $5$ 番目の要素をインデックス配列を用いて取り出せ.
  2. 要素が2で割り切れるものだけをTrueかFalseが入った配列(Bool添え字配列)を用いて取り出せ.
  3. 要素が奇数のものだけを2倍せよ.

問題

  1. さいころを $10$回振ったときの結果の配列を作成し,試行の中で $4$ 以上のものだけを切り出せ.
  2. さいころを $100$回振ったときの結果の配列を作成し,試行の中で $5$ 以上で奇数のものだけを切り出せ.
  3. 平均 $10$,標準偏差 $10$ の正規乱数を $10 \times 10$ の行列に入れ,負のものだけを切り出せ.

すべて (all), いずれか (any), 条件指定 (where)

BOOLインデックス配列を引数とした関数として, all, any, where がある.

  • allは配列の中の要素がすべてTrueのときTrueを返す.
np.all([True,True,False])   => False
  • anyは配列の中の要素に1つでもTrueのものがあればTrueを返す.
np.any([True,True,False])   => True
  • whereは配列がTrueのインデックスの配列を返す.
np.where([True,False,True]) => array([0, 2])

whereは3項演算子 where(BOOLインデックス配列, Trueのときの値, Falseのときの値) としても使える.

x = np.array([1,2,3,4,5])
np.where( x%2, "odd", "even")

は以下の配列を返す.

['odd', 'even', 'odd', 'even', "odd"]
x = np.random.randint(low=1, high=100, size=10)  # [1,100]の一様乱数による配列
print(x)
np.all(x <= 90)  # すべてがTrueのときTrueを返す
[70 36 44 87 34 74 50 39 32 19]
True
np.any(x <= 10)  # いずれかの要素がTrueのときTrueを返す
False

問題

  1. 平均 $100$,標準偏差 $10$ の需要を $10$個($10$日分)入れた配列 demand を作り, 需要が仕入れ量 $110$ より多い場合には品切れ費用, 少ない場合には在庫費用を入れた配列 costを作れ. ただし,1個あたりの品切れ費用は $100$円, 在庫費用は $1$円とする.

  2. $10$日間の費用の合計を計算せよ.

  3. 仕入れ量を変えた場合の費用を計算し,適切な仕入れ量を求めよ.

問題

次のプログラムは入力した正の整数 x に対して何を返すためのプログラムか? 推測せよ.

not np.any([x%i == 0 for i in range(2, x)])

数列の生成法

NumPyには様々な数列を生成するための関数が準備されている.

通常のリストと同じように整数の配列を生成するためには,arange関数を用いる.

引数は,開始,終了,ステップ(階差)である.

x = np.arange(1,11,2) # array([1, 3, 5, 7, 9])

配列のサイズ(要素の個数)を指定して階差が等しくなるような数列を生成するには,linspace関数を用いる.

引数は,開始,終了,要素数であるが,既定値だと最後の要素が終了と一致するように生成することに注意されたい.

x = np.linspace(1,6,10)
>>> 結果
array([ 1.        ,  1.55555556,  2.11111111,  2.66666667,  3.22222222,
        3.77777778,  4.33333333,  4.88888889,  5.44444444,  6.        ])

linspace関数はmatplotlibで関数の図を描画するときによく用いられる.以下に $y=x^2$ のグラフを描画した例を示す.

x = np.linspace(0, 100, 100)  # 0から1まで均等割りした100個の要素から成る配列
y = x**2  # 同じ長さの配列に xの2乗を入れる
import matplotlib.pyplot as plt  # 図を描画する準備

%matplotlib inline
plt.plot(x, y)  # 描画
[<matplotlib.lines.Line2D at 0x7fad605cb250>]

問題

  1. $x$ に1から100まで公差3の数列を代入し,$x$ の平方根のグラフをプロットせよ.ただし,平方根はNumPyのsqrt関数で計算せよ.(NumPyモジュールのsqrt関数はユニバーサル関数である.mathモジュールのsqrt関数を用いると,スカラーしか計算できないことに注意.)

  2. 0.1から10まで1000個の数値を均等に配置した $x$ 座標に対して $y=1/x$ の関数を描画せよ.

import matplotlib.pyplot as plt  # 図を描画する準備

%matplotlib inline

まとめ問題

  1. あるスーパーで販売しているある商品の在庫を管理したい. いま,$t-1$日の在庫量の営業後に在庫を調べて,発注量 $x[t]$ を決める.発注した商品は翌日の朝に届き, それによって $t$ 日の営業後の在庫 $I[t]$ は,以下のように決まるものとする.
$$ I[t] = I[t-1]+x[t] –D[t] $$

ここで $D[t]$ は$t$日の需要量であり,需要が平均 $100$,標準偏差 $10$ の正規分布とする. また,最初の日($t=0$)の営業後の在庫は $0$ であると仮定する.

以下の在庫管理方策のプログラムを作成せよ.

  • (s,S)方策:在庫量 $I[t-1]$ が $s$ 未満であるなら,翌朝の在庫が $S$ になるように発注する.
  1. $10$箇所のスーパーがあり,同じ平均と標準偏差の需要があるとする. NumPyの多次元配列を使って,一度にすべての店舗に対するシミュレーションを行うプログラムを作成せよ.
CPU times: user 723 ms, sys: 94.8 ms, total: 818 ms
Wall time: 818 ms