NumPy 入門

NumPy 入門#

本章では、Python で数値計算を高速に行うためのライブラリ [1] である NumPy の使い方を学びます。本章の目標は、単回帰分析と重回帰分析で解説されている重回帰分析を行うアルゴリズムをNumPy を用いて実装することです。

NumPy による 多次元配列（multidimensional array） の扱い方を知ることは、他の様々なライブラリを利用する際に役立ちます。例えば、様々な機械学習手法を統一的なインターフェースで利用できる scikit-learn や、ニューラルネットワークの記述・学習を行うためのフレームワークである PyTorch は、NumPy に慣れておくことでとても使いやすくなります。

それでは、まず NumPy の基礎的な使用方法を説明します。

NumPy を使う準備#

NumPy は Google Colaboratory（以下 Colab）上のノートブックにはデフォルトでインストールされているため、ここではインストールの方法は説明しません。自分のコンピュータに NumPy をインストールしたい場合は、こちらを参照してください。：Installing packages

Colab 上ではインストール作業は必要ないものの、ノートブックを開いた時点ではまだ numpy モジュールが読み込まれていません。ライブラリの機能を利用するには、そのライブラリが提供するモジュールを読み込む必要があります。

例えば A というモジュールを読み込みたいとき、一番シンプルな記述方法は import A です。ただ、もし A というモジュール名が長い場合は、import A as B のようにして別名を付けることができます。 as を使って別名が与えられると、以降そのモジュールはその別名を用いて利用することができます。 import A as B と書くと、A というモジュールは B という名前で利用することができます。これは Python の機能なので NumPy 以外のモジュールを読み込みたい場合にも使用可能です。

慣習的に、numpy にはしばしば np という別名が与えられます。コード中で頻繁に使用するモジュールには、短い別名をつけて定義することがよく行われます。

それでは、numpy を np という名前で import してみましょう。

import numpy as np

多次元配列を定義する#

ベクトル・行列・テンソルなどは、プログラミング上は多次元配列により表現でき、NumPy では ndarray というクラスで多次元配列を表現します。早速、これを用いてベクトルを定義してみましょう。

# ベクトルの定義
a = np.array([1, 2, 3])

a

array([1, 2, 3])

このように、Python リスト [1, 2, 3] を np.array() に渡すことで、\([1, 2, 3]\) というベクトルを表す ndarray オブジェクトを作ることができます。 ndarray オブジェクトは shape という属性（attribute） を持っており、その多次元配列の形（shape） が保存されています。上で定義した a という ndarray オブジェクトの形を調べてみましょう。

a.shape

(3,)

(3,) という要素数が 1 の Python のタプルが表示されています。 ndarray の形は、要素が整数のタプルで表され、要素数はその多次元配列の次元数（dimensionality, number of dimensions） を表します。形は、その多次元配列の各次元の大きさを順に並べた整数のタプルになっています。

次元数は、ndarray の ndim という属性に保存されています。

a.ndim

これは、len(a.shape) と同じ値になります。今、a という ndarray は 1 次元配列なので、a.shape は要素数が 1 のタプルで、ndim の値は 1 でした。

では次に、\(3 \times 3\) 行列を定義してみましょう。

# 行列の定義
b = np.array(
    [[1, 2, 3],
     [4, 5, 6],
     [7, 8, 9]]
)

b

array([[1, 2, 3],
       [4, 5, 6],
       [7, 8, 9]])

形と次元数を調べます。

print('Shape:', b.shape)
print('Rank:', b.ndim)

Shape: (3, 3)
Rank: 2

ここで、size という属性も見てみましょう。

b.size

これは、b という ndarray が持つ要素の数を表しています。 b は \(3 \times 3\) 行列なので、要素数は 9 です。 「形」「次元数」「サイズ」という言葉がそれぞれ意味するものの違いを確認してください。

NumPy の ndarray の作成方法には、np.array() を用いて Python のリストから多次元配列を作る方法以外にも、色々な方法があります。以下に代表的な例をいくつか紹介します。

