Kei Minagawa's Blog

皆川圭(@keimina)のブログ、Pythonで試したことを書いていきます

MacBook Pro の Radeon GPU で Keras を使って MNIST を高速に学習させる方法

2019年12月現在、MacBook ProRadeon GPU で Keras を使って MNIST を高速に学習させることができました。忘れないよう環境構築手順等を記載します。Plaid-ML という機械学習ライブラリをインストールし Keras のバックエンドにそれを指定するとできることを確認しました。

1 参考文献

以下の記事を参考にしています。

  • Tensorflowのdockerを使ってみる(macOS

qiita.com

  • PlaidML Kerasでやっていく #TokyoR 73

www.slideshare.net

plaidml.github.io

Plaid ML の詳細はこちら

github.com

2 環境構築

環境構築手順について記載します。 Plaid ML Installation Instructions macOS に記載の手順の通りです。

2.1 前提条件

実際に試した環境は Anaconda の Python3.7 が入っているバージョンをインストールした Macbook Pro です。以下の手順ではコマンドが ~/anaconda3/bin/pip のように長くなっていますが、これは私のPCで Anaconda へのパスを設定していない状態のためです。

2.2 "plaidml-keras" のインストール

以下のコマンドを実行します。

~/anaconda3/bin/pip install -U plaidml-keras
Collecting plaidml-keras
  Downloading https://files.pythonhosted.org/packages/49/a4/cb4c18eb0d5ec72e69b0cda53a8dcf9310d02060d2398672d0ac52deb394/plaidml_keras-0.6.4-py2.py3-none-any.whl
Collecting plaidml (from plaidml-keras)
  Downloading https://files.pythonhosted.org/packages/b6/84/2ba8babe2f0a3eb5ec95f5ac569fa9b8299e76b4dcccf197bf4336032142/plaidml-0.6.4-py2.py3-none-macosx_10_10_x86_64.whl (31.0MB)
     |████████████████████████████████| 31.0MB 3.4MB/s 
Collecting keras==2.2.4 (from plaidml-keras)
  Downloading https://files.pythonhosted.org/packages/5e/10/aa32dad071ce52b5502266b5c659451cfd6ffcbf14e6c8c4f16c0ff5aaab/Keras-2.2.4-py2.py3-none-any.whl (312kB)
     |████████████████████████████████| 317kB 5.4MB/s 
Requirement already satisfied, skipping upgrade: six in ./anaconda3/lib/python3.7/site-packages (from plaidml-keras) (1.12.0)
Requirement already satisfied, skipping upgrade: numpy in ./anaconda3/lib/python3.7/site-packages (from plaidml->plaidml-keras) (1.17.2)
Requirement already satisfied, skipping upgrade: cffi in ./anaconda3/lib/python3.7/site-packages (from plaidml->plaidml-keras) (1.12.3)
Collecting enum34>=1.1.6 (from plaidml->plaidml-keras)
  Downloading https://files.pythonhosted.org/packages/af/42/cb9355df32c69b553e72a2e28daee25d1611d2c0d9c272aa1d34204205b2/enum34-1.1.6-py3-none-any.whl
Requirement already satisfied, skipping upgrade: h5py in ./anaconda3/lib/python3.7/site-packages (from keras==2.2.4->plaidml-keras) (2.9.0)
Collecting keras-preprocessing>=1.0.5 (from keras==2.2.4->plaidml-keras)
  Downloading https://files.pythonhosted.org/packages/28/6a/8c1f62c37212d9fc441a7e26736df51ce6f0e38455816445471f10da4f0a/Keras_Preprocessing-1.1.0-py2.py3-none-any.whl (41kB)
     |████████████████████████████████| 51kB 7.8MB/s 
Collecting keras-applications>=1.0.6 (from keras==2.2.4->plaidml-keras)
  Downloading https://files.pythonhosted.org/packages/71/e3/19762fdfc62877ae9102edf6342d71b28fbfd9dea3d2f96a882ce099b03f/Keras_Applications-1.0.8-py3-none-any.whl (50kB)
     |████████████████████████████████| 51kB 9.4MB/s 
Requirement already satisfied, skipping upgrade: pyyaml in ./anaconda3/lib/python3.7/site-packages (from keras==2.2.4->plaidml-keras) (5.1.2)
Requirement already satisfied, skipping upgrade: scipy>=0.14 in ./anaconda3/lib/python3.7/site-packages (from keras==2.2.4->plaidml-keras) (1.3.1)
Requirement already satisfied, skipping upgrade: pycparser in ./anaconda3/lib/python3.7/site-packages (from cffi->plaidml->plaidml-keras) (2.19)
Installing collected packages: enum34, plaidml, keras-preprocessing, keras-applications, keras, plaidml-keras
Successfully installed enum34-1.1.6 keras-2.2.4 keras-applications-1.0.8 keras-preprocessing-1.1.0 plaidml-0.6.4 plaidml-keras-0.6.4

2.3 初期設定コマンド

以下のコマンドを実行します。

~/anaconda3/bin/plaidml-setup
PlaidML Setup (0.6.4)

Thanks for using PlaidML!

Some Notes:
  * Bugs and other issues: https://github.com/plaidml/plaidml
  * Questions: https://stackoverflow.com/questions/tagged/plaidml
  * Say hello: https://groups.google.com/forum/#!forum/plaidml-dev
  * PlaidML is licensed under the Apache License 2.0
 

Default Config Devices:
   metal_intel(r)_uhd_graphics_630.0 : Intel(R) UHD Graphics 630 (Metal)
   metal_amd_radeon_pro_5500m.0 : AMD Radeon Pro 5500M (Metal)

Experimental Config Devices:
   llvm_cpu.0 : CPU (LLVM)
   opencl_amd_radeon_pro_5500m_compute_engine.0 : AMD AMD Radeon Pro 5500M Compute Engine (OpenCL)
   opencl_intel_uhd_graphics_630.0 : Intel Inc. Intel(R) UHD Graphics 630 (OpenCL)
   opencl_cpu.0 : Intel CPU (OpenCL)
   metal_intel(r)_uhd_graphics_630.0 : Intel(R) UHD Graphics 630 (Metal)
   metal_amd_radeon_pro_5500m.0 : AMD Radeon Pro 5500M (Metal)

Using experimental devices can cause poor performance, crashes, and other nastiness.

以下の画面で y を選択します。

Enable experimental device support? (y,n)[n]:y

以下の画面で 6 を選択します。

Multiple devices detected (You can override by setting PLAIDML_DEVICE_IDS).
Please choose a default device:

   1 : llvm_cpu.0
   2 : opencl_amd_radeon_pro_5500m_compute_engine.0
   3 : opencl_intel_uhd_graphics_630.0
   4 : opencl_cpu.0
   5 : metal_intel(r)_uhd_graphics_630.0
   6 : metal_amd_radeon_pro_5500m.0

Default device? (1,2,3,4,5,6)[1]:6

以下の画面で y を選択します。

Selected device:
    metal_amd_radeon_pro_5500m.0

Almost done. Multiplying some matrices...
Tile code:
  function (B[X,Z], C[Z,Y]) -> (A) { A[x,y : X,Y] = +(B[x,z] * C[z,y]); }
Whew. That worked.

Save settings to /Users/kei/.plaidml? (y,n)[y]:y

Success! が表示されることを確認します。

Success!

3 Keras で GPU を使用するために追加するコード

GPU を使用する場合は以下のように、コードを先頭に追加します。

import os
os.environ["KERAS_BACKEND"] = "plaidml.keras.backend"
import keras

4 MNIST の学習

動作確認として Keras 公式サイトの MNIST のニューラルネットを学習させます。

https://keras.io/examples/mnist_cnn/

上記サイトに記載のコードの一番上に、 Keras で GPU を使用するために追加するコード を追加しただけです。ソースコード全体は以下になります。

import os
os.environ["KERAS_BACKEND"] = "plaidml.keras.backend"
import keras
from keras.datasets import mnist
from keras.models import Sequential
from keras.layers import Dense, Dropout, Flatten
from keras.layers import Conv2D, MaxPooling2D
from keras import backend as K

batch_size = 128
num_classes = 10
epochs = 12

# input image dimensions
img_rows, img_cols = 28, 28

# the data, split between train and test sets
(x_train, y_train), (x_test, y_test) = mnist.load_data()

if K.image_data_format() == 'channels_first':
        x_train = x_train.reshape(x_train.shape[0], 1, img_rows, img_cols)
        x_test = x_test.reshape(x_test.shape[0], 1, img_rows, img_cols)
        input_shape = (1, img_rows, img_cols)
else:
        x_train = x_train.reshape(x_train.shape[0], img_rows, img_cols, 1)
        x_test = x_test.reshape(x_test.shape[0], img_rows, img_cols, 1)
        input_shape = (img_rows, img_cols, 1)

x_train = x_train.astype('float32')
x_test = x_test.astype('float32')
x_train /= 255
x_test /= 255
print('x_train shape:', x_train.shape)
print(x_train.shape[0], 'train samples')
print(x_test.shape[0], 'test samples')

# convert class vectors to binary class matrices
y_train = keras.utils.to_categorical(y_train, num_classes)
y_test = keras.utils.to_categorical(y_test, num_classes)

model = Sequential()
model.add(Conv2D(32, kernel_size=(3, 3),
                                 activation='relu',
                                 input_shape=input_shape))
model.add(Conv2D(64, (3, 3), activation='relu'))
model.add(MaxPooling2D(pool_size=(2, 2)))
model.add(Dropout(0.25))
model.add(Flatten())
model.add(Dense(128, activation='relu'))
model.add(Dropout(0.5))
model.add(Dense(num_classes, activation='softmax'))

model.compile(loss=keras.losses.categorical_crossentropy,
                          optimizer=keras.optimizers.Adadelta(),
                          metrics=['accuracy'])

model.fit(x_train, y_train,
                  batch_size=batch_size,
                  epochs=epochs,
                  verbose=1,
                  validation_data=(x_test, y_test))
score = model.evaluate(x_test, y_test, verbose=0)
print('Test loss:', score[0])
print('Test accuracy:', score[1])

出力は以下になります。処理速度は 15秒/Epoch くらいです。

Downloading data from https://s3.amazonaws.com/img-datasets/mnist.npz
11493376/11490434 [==============================] - 5s 0us/step
x_train shape: (60000, 28, 28, 1)
60000 train samples
10000 test samples
INFO:plaidml:Opening device "metal_amd_radeon_pro_5500m.0"
Train on 60000 samples, validate on 10000 samples
Epoch 1/12
60000/60000 [==============================] - 17s 282us/step - loss: 0.2692 - acc: 0.9164 - val_loss: 0.0545 - val_acc: 0.9825
Epoch 2/12
60000/60000 [==============================] - 13s 219us/step - loss: 0.0891 - acc: 0.9738 - val_loss: 0.0398 - val_acc: 0.9864
Epoch 3/12
60000/60000 [==============================] - 13s 219us/step - loss: 0.0648 - acc: 0.9807 - val_loss: 0.0351 - val_acc: 0.9886
Epoch 4/12
60000/60000 [==============================] - 13s 219us/step - loss: 0.0549 - acc: 0.9836 - val_loss: 0.0340 - val_acc: 0.9878
Epoch 5/12
60000/60000 [==============================] - 13s 217us/step - loss: 0.0469 - acc: 0.9859 - val_loss: 0.0317 - val_acc: 0.9891
Epoch 6/12
60000/60000 [==============================] - 13s 219us/step - loss: 0.0405 - acc: 0.9875 - val_loss: 0.0298 - val_acc: 0.9904
Epoch 7/12
60000/60000 [==============================] - 13s 219us/step - loss: 0.0364 - acc: 0.9887 - val_loss: 0.0276 - val_acc: 0.9900
Epoch 8/12
60000/60000 [==============================] - 13s 216us/step - loss: 0.0338 - acc: 0.9899 - val_loss: 0.0314 - val_acc: 0.9898
Epoch 9/12
60000/60000 [==============================] - 13s 216us/step - loss: 0.0299 - acc: 0.9908 - val_loss: 0.0272 - val_acc: 0.9913
Epoch 10/12
60000/60000 [==============================] - 13s 215us/step - loss: 0.0276 - acc: 0.9912 - val_loss: 0.0263 - val_acc: 0.9914
Epoch 11/12
60000/60000 [==============================] - 13s 218us/step - loss: 0.0266 - acc: 0.9917 - val_loss: 0.0268 - val_acc: 0.9910
Epoch 12/12
60000/60000 [==============================] - 13s 219us/step - loss: 0.0256 - acc: 0.9920 - val_loss: 0.0272 - val_acc: 0.9922
Test loss: 0.02717138020992279
Test accuracy: 0.9922

5 まとめ

Plaid ML をインストール、初期設定を行いました。Keras のバックエンドとして Plaid ML を指定することで Macbook ProRadeon GPU を使用することができました。

以上です。おやすみなさい。

ボロノイ図を書いて、マウスクリックした時に一番近いデータポイントをを表示するプログラムを実装する

ボロノイ図を書くプログラムと、マウスクリックした時に一番近いデータポイントをを表示するプログラムを実装します。

1 ボロノイ図とは

ボロノイ図の詳細は Wikipedia を参照してください。

ja.wikipedia.org

簡単な例をあげると、xy平面上に点Aと点Bがあったとして、点Xが点Aに近ければ赤、点Bに近ければ青というように、点Xに色を塗るとします。xy平面上にある点全ての点Xに前者の操作を繰り返していくと、全ての領域が赤または青に塗られると思います。この例では、点Aと点Bの二つしかありませんでしたが、多数の点でこのようなことを行うと、色のついた領域が四角形や六角形などの多角形になり、独特な図になります。このような図形をボロノイ図と呼ぶと認識しています。

2 ボロノイ図を書くプログラム実装する

2.1 コード

ボロノイ図を書くコードを実装します。コードは以下の通りです。

# -*-coding:utf-8 -*-
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

np.random.seed(0)

# グリッドに関する設定
GRID_N = 100
GRID_M= 100

STEP_x = 1.0/GRID_N
STEP_y = 1.0/GRID_M
mx, my = np.meshgrid(np.linspace(0, 1, GRID_N), np.linspace(0, 1, GRID_M))

# データに関する設定
DATA_N = 10
df = pd.DataFrame(np.random.random((DATA_N, 2)), columns=['x', 'y'])

# データをブロードキャスト(データを3階のテンソルへ変換)
x_broad = np.broadcast_to(df.x[:, np.newaxis, np.newaxis], (DATA_N, GRID_N, GRID_N))
y_broad = np.broadcast_to(df.y[:, np.newaxis, np.newaxis], (DATA_N, GRID_M, GRID_N))

# 2乗距離が最小の平面のインデックスを求める
dx = x_broad - mx
dy = y_broad - my

r = dx**2 + dy**2
z = r.argmin(axis=0)

# 描画
fig, ax = plt.subplots()
ax.plot(df.x, df.y, lw=0, marker='o')
ax.pcolormesh(mx, my, z)

# アスペクト比を 1:1 にする
ax.set_aspect('equal')
plt.show(True)

2.2 コードの説明

コードの説明は以下の通りです。

説明
21行目
x_broad = np.broadcast_to(df.x[:, np.newaxis, np.newaxis], (DATA_N, GRID_N, GRID_N))
np.broadcast_to 関数を使用して、配列の 1D -> 3D へのブロードキャストを行っています。np.broadcast_to 関数については、以下を参照してください。

keimina.hatenablog.jp

25
dx = x_broad - mx
np.meshgrid で作成した行列 x からデータのポイントをブロードキャストして作成した配列を引き算しています。 3次元と2次元の計算ですので、これもブロードキャストの演算が行われます。イメージ的には0軸が上からしたに向かって伸びており、そこから横に、枝のように、1軸と2軸からなる平面がぶら下がってる状態で、そのぶら下がってる平面それぞれに対して一つ、一つに行列 mx の平面と引き算を行うという感じです。何行ってるかわかんないですね。
28
r = dx**2 + dy**2
3次元の r 配列の要素である行列はデータポイントとの距離を表しています。
29
z = r.argmin(axis=0)
3次元の r 配列の0軸方向で一番値の小さいインデックスを求めています。これは、xy座標上のある領域において、その領域では、どのデータポイントが一番近いか求めているということです。
33
ax.plot(df.x, df.y, lw=0, marker='o')
plot 関数により、データポイントの描画を行います。 scatter でも良いでしょう。
34
ax.pcolormesh(mx, my, z)
pcolormesh 関数により、領域に色を付けます

3 出力結果

出力結果は以下になります。青い点がある領域内における点に一番近い青い点は、その領域内にある青い点になっていることがわかります。

f:id:keimina:20191215203434p:plain

図1. 出力結果(ボロノイ図)

4 マウスクリックしたときに一番近いデータポイントを表示するプログラムを実装する

Matplotlib で何らかのデータの点の分布を描画したときに、一番近いデータポイントを表示したい時があります。これは、ボロノイ図アルゴリズムを応用し簡単に実現できます。以下にコードだけ貼り付けておきます。

4.1 コード

# -*-coding:utf-8 -*-
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

np.random.seed(0)

# グリッドに関する設定
GRID_M = 100
GRID_N = 100

STEP_x = 1.0/GRID_N
STEP_y = 1.0/GRID_M
mx, my = np.meshgrid(np.linspace(0, 1, GRID_N), np.linspace(0, 1, GRID_M))

# データに関する設定
DATA_N = 10
df = pd.DataFrame(np.random.random((DATA_N, 2)), columns=['x', 'y'])

print(df)

# データをブロードキャスト(データを3階のテンソルへ変換)
x_broad = np.broadcast_to(df.x[:, np.newaxis, np.newaxis], (DATA_N, GRID_M, GRID_N))
y_broad = np.broadcast_to(df.y[:, np.newaxis, np.newaxis], (DATA_N, GRID_M, GRID_N))

# 2乗距離が最小の平面のインデックスを求める
dx = x_broad - mx
dy = y_broad - my

r = dx**2 + dy**2
z = r.argmin(axis=0)

# 描画
fig, ax = plt.subplots()
ax.plot(df.x, df.y, lw=0, marker='o')
ax.pcolormesh(mx, my, z)

def onclick(event):
        if (event.xdata is None) or (event.ydata is None):
                return
        gx = int(event.xdata*GRID_N)
        gy = int(event.ydata*GRID_M)
        index = int(z[gy, gx])
        print(df.loc[[index]])

cid = fig.canvas.mpl_connect('button_press_event', onclick)

# アスペクト比を 1:1 にする
ax.set_aspect('equal')
ax.set_xlim(0,1)
ax.set_ylim(0,1)
plt.show(True)

f:id:keimina:20191215203255g:plain

図2. マウスクリックした座標に一番近いデータポイントの座標を表示するプログラム

5 まとめ

ボロノイ図を書くプログラムを実装しました。またその応用としてマウスクリックした座標に一番近いデータポイントの座標を表示するプログラムを作りました。

以上です。

Author: MINAGAWA KEI

Created: 2019-12-15 日 20:24

Validate

`np.broadcast_to`関数を使用して1次元配列を3次元配列に変換する

今回は、1次元配列を3次元配列に変換する方法について書きます。実装には、 `np.broadcast_to` 関数を使用しました。

1 今回、作成したい配列

どのような配列の変換をしたいかというと、変換後の配列の各要素に、変換前の要素の値がいい感じに入っている配列です。"いい感じ"が言葉では伝わらないと思いますので、作成したい配列について絵にしました。同じ色の領域は同じ値が入っています。このように、値を保ちつつ自然に1次元の配列を3次元に拡張するのが今回の目的です。

f:id:keimina:20191215135557j:plain

図1. 作成したい配列

2 作成する方法 その1

作成する方法の一つを書きます。以下の手順で変換するのが一般的かと思います。

  1. 1次元配列の一つの要素を取り出して、その取り出した値を要素とする行列を作成する
  2. 作成した行列をリストにappend
  3. 1次元配列の全ての要素について 上記(1)(2)を実行する
  4. (2)のリストが変換結果です

3 作成する方法 その2

今回の配列では、作成する方法はもう一つあります。以下の、方法です。

  1. np.broadcast_to 関数を使用する

今回は、こちらで実装します。もちろん、作成する方法 その1 を、そのまま Python で実装しても良いです(むしろそっちの方が可読性が良いかも)。

4 コード

コードは以下の通りです。

# -*-coding:utf-8 -*-
import numpy as np
M, N = 4, 4
# 入力
arr1d = np.array([1, 2, 3, 4])
# 出力
arr3d = np.broadcast_to(arr1d[:,np.newaxis,np.newaxis], (arr1d.shape[0], M, N))
print(arr3d)

5 出力結果

出力結果は以下の通りです。arr3d[0] が値が1の要素からなる行列になっています。arr3d[1], arr3d[2], arr3d[3] も同様です。意図した配列が作成されていることがわかります。

array([[[1, 1, 1, 1],
        [1, 1, 1, 1],
        [1, 1, 1, 1],
        [1, 1, 1, 1]],

       [[2, 2, 2, 2],
        [2, 2, 2, 2],
        [2, 2, 2, 2],
        [2, 2, 2, 2]],

       [[3, 3, 3, 3],
        [3, 3, 3, 3],
        [3, 3, 3, 3],
        [3, 3, 3, 3]],

       [[4, 4, 4, 4],
        [4, 4, 4, 4],
        [4, 4, 4, 4],
        [4, 4, 4, 4]]])

6 コードの説明

`np.broadcast_to` 関数の第1引数に与えている `arr1d[:,np.newaxis,np.newaxis]` ですが、これは arr1d を配列を要素数をそのままに次元数を3次元にを増やした配列を表しています。np.newaxis は None で定義されていますが、このように使用すると、配列に新しい軸を与え、次元を増やすことができます。第2引数は、変換したい配列の shape を与えています。

7 まとめ

1次元配列を3次元配列に変換する方法、np.broadcast_to 関数について書きました。

以上です。

Author: MINAGAWA KEI

Created: 2019-12-15 日 13:43

Validate

pd.read_csv でデータを読み込んだ時の列の型

自分が主催をしている Pandas 勉強会で read_csv について学んだことの一部をここに記載します。結論を言ってしまうと、 csv を受け取る時は列の型を聞いておく必要があるよねという話です。最後に、勉強会で、 pd.read_csv の引数である keep_date_col についても勉強したのでそれも記載します。

本題に入る前に ー テキストを csv ファイル形式で扱う方法

本題に入る前にテキストを csv ファイル形式で扱う方法について記載します。これは Pandas ではなく標準 Python のモジュールの話です。csv ファイルの読み込みの練習をする時には、ソースコード上に直接記載されている csv を仮想的にファイルとして扱えるようにできると便利です。以下のように io モジュールを使用するとメモリ上のデータをあたかもファイル上のデータであるかのように扱うことができます。(ただし、全く同じではありませんので注意してください)

import pandas as pd
import io

t = '''col1,col2,col3,col4,col5,col6
2019,10/1,1,2,3,NA
2020,11/1,,6,7,A
2021,12/1,9,10,11,X'''

f = io.StringIO(t)
print(f.read())
col1,col2,col3,col4,col5,col6
2019,10/1,1,2,3,NA
2020,11/1,,6,7,A
2021,12/1,9,10,11,X

Pandas で csv ファイルを読みこむ方法

ここから本題ですが、Pandas で csv ファイルを読みこむ方法は以下になります。pd.read_csv 関数でファイルを読み込みます。 変数 f は io.StringIO で作成したファイルオブジェクトのようなものです。

f = io.StringIO(t)
df = pd.read_csv(f)
print(df)
   col1  col2  col3  col4  col5 col6
0  2019  10/1   1.0     2     3  NaN
1  2020  11/1   NaN     6     7    A
2  2021  12/1   9.0    10    11    X

読み込んだ列の型(dtype)について

print 関数で df を 表示しました。ここで気になるのは、例えば 1列目が文字列なのか日付型なのか、はたまた整数型なのかということです。これを調べるために以下を実行します。 DataFrame の applymap メソッドを使用し全ての要素に type 関数を適用して型を調べています。

# 全ての要素の型を調べる
df.applymap(lambda x: type(x))
col1 col2 col3 col4 col5 col6
0 <class 'int'> <class 'str'> <class 'float'> <class 'int'> <class 'int'> <class 'float'>
1 <class 'int'> <class 'str'> <class 'float'> <class 'int'> <class 'int'> <class 'str'>
2 <class 'int'> <class 'str'> <class 'float'> <class 'int'> <class 'int'> <class 'str'>

実行結果を見ると、 col1 の先頭行の要素は 2019 となっており整数型になっていることがわかります。 col3 の要素は整数かと思いきや 浮動少数点型になっています(未入力の要素がありそれが NaN (浮動小数点型)になるため、その列の他の要素もそれになっているのかもしれません)。このように、 pd.read_csv 関数をデフォルトの状態で使用してデータを読み込んだ場合、想定する型と違うことになっていることがあるかもしれないということに注意が必要です。これは、csv は型の情報を含まないデータフォーマットであるためです。csv をちゃんと処理したい場合は、列の型の情報は知っていなければならないということです。csv のデータを他の人から受け取る時には列の型について聞いておくようにします。これは Pandas や Python に関係のない話です。

keep_date_col について

ここからの話は、おまけです。

勉強会では教科書(Python for Data Analysis)に記載されている引数について一つづつ確認して行ったのですが、ここでは、他の人が記事あまりにしていないもの keep_date_col について紹介しようと思います(需要はあまりないと思いますが)。以下のように parse_dates 引数と一緒に使用します(parse_dates引数に2次元リストを渡していることに注目)。 parse_dates 引数に2次元リストを渡すことで、日付列の集合を併合できます。この時、併合前の列を消す(デフォルトの動作)か保持するかを keep_date_col 引数で決めることができます。

keep_date_col=False (デフォルト)の時

f = io.StringIO(t)
df = pd.read_csv(f, parse_dates=[[0, 1]], keep_date_col=False)
df
col1_col2 col3 col4 col5 col6
0 2019-10-01 1.0 2 3 NaN
1 2020-11-01 NaN 6 7 A
2 2021-12-01 9.0 10 11 X

まず、 keep_date_col=False (デフォルト)の場合 parse_dates 引数で指定した列番号が日付パーサーで一つの列として併合されます。 parse_dates 引数で指定する列は複数列を指定する時は2次元配列として指定することに注意してください。

keep_date_col=True の時

f = io.StringIO(t)
df = pd.read_csv(f, parse_dates=[[0, 1]], keep_date_col=True)
df
col1_col2 col1 col2 col3 col4 col5 col6
0 2019-10-01 2019 10/1 1.0 2 3 NaN
1 2020-11-01 2020 11/1 NaN 6 7 A
2 2021-12-01 2021 12/1 9.0 10 11 X

keep_date_col=True の時は併合前の日付列が保持されていることがわかります。

date_parser で複数の列の集合を併合する

parse_dates 引数で指定する列は複数列を指定する時は2次元配列を指定するのでした。これは以下のように複数の列の集合を併合できるようしているからです。

まずテーブルを定義して、普通に読み込んだものを表示します。

t = '''col1,col2,col3,col4,col5,col6,col7,col8,col9,col10,col11,co12
2019,10/1,1.0,2,3,,2019,1/1,1.0,2,3,
2020,11/1,,6,7,A,2020,2/1,,6,7,A
2021,12/1,9.0,10,11,X,2021,3/1,9.0,10,11,X
'''

f = io.StringIO(t)
df = pd.read_csv(f)
df
col1 col2 col3 col4 col5 col6 col7 col8 col9 col10 col11 co12
0 2019 10/1 1.0 2 3 NaN 2019 1/1 1.0 2 3 NaN
1 2020 11/1 NaN 6 7 A 2020 2/1 NaN 6 7 A
2 2021 12/1 9.0 10 11 X 2021 3/1 9.0 10 11 X

次に、 parse_dates で複数の列の集合を指定してみます。

f = io.StringIO(t)
df = pd.read_csv(f, parse_dates=[[0, 1], [6, 7]], keep_date_col=False)
df
col1_col2 col7_col8 col3 col4 col5 col6 col9 col10 col11 co12
0 2019-10-01 2019-01-01 1.0 2 3 NaN 1.0 2 3 NaN
1 2020-11-01 2020-02-01 NaN 6 7 A NaN 6 7 A
2 2021-12-01 2021-03-01 9.0 10 11 X 9.0 10 11 X

parse_dates 引数で指定した列の集合それぞれに対して併合が行われていることがわかります。

以上です。おやすみなさい。

Pandas の Series オブジェクトと DataFrame オブジェクト の四則演算(足し算)した時の挙動についていろいろ試してみた

Pandas の Series オブジェクトと DataFrame オブジェクト の四則演算(足算)の挙動を調べてみました。結果、 Series オブジェクト の index のラベルと DataFrame オブジェクト の columns のラベルがマッチする時、計算が行われることがわかりました。また計算は numpy 風に行われることがわかりました。説明が難しいため、コードだけ貼り付けておきます。

import pandas as pd
import numpy as np

np.random.seed(777)
s_data  = np.random.randint(0,5,5)
s_index = list('ABCDE')
s = pd.Series(s_data, index=s_index)

s
# In [243]: s
# Out[243]: 
# A    3
# B    1
# C    4
# D    1
# E    2
# dtype: int64

d_data  = np.random.randint(0,10,(3, 5))
d_columns = list('ABCDE')
d = pd.DataFrame(d_data, columns=d_columns)

d
# In [244]: d
# Out[244]: 
#    A  B  C  D  E
# 0  8  7  2  0  1
# 1  2  4  5  7  1
# 2  7  2  2  7  4

s + d
# In [245]: s + d
# Out[249]: 
#     A  B  C  D  E
# 0  11  8  6  1  3
# 1   5  5  9  8  3
# 2  10  3  6  8  6

d + s
# In [250]: d + s
# Out[251]: 
#     A  B  C  D  E
# 0  11  8  6  1  3
# 1   5  5  9  8  3
# 2  10  3  6  8  6

# s の index を並べ替えたものを s2 とする
s_index = np.random.choice(s.index.values, replace=False, size=5)
s2 = s.reindex(s_index)

s2
# In [439]: s2
# Out[439]: 
# B    1
# E    2
# D    1
# C    4
# A    3

s + s2
# In [440]: s + s2
# Out[445]: 
# A    6
# B    2
# C    8
# D    2
# E    4

s2 + d
# In [446]: s2 + d
# Out[452]: 
#     A  B  C  D  E
# 0  11  8  6  1  3
# 1   5  5  9  8  3
# 2  10  3  6  8  6

s_data  = np.concatenate(([-100], s2.values, [100]))
s_index = np.concatenate((['X'] , s2.index , ['Y']))
s3 = pd.Series(s_data, s_index)

s3
# In [495]: s3
# Out[495]: 
# X   -100
# B      1
# E      2
# D      1
# C      4
# A      3
# Y    100

s + s3
# In [496]: s + s3
# Out[505]: 
# A    6.0
# B    2.0
# C    8.0
# D    2.0
# E    4.0
# X    NaN
# Y    NaN


d + s3
# In [506]: d + s3
# Out[514]: 
#     A  B  C  D  E   X   Y
# 0  11  8  6  1  3 NaN NaN
# 1   5  5  9  8  3 NaN NaN
# 2  10  3  6  8  6 NaN NaN

# 上記の ABCDE列 については計算されていることを確認
(d + s3).loc[:, list('ABCDE')] == d + s2
# In [590]: (d + s3).loc[:, list('ABCDE')] == d + s2
# Out[612]: 
#       A     B     C     D     E
# 0  True  True  True  True  True
# 1  True  True  True  True  True
# 2  True  True  True  True  True

#########################################################################
# A,B,C,D,E の列は共通のラベルなので NaN にはならないが、 X,Y は共通で
# はないので NaN になる。以下のように共通のラベルがない時は全て NaN に
# なる。
#########################################################################

s_data  = np.random.randint(0,5,5)
s_index = list('PQRST')
s4 = pd.Series(s_data, index=s_index)

s4
# In [581]: s4
# Out[583]: 
# P    4
# Q    0
# R    0
# S    0
# T    3

s + s4
# In [584]: s + s4
# Out[586]: 
# A   NaN
# B   NaN
# C   NaN
# D   NaN
# E   NaN
# P   NaN
# Q   NaN
# R   NaN
# S   NaN
# T   NaN

d + s4
# In [587]: d + s4
# Out[589]: 
#     A   B   C   D   E   P   Q   R   S   T
# 0 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
# 1 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
# 2 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN

# 行ラベルのみに Series の index と同じ ラベルがあった場合
d.T
# In [621]: d.T
# Out[621]: 
#    0  1  2
# A  8  2  7
# B  7  4  2
# C  2  5  2
# D  0  7  7
# E  1  1  4

d.T.index
# In [622]: d.T.index
# Out[631]: Index(['A', 'B', 'C', 'D', 'E'], dtype='object')

d.T.columns
# In [632]: d.T.columns
# Out[635]: RangeIndex(start=0, stop=3, step=1)

d + s
# In [649]: d + s
# Out[649]: 
#     A  B  C  D  E
# 0  11  8  6  1  3
# 1   5  5  9  8  3
# 2  10  3  6  8  6

d.T + s
# In [650]: d.T + s
# Out[652]: 
#     0   1   2   A   B   C   D   E
# A NaN NaN NaN NaN NaN NaN NaN NaN
# B NaN NaN NaN NaN NaN NaN NaN NaN
# C NaN NaN NaN NaN NaN NaN NaN NaN
# D NaN NaN NaN NaN NaN NaN NaN NaN
# E NaN NaN NaN NaN NaN NaN NaN NaN

# 上記よりデフォルトでは常に DataFrame の列ラベルに対してマッチングが行われることがわかる

####################################################################
# おまけ(1)
# DataFrame オブジェクト の columns ではなく、 index にマッチして計算
# して欲しい場合は、 メソッドを使用し引数に axis='index' を指定すると
# できる
####################################################################

d.T.add(s, axis='index')
# In [745]: d.T.add(s, axis='index')
# Out[746]: 
#     0  1   2
# A  11  5  10
# B   8  5   3
# C   6  9   6
# D   1  8   8
# E   3  3   6

d.T.add(s, axis='index').T == d + s
# In [776]: d.T.add(s, axis='index').T == d + s
# Out[781]: 
#       A     B     C     D     E
# 0  True  True  True  True  True
# 1  True  True  True  True  True
# 2  True  True  True  True  True

####################################################################
# おまけ(2)
# index などを無視して、何がなんでも足したい場合は、 values に対して計
# 算を行う(この場合 numpy の計算なので、 numpy broadcastng rule が適用
# 可能な形状にかぎる))
####################################################################

s4
# In [782]: s4
# Out[814]: 
# P    4
# Q    0
# R    0
# S    0
# T    3

d
# In [815]: d
# Out[823]: 
#    A  B  C  D  E
# 0  8  7  2  0  1
# 1  2  4  5  7  1
# 2  7  2  2  7  4

dd = d.copy()
dd.loc[:,:] = dd.values + s4.values
# In [931]: dd
# Out[931]: 
#     A  B  C  D  E
# 0  12  7  2  0  4
# 1   6  4  5  7  4
# 2  11  2  2  7  7

共通ラベルにない値は 「NaN + 値」 で 「NaN」になってしまうが、 「NaN + 値」を 「0 + 値」、 「NaN + NaN」 を 「0 + 0」 のように計算して欲しい場合は、以下のように fill_value 引数を使用する。(ただし、DataFrame と Series 同士の演算で fill_value を使用するとエラーになってしまったので注意)

s
# In [243]: s
# Out[243]: 
# A    3
# B    1
# C    4
# D    1
# E    2
# dtype: int64

s3
# In [495]: s3
# Out[495]: 
# X   -100
# B      1
# E      2
# D      1
# C      4
# A      3
# Y    100

s + s3
# In [496]: s + s3
# Out[505]: 
# A    6.0
# B    2.0
# C    8.0
# D    2.0
# E    4.0
# X    NaN
# Y    NaN

s.add(s3, fill_value=0)
# In [938]: s.add(s3, fill_value=0)
# Out[971]: 
# A      6.0
# B      2.0
# C      8.0
# D      2.0
# E      4.0
# X   -100.0
# Y    100.0


d + s3
# In [506]: d + s3
# Out[514]: 
#     A  B  C  D  E   X   Y
# 0  11  8  6  1  3 NaN NaN
# 1   5  5  9  8  3 NaN NaN
# 2  10  3  6  8  6 NaN NaN

d.add(s3, fill_value=0.0)
# NotImplementedError: fill_value 0.0 not supported.
# 上のようにDataFrame と Series で fill_value 引数を指定するとエラーになる
# パッと思いつたところでは以下のようにかけば実現できるのではないかと思います
d.reindex(columns=s3.index).fillna(0) + s3.fillna(0)

# ただこれだと、元々持っていた d の列が消えてしまうのでそれを保持するには
# 以下のようにするとできるのではないかと思います
dd = d.reindex(columns=d.columns|s3.index)
dd.loc[:, s3.index] = d.reindex(columns=s3.index).fillna(0) + s3.fillna(0)

以上です。おやすみなさい。

Pandas の obj[i] と obj.loc[i] などの違いについて

Pandas の Series オブジェクトの話になります。 Series オブジェクトを obj とした時、 obj[i] と obj.loc[i] と obj.iloc[i] の違い、さらに obj[s:e] のようにスライスした時の動作の違いを理解するために、コードを書き実行して動作を確かめました。確認のためのコードは以下になります。Series オブジェクトの index が整数のみからなるか否かで、 obj[i] の振る舞いが変わるのには注意が必要だと思いました。コードの説明はコメントとして記載しました。

############################################################################
# index が整数型 Seires の場合(index の値が整数のみからなる Series の場合)
############################################################################
import pandas as pd

s1 = pd.Series([15,20,25,30,40], index=[3,4,5,6,7])
s2 = pd.Series([15,20,25,30,40])

# (1) s1[-1] を取得する
s1[-1]
# KeyError: -1

# (2) s2[-1] を取得する
s2[-1]
# KeyError: -1

# 上記よりindex が整数のみからなる場合は、obj[i] の i はラベルのように認
# 識されることがわかる。
# ただし、スライスの場合は obj[s:e] の s, e はインデックスとして扱われる

# (3) s1[3:5] を取得する
s1[3:5]
# 6    30
# 7    40
# dtype: int64

# (4) s1.loc[-1] を取得する
s1.loc[-1]
# KeyError: 'the label [-1] is not in the [index]'

# (5) s1.loc[3:5] で loc にスライスを使用して取得する
s1.loc[3:5]
# 3    15
# 4    20
# 5    25
# dtype: int64

# (6) s1.iloc[-1] を使った場合
s1.iloc[-1]
# 40

# (7) s1.iloc[3:5] で iloc にスライスを使用して取得する
s1.iloc[3:5]
# 6    30
# 7    40
# dtype: int64

################################################################
# index がオブジェクト型の  Series の場合
################################################################

s3 = pd.Series([15,20,25,30,40], index=[3,4,5,6,'X'])

# (8) s3[-1] を取得する
s3[-1]
# 40

# (9) s3[3] を取得する
s3[3]
# 15

# index がオブジェクト型の場合は、 obj[i] はインデックスとして扱われる
# loc, iloc, スライス を使用した際の挙動は (3)(4)(5)(6)(7) と同様になると思われる(未確認)

################################################################
# index が重複している時の Series の場合
################################################################

# インデックスが重複しているとき
s4 = pd.Series([10, 100, 1000, 10000], index=['a', 'a', 'b', 'b'])
s5 = pd.Series([200, 2000], index=['a', 'b'])
s6 = pd.Series([200, 200], index=['a', 'a'])

# (10) 値を取得する
s4['a']
# a     10
# a    100
# dtype: int64

# (11) 値を取得する
s4[['a','b']]
# a       10
# a      100
# b     1000
# b    10000
# dtype: int64

# (12) 重複していないインデックスを重複するように reindex する
s5.reindex(s4.index)
# a     200
# a     200
# b    2000
# b    2000
# dtype: int64

# (13) インデックスが重複しているとき reindex する
s6.reindex(['x', 'y'])
# ValueError: cannot reindex from a duplicate axis

以上です。

三面体のサイコロを2000回投げたとき、それぞれの面の出る回数の分布

1. はじめに

(注意:一般的に三面体のサイコロというものは現実には存在しません。説明の便宜上、3つの面をもつ空想のサイコロのことを三面体のサイコロと表現しています)
Pandas 勉強会で、サンプルの状態が割とよく発生することなのか、稀にしか発生しないことなのかを知りたいと言う話がありました。それに関連するような、確率の問題を自分で作って、シミュレーションを行いました。結果、仮定がデータと紐づいているのであれば、割とよく発生することなのか、稀にしか発生しないことなのかがわかることがわかりました。一方で仮定がデータと紐づかないのであれば、仮定に基づいて議論を進めたところであまり意味のない議論になってしまうこともわかりました。

2. 問題

以下のように問題を定義しました。

2000個のサンプルがあったとして、 3つの領域 A, B, C の内どれかに分類されるとした時、観測したそのサンプルが領域 A に分類された個数、 B に分類された個数、 C に分類された個数の状態が、割とよくおこるものなのか、そうでないのか知りたい。

3. 問題を解くための仮定

まず、Aに分類されるか、B または C に分類されるかは一定の確率で決まることとします。議論を進めるため、その確率を以下のように定義します。

A の出る確率: 1/3
B の出る確率: 1/3
C の出る確率: 1/3

これは、簡単に言えば三面体のサイコロのそれぞれの面がでる確率と同じです。サイコロの方が直感的にイメージしやすいため、問題を3面体のサイコロをを2000回投げた時、それぞれの面のでる回数の分布を求める問題に置き換えます。A に分類されるか、B または C に分類されるかの確率が決まっていますので、この確率を前提として、シミュレーションを行います。

4. シミュレーションを行うコード

シミュレーションを行うコードは以下になります。

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import scipy.stats

# 1. 
samples = np.random.choice(['A', 'B', 'C'], p=[1/3, 1/3, 1/3], size=(10000, 2000))
# 2.
df = pd.DataFrame(samples)
# 3.
df = df.T.apply(lambda x: x.value_counts())
# 4.
plt.hist(df.loc['A', :], bins=50)
plt.show()

# 5.
mean = df.loc['A', :].mean()
sigma = df.loc['A', :].std()
print('平均:{}回\n標準偏差:{}'.format(mean, sigma))

# 6.
n, bins, patches = plt.hist(df.loc['A', :], bins=50, density=True) 
# 7.
y = scipy.stats.norm.pdf(loc=mean, scale=sigma, x=bins)
plt.plot(bins, y)
plt.show()

# 8.
print("'A'のでる回数xは...")
print("68% は {:.0f} < x < {:.0f}となる".format(mean - sigma, mean + sigma))
print("95% は {:.0f} < x < {:.0f}となる".format(mean - 2*sigma, mean + 2*sigma))
print("99.7% は{:.0f} < x < {:.0f}となる".format(mean - 3*sigma, mean + 3*sigma))

5. コードの説明

コードの説明は以下の通りです。

1. まず、 `np.random.choice` 関数で 'A', 'B', 'C' をそれぞれ 1/3 の確率で 2000個 出現する配列を作成します。 第2引数は第1引数で指定したリストの要素が出現する確率を指定しています 。今回はさらに、 それを10000回行ったこととし、 size 引数に (10000, 2000) を渡すことで 10000行 の行列を作成します。
2. 次に 1. で作成した numpy配列を DataFrame に変換します。
3. df を転置させたものの全ての列に対して、 `Series.value_counts` メソッドを適用します。
4. 3. で求めたものをヒストグラムで可視化します。形状が正規分布となっていることがわかります。
5. 平均と標準偏差を求めます。
6. `scipy.stats.norm.pdf` 関数のパラメータに先ほどの平均と標準偏差と 3. で表示したヒストグラムのバーのx座標の値を渡し、その母数での正規分布のx座標に対するy座標の値を得ます。
7. 3. で表示したヒストグラムを面積が1となるように正規化したグラフと、6. で求めた正規分布のグラフを重ねて表示します。
8. 正規分布の性質(平均値から±標準偏差の3倍までの面積は 0.997 、2倍までの面積は 0.95 、1倍までの面積は 0.68) から、範囲に対する確率を求め、表示します。

6. 出力結果

出力結果は以下になります。

平均:666.8196回
標準偏差:21.023074294814155
'A'のでる回数xは...
68% は 646 < x < 688となる
95% は 625 < x < 709となる
99.7% は604 < x < 730となる

f:id:keimina:20191005012451p:plain
ヒストグラム1

f:id:keimina:20191005012533p:plain
ヒストグラム1(正規化済)&正規分布

7. 考察・まとめ

問題を定義し、仮定を定義しました。三面体のサイコロ投げで面が出る回数について10000回シミュレーションを行い、その確率の分布を求めました。サイコロの目が出る回数を数えたら結局は正規分布になりました。今回行ったシミュレーションから 99.7% の割合で 'A' が出現する回数xは 604 < x < 730 の範囲となると推測しました。 この範囲外にあった場合、それは割とよく発生するものでない(異常な状態)ということがわかります。ただし、この推測が役に立つのは、あくまで仮定が成立する時の話で、この仮定とデータが何ら紐づかない場合は実際あまり意味のない議論になります。しかし、データと仮定が紐づかない場合でも、このような条件が議論する上で重要であるということがわかっているのであれば、何の指標も方向性も定まっていない状態より良いのではないかと思います。

以上です。