# 形を指定して、要素が全て 0 で埋められた ndarray を作る
a = np.zeros((3, 3))

a

array([[0., 0., 0.],
       [0., 0., 0.],
       [0., 0., 0.]])

# 形を指定して、要素が全て 1 で埋められた ndarray を作る
b = np.ones((2, 3))

b

array([[1., 1., 1.],
       [1., 1., 1.]])

# 形と値を指定して、要素が指定した値で埋められた ndarray を作る
c = np.full((3, 2), 9)

c

array([[9, 9],
       [9, 9],
       [9, 9]])

# 指定された大きさの単位行列を表す ndarray を作る
d = np.eye(5)

d

array([[1., 0., 0., 0., 0.],
       [0., 1., 0., 0., 0.],
       [0., 0., 1., 0., 0.],
       [0., 0., 0., 1., 0.],
       [0., 0., 0., 0., 1.]])

# 形を指定して、 0 ~ 1 の間の乱数で要素を埋めた ndarray を作る
e = np.random.random((4, 5))

e

array([[0.49165323, 0.36240804, 0.38689295, 0.21150064, 0.52094163],
       [0.32896834, 0.58636557, 0.1148585 , 0.75942012, 0.50527133],
       [0.89393832, 0.89432732, 0.39991541, 0.27745812, 0.33116905],
       [0.19885229, 0.86450159, 0.71562185, 0.59215422, 0.46585769]])

# 3 から始まり 10 になるまで 1 ずつ増加する数列を作る（10 は含まない）
f = np.arange(3, 10, 1)

f

array([3, 4, 5, 6, 7, 8, 9])

多次元配列の要素を選択する#

前節では NumPy を使って多次元配列を定義するいくつかの方法を紹介しました。本節では、作成した ndarray のうちの特定の要素を選択して、値を取り出す方法を紹介します。最もよく行われる方法は [] を使った添字表記（subscription） による要素の選択です。

整数による要素の選択#

例えば、上で作成した e という \(4 \times 5\) 行列を表す多次元配列から、1 行 2 列目の値を取り出すには、以下のようにします。

val = e[0, 1]

val

0.36240803558837553

「1 行 2 列目」を指定するのに、インデックスは [0, 1] でした。これは、NumPy の ndarray の要素は Python リストと同じく、添字が 0 から始まるゼロベースインデックス（zero-based index） が採用されているためです。つまり、この行列の i 行 j 列目の値は、[i - 1, j - 1] で取り出すことができます。

スライスによる要素の選択#

NumPy の ndarray に対しても、Python のリストと同様にスライス表記（slicing） を用いて選択したい要素を範囲指定することができます。 ndarray はさらに、カンマ区切りで複数の次元に対するスライスを指定できます。

# 4 x 5 行列 e の真ん中の 2 x 3 = 6 個の値を取り出す
center = e[1:3, 1:4]

center

array([[0.58636557, 0.1148585 , 0.75942012],
       [0.89432732, 0.39991541, 0.27745812]])

前節最後にある e の出力を見返すと、ちょうど真ん中の部分の \(2 \times 3\) 個の数字が取り出せていることが分かります。ここで、e の中から [1, 1] の要素を起点として 2 行 3 列を取り出して作られた center の形を、e の形と比較してみましょう。

print('Shape of e:', e.shape)
print('Shape of center:', center.shape)

Shape of e: (4, 5)
Shape of center: (2, 3)

また、インデックスを指定したり、スライスを用いて取り出した ndarray の一部に対し、値を代入することもできます。

# 先程の真ん中の 6 個の値を 0 にする
e[1:3, 1:4] = 0

e

array([[0.49165323, 0.36240804, 0.38689295, 0.21150064, 0.52094163],
       [0.32896834, 0.        , 0.        , 0.        , 0.50527133],
       [0.89393832, 0.        , 0.        , 0.        , 0.33116905],
       [0.19885229, 0.86450159, 0.71562185, 0.59215422, 0.46585769]])

整数配列による要素の選択#

ndarray の [] には、整数やスライスの他に、整数配列を渡すこともできます。整数配列とは、ここでは整数を要素とする Python リストまたは ndarray のことを指しています。

具体例を示します。まず、\(3 \times 3\) 行列を表す a という ndarray を定義します。

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

a

array([[1, 2, 3],
       [4, 5, 6],
       [7, 8, 9]])

この ndarray から、

1 行 2 列目：a[0, 1]
3 行 2 列目：a[2, 1]
2 行 1 列目：a[1, 0]

の 3 つの要素を選択して並べ、形が (3,) であるような ndarray を作りたいとします。

これは、以下のように、順に対象の要素を指定して並べて新しい ndarray にすることでももちろん実現できます。

np.array([a[0, 1], a[2, 1], a[1, 0]])

array([2, 8, 4])

しかし、同じことが選択したい行、選択したい列を、順にそれぞれリストとして与えることでも行えます。

a[[0, 2, 1], [1, 1, 0]]

array([2, 8, 4])

選択したい 3 つの値がどの行にあるかだけに着目すると、それぞれ 1 行目、3 行目、2 行目にある要素です。
ゼロベースインデックスでは、それぞれ 0, 2, 1 行目です。
これが a の [] に与えられた 1 つ目のリスト [0, 2, 1] の意味です。

同様に、列に着目すると、ゼロベースインデックスでそれぞれ 1, 1, 0 列目の要素です。
これが a の [] に与えられた 2 つ目のリスト [1, 1, 0] の意味です。

ndarray のデータ型#

1 つの ndarray の要素は、全て同じ型を持ちます。 NumPy では様々なデータ型を使うことができますが、ここでは一部だけを紹介します。 NumPy は Python リストを渡して ndarray を作る際などには、その値からデータ型を推測します。 ndarray のデータ型は、dtype という属性に保存されています。

# 整数（Python の int 型）の要素をもつリストを与えた場合
x = np.array([1, 2, 3])

x.dtype

dtype('int64')

# 浮動小数点数（Python の float 型）の要素をもつリストを与えた場合
x = np.array([1., 2., 3.])

x.dtype

dtype('float64')

以上のように、Python の int 型は自動的に NumPy の int64 型になりました。また、Python の float 型は自動的に NumPy の float64 型になりました。 Python の int 型は NumPy の int_ 型に対応づけられており、Python の float 型は NumPy の float_ 型に対応づけられています。この int_ 型はプラットフォームによって int64 型と同じ場合と int32 型と同じ場合があります。 float_ 型についても同様で、プラットフォームによって float64 型と同じ場合と float32 型と同じ場合があります。

特定の型を指定して ndarray を作成するには、以下のようにします。

x = np.array([1, 2, 3], dtype=np.float32)

x.dtype

dtype('float32')

このように、dtype という引数に NumPy の dtype オブジェクトを渡します。これは 32 ビット浮動小数点数型を指定する例です。同じことが、文字列で指定することによっても行えます。

x = np.array([1, 2, 3], dtype='float32')

x.dtype

dtype('float32')

これはさらに、以下のように短く書くこともできます。

x = np.array([1, 2, 3], dtype='f')

x.dtype

dtype('float32')

一度あるデータ型で定義した配列のデータ型を別のものに変更するには、astype を用いて変換を行います。

x = x.astype(np.float64)

x.dtype

dtype('float64')

多次元配列を用いた計算#

ndarray を使って行列やベクトルを定義して、それらを用いていくつかの計算を行ってみましょう。

ndarray として定義されたベクトルや行列同士の要素ごとの加減乗除は、Python の数値同士の四則演算に用いられる +、-、*、/ という記号を使って行えます。

それでは、同じ形の行列を 2 つ定義し、それらの要素ごとの加減乗除を実行してみましょう。

# 同じ形 (3 x 3) の行列を 2 つ定義する
a = np.array([
    [0, 1, 2],
    [3, 4, 5],
    [6, 7, 8]
])

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

# 足し算
c = a + b

c

array([[ 1,  3,  5],
       [ 7,  9, 11],
       [13, 15, 17]])

# 引き算
c = a - b

c

array([[-1, -1, -1],
       [-1, -1, -1],
       [-1, -1, -1]])

# 掛け算
c = a * b

c

array([[ 0,  2,  6],
       [12, 20, 30],
       [42, 56, 72]])

# 割り算
c = a / b

c

array([[0.        , 0.5       , 0.66666667],
       [0.75      , 0.8       , 0.83333333],
       [0.85714286, 0.875     , 0.88888889]])

NumPy では、与えられた多次元配列に対して要素ごとに計算を行う関数が色々と用意されています。以下にいくつかの例を示します。

# 要素ごとに平方根を計算する
c = np.sqrt(b)

c

array([[1.        , 1.41421356, 1.73205081],
       [2.        , 2.23606798, 2.44948974],
       [2.64575131, 2.82842712, 3.        ]])

# 要素ごとに値を n 乗する
n = 2
c = np.power(b, n)

c

array([[ 1,  4,  9],
       [16, 25, 36],
       [49, 64, 81]])

要素ごとに値を n 乗する計算は、以下のようにしても書くことができます。

c ** n

array([[   1,   16,   81],
       [ 256,  625, 1296],
       [2401, 4096, 6561]])

はじめに紹介した四則演算は、同じ大きさの 2 つの行列同士で行っていました。ここで、\(3 \times 3\) 行列 a と 3 次元ベクトル b という大きさのことなる配列を定義して、それらを足してみましょう。

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

b = np.array([1, 2, 3])

c = a + b

c

array([[ 1,  3,  5],
       [ 4,  6,  8],
       [ 7,  9, 11]])

形が同じ行列同士の場合と同様に計算することができました。

これは NumPy が自動的に ブロードキャスト（broadcast） と呼ばれる操作を行っているためです。これについて次節で説明します。

ブロードキャスト#

行列同士の要素ごとの四則演算は、通常は行列の形が同じでなければ定義できません。しかし、前節の最後では \(3 \times 3\) 行列に 3 次元ベクトルを足す計算が実行できました。

これが要素ごとの計算と同じように実行できる理由は、NumPy が自動的に 3 次元ベクトル b を 3 つ並べてできる \(3 \times 3\) 行列を想定し、a と同じ形に揃える操作を暗黙に行っているからです。この操作を、ブロードキャストと呼びます。

算術演算を異なる形の配列同士で行う場合、NumPy は自動的に小さい方の配列をブロードキャストし、大きい方の配列と形を合わせます。ただし、この自動的に行われるブロードキャストでは、行いたい算術演算が、大きい方の配列の一部に対して繰り返し行われることで実現されるため、実際に小さい方の配列のデータをコピーして大きい配列をメモリ上に作成することは可能な限り避けられます。また、この繰り返しの計算は NumPy の内部の C 言語によって実装されたループで行われるため、高速です。

よりシンプルな例で考えてみましょう。以下のような配列 a があり、この全ての要素を 2 倍にしたいとします。

a = np.array([1, 2, 3])

a

array([1, 2, 3])

このとき、一つの方法は以下のように同じ形で要素が全て 2 である別の配列を定義し、これと要素ごとの積を計算するやり方です。

b = np.array([2, 2, 2])

c = a * b

c

array([2, 4, 6])

しかし、スカラの 2 をただ a に掛けるだけでも同じ結果が得られます。

c = a * 2

c

array([2, 4, 6])

* 2 という計算が、c の 3 つの要素のどの要素に対する計算なのかが明示されていないため、NumPy はこれを全ての要素に対して行うという意味だと解釈して、スカラの 2 を a の要素数 3 だけ引き伸ばしてから掛けてくれます。

形の異なる配列同士の計算がブロードキャストによって可能になるためにはルールがあります。

それは、 「2 つの配列の各次元が同じ大きさになっているか、どちらかが 1 であること」 です。このルールを満たさない場合、NumPy は “ValueError: operands could not be broadcast together with shapes (1 つ目の配列の形) (2 つ目の配列の形)” というエラーを出します。

ブロードキャストされた配列の各次元のサイズは、入力された配列のその次元のサイズの中で最大の値と同じになっています。入力された配列は、各次元のサイズが入力のうち大きい方のサイズと同じになるようブロードキャストされ、その拡張されたサイズで計算されます。

もう少し具体例を見てみましょう。以下のような 2 つの配列 a と b を定義し、足します。

# 0 ~ 9 の範囲の値をランダムに用いて埋められた (2, 1, 3) と (3, 1) という大きさの配列を作る
a = np.random.randint(0, 10, (2, 1, 3))
b = np.random.randint(0, 10, (3, 1))

print('a:\n', a)
print('\na.shape:', a.shape)
print('\nb:\n', b)
print('\nb.shape:', b.shape)

# 加算
c = a + b

print('\na + b:\n', c)
print('\n(a + b).shape:', c.shape)

a:
 [[[6 3 8]]

 [[8 4 4]]]

a.shape: (2, 1, 3)

b:
 [[6]
 [7]
 [5]]

b.shape: (3, 1)

a + b:
 [[[12  9 14]
  [13 10 15]
  [11  8 13]]

 [[14 10 10]
  [15 11 11]
  [13  9  9]]]

(a + b).shape: (2, 3, 3)

a の形は (2, 1, 3) で、b の形は (3, 1) でした。この 2 つの配列の末尾次元 (trailing dimension) はそれぞれ 3 と 1 なので、ルールにあった「次元が同じサイズであるか、どちらかが 1 であること」を満たしています。

次に、各配列の第 2 次元に注目してみましょう。それぞれ 1 と 3 です。これもルールを満たしています。

ここで、a は 3 次元配列ですが、b は 2 次元配列です。つまり、次元数が異なっています。このような場合は、b は一番上の次元にサイズが 1 の次元が追加された形 (1, 3, 1) として扱われます。そして 2 つの配列の各次元ごとのサイズの最大値をとった形 (2, 3, 3) にブロードキャストされ、足し算が行われます。

このように、もし 2 つの配列のランクが異なる場合は、次元数が小さい方の配列が大きい方と同じ次元数になるまでその形の先頭に新たな次元が追加されます。サイズが 1 の次元がいくつ追加されても、要素の数は変わらないことに注意してください。要素数（size 属性で取得できる値）は、各次元のサイズの掛け算になるので、1 を何度かけても値は変わらないことから、これが成り立つことが分かります。

NumPy がブロードキャストのために自動的に行う新しい次元の挿入は、[] を使った以下の表な表記を用いることで手動で行うこともできます。

print('Original shape:', b.shape)

b_expanded = b[np.newaxis, :, :]

print('Added new axis to the top:', b_expanded.shape)

b_expanded2 = b[:, np.newaxis, :]

print('Added new axis to the middle:', b_expanded2.shape)

Original shape: (3, 1)
Added new axis to the top: (1, 3, 1)
Added new axis to the middle: (3, 1, 1)

np.newaxis が指定された位置に、新しい次元が挿入されます。配列が持つ数値の数は変わっていません。そのため、挿入された次元のサイズは必ず 1 になります。

array([[6],
       [7],
       [5]])

b_expanded

array([[[6],
        [7],
        [5]]])

b_expanded2

array([[[6]],

       [[7]],

       [[5]]])

NumPy のブロードキャストは慣れるまで直感に反するように感じる場合があるかもしれません。しかし、使いこなすと同じ計算が Python のループを使って行うよりも高速に行えるため、ブロードキャストを理解することは非常に重要です。一つ具体例を見てみます。

\(5 \times 5\) 行列 a に、3 次元ベクトル b を足します。まず、a、b および結果を格納する配列 c を定義します。

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

b = np.array([1, 2, 3, 4, 5])

# 結果を格納する配列を先に作る
c = np.empty((5, 5))

%%timeit という Jupyter Notebook で使用できるそのセルの実行時間を計測するためのマジックを使って、a の各行（1 次元目）に b の値を足していく計算を Python のループを使って 1 行ずつ処理していくコードの実行時間を測ってみます。

%%timeit
for i in range(a.shape[0]):
    c[i, :] = a[i, :] + b

7.35 µs ± 172 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)

array([[ 1.,  3.,  5.,  5.,  5.],
       [ 4.,  6.,  8.,  8.,  8.],
       [ 7.,  9., 11., 11., 11.],
       [ 4.,  6.,  8.,  8.,  9.],
       [ 1.,  3.,  5.,  5.,  5.]])

次に、NumPy のブロードキャストを活用した方法で同じ計算を行ってみます。

%%timeit
c = a + b

1.36 µs ± 24.7 ns per loop (mean ± std. dev. of 7 runs, 1,000,000 loops each)

array([[ 1.,  3.,  5.,  5.,  5.],
       [ 4.,  6.,  8.,  8.,  8.],
       [ 7.,  9., 11., 11., 11.],
       [ 4.,  6.,  8.,  8.,  9.],
       [ 1.,  3.,  5.,  5.,  5.]])

計算結果は当然同じになります。しかし、実行時間が数倍短くなっています。

このように、ブロードキャストを理解して活用することで、記述が簡単になるだけでなく、実行速度という点においても有利になります。

行列積#

行列の要素ごとの積は * を用いて計算できました。一方、通常の行列同士の積（行列積）の計算は、* ではなく、別の方法で行います。方法は 2 種類あります。

1つは、np.dot() 関数を用いる方法です。 np.dot() は 2 つの引数をとり、それらの行列積を計算して返す関数です。今、A という行列と B という行列があり、行列積 AB を計算したいとします。これは np.dot(A, B) と書くことで計算できます。もし BA を計算したい場合は、np.dot(B, A) と書きます。

もう 1 つは、ndarray オブジェクトが持つ dot() メソッドを使う方法です。これを用いると、同じ計算が A.dot(B) と書くことによって行えます。

# 行列 A の定義
A = np.array([
    [0, 1, 2],
    [3, 4, 5],
    [6, 7, 8]
])

# 行列 B の定義
B = np.array([
    [1, 2, 3],
    [4, 5, 6],
    [7, 8, 9]
])

実際にこの \(3 \times 3\) の 2 つの行列の行列積を計算してみましょう。

# 行列積の計算 (1)
C = np.dot(A, B)

C

array([[ 18,  21,  24],
       [ 54,  66,  78],
       [ 90, 111, 132]])

同じ計算をもう一つの記述方法で行ってみます。

C = A.dot(B)

C

array([[ 18,  21,  24],
       [ 54,  66,  78],
       [ 90, 111, 132]])

# データ型の確認（整数値）
a.dtype

dtype('int64')

基本的な統計量の求め方#

本節では、多次元配列に含まれる値の平均・分散・標準偏差・最大値・最小値といった統計値を計算する方法を紹介します。 \(8 \times 10\) の行列を作成し、この中に含まれる値全体に渡るこれらの統計値を計算してみましょう。

x = np.random.randint(0, 10, (8, 10))

x

array([[0, 6, 2, 1, 9, 7, 1, 4, 3, 7],
       [5, 3, 8, 9, 6, 7, 4, 0, 4, 6],
       [0, 9, 1, 1, 3, 1, 5, 0, 4, 7],
       [1, 7, 1, 8, 8, 2, 4, 6, 7, 7],
       [1, 8, 2, 1, 5, 3, 7, 1, 5, 4],
       [9, 0, 1, 8, 9, 5, 0, 7, 3, 3],
       [2, 9, 0, 9, 1, 1, 2, 0, 5, 0],
       [0, 9, 8, 9, 3, 7, 3, 2, 2, 5]])

# 平均値
x.mean()

4.1625

# 分散
x.var()

9.31109375

# 標準偏差
x.std()

3.0514084862567974

# 最大値
x.max()

# 最小値
x.min()

ここで、x は 2 次元配列なので、各次元に沿ったこれらの統計値の計算も行えます。例えば、最後の次元内だけで平均をとると、8 個の平均値が得られるはずです。平均を計算したい軸（何次元目に沿って計算するか）を axis という引数に指定します。

x.mean(axis=1)

array([4. , 5.2, 3.1, 5.1, 3.7, 4.5, 2.9, 4.8])

これは、以下のように 1 次元目の値の平均を計算していったものを並べているのと同じことです。（ゼロベースインデックスで考えています。x の形は (8, 10) なので、0 次元目のサイズが 8、1 次元目のサイズが 10 です。）

np.array([
    x[0, :].mean(),
    x[1, :].mean(),
    x[2, :].mean(),
    x[3, :].mean(),
    x[4, :].mean(),
    x[5, :].mean(),
    x[6, :].mean(),
    x[7, :].mean(),
])

array([4. , 5.2, 3.1, 5.1, 3.7, 4.5, 2.9, 4.8])

NumPy を用いた重回帰分析#

単回帰分析と重回帰分析で解説されている重回帰分析を NumPy を用いて行いましょう。

4 つのデータをまとめた、以下のようなデザイン行列が与えられたとします。

# Xの定義
X = np.array([
    [2, 3],
    [2, 5],
    [3, 4],
    [5, 9],
])

X

array([[2, 3],
       [2, 5],
       [3, 4],
       [5, 9]])

4 章の解説と同様に、切片を重みベクトルに含めて扱うため、デザイン行列の 0 列目に 1 という値を付け加えます。

# データ数（X.shape[0]) と同じ数だけ 1 が並んだ配列
ones = np.ones((X.shape[0], 1))

# concatenate を使い、1 次元目に 1 を付け加える
X = np.concatenate((ones, X), axis=1)

# 先頭に 1 が付け加わったデザイン行列
X

array([[1., 2., 3.],
       [1., 2., 5.],
       [1., 3., 4.],
       [1., 5., 9.]])

また、目標値が以下で与えられたとします。

# t の定義
t = np.array([1, 5, 6, 8])

t

array([1, 5, 6, 8])

重回帰分析は、正規方程式を解くことで最適な 1 次方程式の重みを決定することができました。正規方程式の解は以下のようなものでした。

\[ {\bf w} = ({\bf X}^{{\rm T}}{\bf X})^{\rm -1}{\bf X}^{\rm T}{\bf t} \]

これを、4 つのステップに分けて計算していきます。

まずは、\({\bf X}^{\rm T}{\bf X}\) の計算です。ndarrayに対して .T で転置した配列を得られます。

# Step 1
xx = np.dot(X.T, X)

xx

array([[  4.,  12.,  21.],
       [ 12.,  42.,  73.],
       [ 21.,  73., 131.]])

次に、この逆行列を計算します。

# Step 2
xx_inv = np.linalg.inv(xx)

xx_inv

array([[ 1.76530612, -0.39795918, -0.06122449],
       [-0.39795918,  0.84693878, -0.40816327],
       [-0.06122449, -0.40816327,  0.24489796]])

逆行列の計算は np.linalg.inv() で行うことができます。

次に、\({\bf X}^{\rm T}{\bf t}\) の計算をします。

# Step 3
xt = np.dot(X.T, t)

xt

array([ 20.,  70., 124.])

最後に、求めた xx_inv と xt を掛け合わせます。

# Step 4
w = np.dot(xx_inv, xt)

w

array([-0.14285714,  0.71428571,  0.57142857])

以上の計算は、行列積の演算子 @ を用いて 1 行で書くこともできます。

w_ = np.linalg.inv(X.T @ X) @ X.T @ t

w_

array([-0.14285714,  0.71428571,  0.57142857])

実際には逆行列を陽に求めることは稀で、連立一次方程式を解く、すなわち逆行列を計算してベクトルに掛けるのに等しい計算をひとまとめに行う関数 numpy.linalg.solve を呼ぶ方が速度面でも精度面でも有利です。

w_ = np.linalg.solve(X.T.dot(X), X.T.dot(t))

w_

array([-0.14285714,  0.71428571,  0.57142857])

数式を NumPy による配列の計算に落とし込むことに慣れていくには少し時間がかかりますが、慣れると少ない量のコードで記述できるだけでなく、高速に計算が行なえるため、大きな恩恵があります。