混合分布が確率であることについて

混合分布が確率であることについて書きます。混合分布とは、確率質量関数または確率密度関数を線形結合させたものだと理解しています。本を読んでいたところ、結合係数の和が1のとき、混合分布も確率質量関数または確率密度関数として解釈することができるということでした。個人的にこれがなぜかわからなかったので備忘録としてまとめました。

Table of Contents

• 1. わからなかったこと
• 2. このように考えたらわかるようになった

1 わからなかったこと

まず、わからなかったことについて書きます。「このように考えたらわかるようになった」ということについては次章に書きます。

話を簡単にするため、ここでは、離散的な事象についての確率質量関数に限定して話を進めます。例えば、四面体のサイコロと、いかさまサイコロの二つの確率質量関数の分布が以下だったとします。※サイコロがいかさまか否かは本題ではありません。

普通のサイコロ:

X	1	2	3	4
P1(X)	0.25	0.25	0.25	0.25

いかさまサイコロ:

X	1	2	3	4
P2(X)	0.1	0.2	0.3	0.4

このように、離散的な確率分布であれば、イメージしやすいです。まずは、二つの確率分布がイメージできました。これら二つの確率質量関数をそれぞれ P1(X), P2(X) とします。

※Xは確率変数

この二つの確率分布を用いて混合分布の関数 P3 を作成すると以下になります。

P3(X) = w1*P1(X) + w2*P2(X)

※w1, w2 はスカラで、結合係数

X	1	2	3	4
P3(X)	w1*0.25 + w2*0.1	w1*0.25 + w2*0.2	w1*0.25 + w2*0.3	w1*0.25 + w2*0.4

このとき

w1 + w2 = 1

であれば

P3(1) + P3(2) + P3(3) + P3(4)

が成り立ち、P3 は確率質量関数となる。

上記のがなぜそう言えるのかがわかりませんでした

2 このように考えたらわかるようになった

P3(1) + P3(2) + P3(3) + P3(4)

を計算すると確かに 1 になるので確率なのですが、次のように式を変形し

X	1	2	3	4
P3(X)	w1*0.25	w1*0.25	w1*0.25	w1*0.25
	+	+	+	+
	w2*0.1	w2*0.2	w2*0.3	w2*0.4

このように分解して見ると理解しやすかったです。

そして、

w1*0.25 + w1*0.25 + w1*0.25 + w1*0.25

は

w1*(0.25 + 0.25 + 0.25 + 0.25)

すなわち

w1*1.0

すなわち

w1

結局のところ、「足して1になるベクトルの要素のそれぞれに w1 をかけて、合計を求める」ことは「1にw1をかける」ことと同じということです。

Python でこれを表現すると以下のようになるでしょうか、、、

import numpy as np
from sympy import symbols, solve, Eq, simplify
w1 = symbols("w1")
a1, a2, a3, a4 = symbols("a1 a2 a3 a4") # [0.25, 0.25, 0.25, 0.25] を [a1, a2, a3, a4] としています。
vec = np.array([a1, a2, a3, a4])
s = sum(w1*vec)
print(s.factor())

# a1 + a2 + a3 + a4 = 1 の条件を s に反映して s を求める
left = a1 + a2 + a3 + a4
a1_new = solve(Eq(left, 1), a1)[0]
# s に 条件を反映
print(simplify(s.subs({a1:a1_new})))

w1*(a1 + a2 + a3 + a4)
w1

同様に

w2*0.1 + w2*0.2 + w2*0.3 + w2*0.4

は

w2(0.1 + 0.2 + 0.3 + 0.4)

すなわち

w2*1.0

すなわち

w2

です。

これらをまとめると

P3(1) + P3(2) + P3(3) + P3(4)

は

w1 + w2

と等しくなり

w1 + w2 = 1

であるならば

P3(1) + P3(2) + P3(3) + P3(4) = w1 + w2 = 1

となるため、P1(X) と P2(X) を結合係数 w1, w2(ただし w1 + w2 = 1 のとき) で線形結合させた関数 P3(X)は確率質量関数としてみなせる。

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

if文を自動生成する方法

if文を自動生成する方法について記載します。

Table of Contents

• 1. したいこと
• 2. 解決策 (1)
• 3. 解決策 (2)
• 4. sympy.SOPform 関数について理解を深める
• 5. 参考文献

1 したいこと

以下のような表があったとき、入力 (a, b) を f(a, b) に変換する関数を自動で作成したいと思います。

a	b	f(a,b)
0	0	1
0	1	0
1	0	1
1	1	0

表1.

この関数を手動でプログラムすると以下のようになるかと思います。

コード:

def f(a, b):

        if a==0 and b==0:

                ret = 1

        if a==0 and b==1:

                ret = 0

        if a==1 and b==0:

                ret = 0

        if a==1 and b==1:

                ret = 1

        return ret

上記の関数の出力が表1.の出力と同じになるか確認します。

コード:

from itertools import product

for x in product([0, 1], [0, 1]):

        line = map(str, x + (f(*x),))

        print("|".join(line))

'''

コメント:

itertools.product はドット積のパターンを生成するのに使用しています。

for文を入れ子にしなくて良くなるので便利です。'''

出力:

0|0|1

0|1|0

1|0|0

1|1|1

したいことをまとめると、、、

上記の表1.のような、入力と出力の関係を定義した表が与えられたとき、上記「def f(…):」の中でプログラムしたようなif文を自動生成生成したいということです。

2 解決策 (1)

解決策 (1) では if文を自動生成するプログラムを自作します。以下のように実装しました。

コード:

import pandas as pd
# 入力となる表を DataFrame として作成する

df = pd.DataFrame(

        [[0, 0, 1],

         [0, 1, 0],

         [1, 0, 0],

         [1, 1, 1]],

        columns=['a', 'b', 'ret']

)
# 評価を入力とし用いて、関数fを定義するための文字列を作成する

arg_t = ",".join(df.columns[:-1])

if_t = ""

for l in df.values:

        s = " and ".join("%s==%d"%(col, i) for col, i in zip(df.columns[:-1], l[:-1]))

        if_t += "    if %s:\n"%s

        if_t += "        ret = %d\n"%l[-1]


f_t = "def f({}):\n"\

          "{}"\

          "    return ret".format(arg_t, if_t)
# 上で定義した文字列をインタプリンタで評価し、関数fを上書きする

exec(f_t)
# 関数fが期待通りの動作になっているか確認する

from itertools import product

for x in product([0, 1], [0, 1]):

        line = map(str, x + (f(*x),))

        print("|".join(line))

出力:

0|0|1

0|1|0

1|0|0

1|1|1

if文をテキストとして書いてそれを exec関数で評価して、関数を定義しています。その関数を使用して表と同じ出力を出力しています。

※exec関数は、文字列をインタプリンタの入力として評価します。予期しないバグの温床になるので使用しないようにしましょう。また、任意のコードを実行可能なので実行するのはセキュリティ上の危険が伴います。ソースコードをよく理解した上で使用しましょう。

念のため生成した関数を表示してみます.

コード:

print('====================作成した関数の表示====================')

print(f_t)

出力:

====================作成した関数の表示====================

def f(a,b):

    if a==0 and b==0:

        ret = 1

    if a==0 and b==1:

        ret = 0

    if a==1 and b==0:

        ret = 0

    if a==1 and b==1:

        ret = 1

    return ret

問題なく自動生成されていることがわかります。

3 解決策 (2)

解決策としては 解決策 (1) でも問題ないのですが、もう少しスマートな方法はないでしょうか。調べたところ sympy.SOPform 関数このような問題を解決するのにうってつけな関数でした。sympy.SOPform 関数を使用すると、表1.のような表を与えると、それに紐づく論理式を組み立ててくれます。さらに、冗長な論理式であった場合は冗長でない論理式にしてくれます。

コード:

import pandas as pd

from sympy import symbols

from sympy.logic import SOPform
# 入力となる表を DataFrame として作成する

df = pd.DataFrame(

        [[0, 0, 1],

         [0, 1, 0],

         [1, 0, 0],

         [1, 1, 1]],

        columns=['a', 'b', 'ret']

)
ss = symbols(" ".join(df.columns[:-1]))

idx = df.iloc[:, -1]==1
if any(idx):

        variables = ss

        minterms = df[idx].values[:,:-1].tolist()

        g = SOPform(variables, minterms)# 出力が 1 になるパターンだけを第二引数に設定する


# 評価を入力とし用いて、関数fを定義するための文字列を作成する

arg_t = ",".join(df.columns[:-1])

f_t = "def f({0}):\n"\

          "    args_bak = {0}\n"\

          "    {0} = symbols('{0}')\n"\

          "    expr = {1}\n"\

          "    return int(bool(expr.subs(dict(zip([{0}], args_bak)))))".format(arg_t, str(g))
exec(f_t)
# 関数fが期待通りの動作になっているか確認する

for line in df.values.tolist():

        line = line[:-1] + [f(*line[:-1])]

        print("|".join(str(i) for i in line))

出力:

0|0|1

0|1|0

1|0|0

1|1|1

作成した関数も見てみます。

コード:

print(f_t)

出力:

def f(a,b):

    args_bak = a,b

    a,b = symbols('a,b')

    expr = (a & b) | (~a & ~b)

    return int(bool(expr.subs(dict(zip([a,b], args_bak)))))

「(a & b) | (~a & ~b)」が sympy.SOPform 関数で作成した論理式です。自動で論理式を作成してくれて大変便利です。

次は、もっと大きい表で試してみます。まず、先ほどより大きめの表を作成します。

コード:

import numpy as np

from itertools import product
np.random.seed(1)

ret = np.random.choice([0, 1], 16)

columns = ['a', 'b', 'c', 'd']

df = pd.DataFrame(list(product([0,1], [0,1], [0,1], [0,1])), columns=columns)

df['ret'] = ret

print(df)

出力:

    a  b  c  d  ret

0   0  0  0  0    1

1   0  0  0  1    1

2   0  0  1  0    0

3   0  0  1  1    0

4   0  1  0  0    1

5   0  1  0  1    1

6   0  1  1  0    1

7   0  1  1  1    1

8   1  0  0  0    1

9   1  0  0  1    0

10  1  0  1  0    0

11  1  0  1  1    1

12  1  1  0  0    0

13  1  1  0  1    1

14  1  1  1  0    1

15  1  1  1  1    0

この表に対して、 前の手順と同じことをします。

コード:

from sympy import symbols

from sympy.logic import SOPform
ss = symbols(" ".join(df.columns[:-1]))

idx = df.iloc[:, -1]==1
if any(idx):

        variables = ss

        minterms = df[idx].values[:,:-1].tolist()

        g = SOPform(variables, minterms)
# 評価を入力とし用いて、関数fを定義するための文字列を作成する

arg_t = ",".join(df.columns[:-1])

f_t = "def f({0}):\n"\

          "    args_bak = {0}\n"\

          "    {0} = symbols('{0}')\n"\

          "    expr = {1}\n"\

          "    return int(bool(expr.subs(dict(zip([{0}], args_bak)))))".format(arg_t, str(g))
exec(f_t)
print('========関数fが期待通りの動作になっているか確認する========')

# 関数fが期待通りの動作になっているか確認する

lst = []

for line in df.values.tolist():

        line = line[:-1] + [f(*line[:-1])]

        print("|".join(str(i) for i in line))

        lst.append(line)

print('===================df と lst を比較する====================')

print('df と lst の', '比較OK' if df.values.tolist() == lst else '比較NG')

print('====================作成した関数の表示====================')

print(f_t)

出力:

========関数fが期待通りの動作になっているか確認する========

0|0|0|0|1

0|0|0|1|1

0|0|1|0|0

0|0|1|1|0

0|1|0|0|1

0|1|0|1|1

0|1|1|0|1

0|1|1|1|1

1|0|0|0|1

1|0|0|1|0

1|0|1|0|0

1|0|1|1|1

1|1|0|0|0

1|1|0|1|1

1|1|1|0|1

1|1|1|1|0

===================df と lst を比較する====================

df と lst の 比較OK

====================作成した関数の表示====================

def f(a,b,c,d):

    args_bak = a,b,c,d

    a,b,c,d = symbols('a,b,c,d')

    expr = (b & ~a) | (~a & ~c) | (b & c & ~d) | (b & d & ~c) | (a & c & d & ~b) | (~b & ~c & ~d)

    return int(bool(expr.subs(dict(zip([a,b,c,d], args_bak)))))

作成した関数の出力と、表で定義した関数の出力が同じであることを確認し、問題なく、動作していることがわかります。

4 sympy.SOPform 関数について理解を深める

sympy.SOPform 関数を使用して自動で求まった論理式ですが、以下の sympy.logic.boolalg.to_dnf 関数を通しても同じなります。

コード:

from sympy.logic.boolalg import to_dnf

print(to_dnf(expr, simplify=True))

出力:

(b & ~a) | (~a & ~c) | (b & c & ~d) | (b & d & ~c) | (a & c & d & ~b) | (~b & ~c & ~d)

to_dnf関数とはなんでしょうか？関数のドキュメントをみてみます。

sympy.logic.boolalg.to_dnf(expr, simplify=False)
Convert a propositional logical sentence s to disjunctive normal

form. That is, of the form ((A & ~B & …) | (B & C & …) | …) If

simplify is True, the expr is evaluated to its simplest DNF form using

the Quine-McCluskey algorithm.

引数 simplify=True のとき Quine-McCluskey アルゴリズムで最も単純なDNF形式に変換されるとあります。参考文献No4. によると Quine–McCluskey アルゴリズムはブール関数を最小化するために使用されるアルゴリズムとのことです。また、参考文献No2. No3. から DNF とは論理式を表現する時に使える演算子を制限した論理式のフォーマットのようなものだと解釈しました。このような、分野をブール関数の簡単化(参考文献No6.)というようです。

まとめると、解決策 (2) では sympy.SOPform 関数を使用して、真理値表からブール関数を自動的に構築するだけでなく、そのブール関数の簡単化を行ったということになります。

5 参考文献

参考文献は以下の通りです。

No.	タイトル	URL	内容
1	Logic	https://docs.sympy.org/latest/modules/logic.html	sympy.SOPform関数について
2	連言標準形	https://ja.wikipedia.org/wiki/連言標準形	CNFについて
3	選言標準形	https://ja.wikipedia.org/wiki/選言標準形	DNFについて
4	Quine–McCluskey algorithm	https://en.wikipedia.org/wiki/Quine–McCluskey_algorithm	Quine–McCluskey algorithmについて
5	クワイン・マクラスキー法	https://ja.wikipedia.org/wiki/クワイン・マクラスキー法	クワイン・マクラスキー法について
6	ブール関数	https://ja.wikipedia.org/wiki/ブール関数	ブール関数について、簡単化について

以上です。

sympyで「... = ...」の式を「... = 0」の式にする方法

sympyで「… = …」の式を「… = 0」の式にする方法を記載します。

Table of Contents

• 1. 行いたいこと
• 2. 解決策
• 3. 練習問題
  • 3.1. 【問 1】
  • 3.2. 【問 2】

1 行いたいこと

まず、行いたいことを明確にします。例えば、以下の式があったとします。

a = b

上記の式の両辺から b を引く(b を左辺に移項する)と以下のようになります。

a - b = 0

今回行いたいことは、上記のように両辺から右辺の式を引く、または左辺の式を引くことにより「… = 0」の式を作成することです。

2 解決策

以下が解決策です。まず、 a = b を表す等式を作成する方法は以下になります。

from sympy import Symbol, Eq
a = Symbol('a')
b = Symbol('b')
left = a
right = b
e = Eq(left, right)
print(e)

Eq(a, b)

上記は、a と b を表す Symbol オブジェクトを left, right 変数に格納しています。a = b を表す等式を変数 e に格納するということを行っています。等式の作成には sympy.Eq 関数を使用しています。

両辺から右辺の式を引き算すれば「… = 0」の式になります。これを、sympy で書くと以下になります。

# 両辺から右辺の式 right を引く
e2 = Eq(left - right, right - right)
# 上記プログラムは以下と同じ
e2 = Eq(left - right, 0)
print(e2)

Eq(a - b, 0)

sympyで「… = …」の式を「… = 0」の式にする方法の説明は以上です。物足りないので練習問題を追加しておきます。

3 練習問題

練習として以下を解くことにします。

3.1 【問 1】

問:

以下の等式を移項して「... = 0」の等式にするにはどうすれば良いでしょうか？

　　a = (b - 1.08614074 * a)**2

答え:

# Symbol の定義
a = Symbol('a')
b = Symbol('b')
# 右辺、左辺の定義
left = a
right = (b - 1.08614074 * a)**2
# (変形前)等式の定義
e = Eq(left, right)
print("変形前: ", e)
# (変形後)等式の定義
e2 = Eq(left - right, 0)
print("変形後: ", e2)

変形前:  Eq(a, (-1.08614074*a + b)**2)
変形後:  Eq(a - (-1.08614074*a + b)**2, 0)

3.2 【問 2】

問:

【問 1】の式を b について解く(b = f(a) の形に変形する)にはどうす
れば良いでしょうか？

※式を解くには sympy.solve 関数を使用します。第一引数に等式を指定します。
【問 1】のソースコードの e, e2 どちらを使用しても問題ありません)第二引
数にどのSymbolについて解くかを指定します。戻り値は要素が解のリストです。

答え:

from sympy import solve
# sympy.solve 関数を使用して b について解く
answer = solve(e, b)
# 上記プログラムは以下と同じ
answer = solve(e2, b)
print(answer)

[-sqrt(a) + 1.08614074*a, sqrt(a) + 1.08614074*a]

以上です。

ヘビサイド関数を使って矩形関数を作る方法(Sympy)

0 or 1 をとるシンプルな関数として、ヘビサイド関数があります。本記事は、 Python のライブラリである Sympy の Heaviside 関数を使って、矩形関数を実装したときの記事です。ヘビサイド関数を用いて矩形関数を作ることができました。

Table of Contents

• 1. 参考文献
• 2. ヘビサイド関数とは
• 3. 矩形関数とは
• 4. 矩形関数をヘビサイド関数で表すには
• 5. コード
• 6. 解説
• 7. 出力
• 8. まとめ

1 参考文献

以下のサイトが参考になりました。

oeis.org

2 ヘビサイド関数とは

ヘビサイド関数とは以下のようなグラフを出力する関数です。

図1. ヘビサイド関数

3 矩形関数とは

矩形関数とは以下のようなグラフを出力する関数です。

図2. 矩形関数

4 矩形関数をヘビサイド関数で表すには

矩形関数をヘビサイド関数で表すには、簡単に言ってしまうと、上に凸で正の値をとる２次関数の値をヘビサイド関数に与えると、矩形関数になります。以下にこれを実装したコードを記載します。

5 コード

矩形関数をヘビサイド関数で実装したコードは以下になります。

# ライブラリ読み込み(数値計算、グラフ出力用)
import numpy as np
import matplotlib.pyplot as plt

# ライブラリ(数式モデル構築用)
from sympy import Heaviside

# ヘビサイド関数は sympy.Heaviside を使用する
# numpy のベクトルに関数を適用できるように np.vectorize を使用する
f = lambda x: Heaviside(-(x - 2)*(x + 2), 0)
f = np.vectorize(f)

# 代入する x の値を定義する
xdata = np.linspace(-6, 6, 100)

# 代入する y の値を求める
ydata = f(xdata)

# 求めた xdata, ydata をプロットする
plt.plot(xdata, ydata)
plt.show(False)

6 解説

ヘビサイド関数は sympy.Heaviside を使用します。 sympy.Heaviside ではx=0 の時の戻り値を第２引数で設定することができます。上記コードでは第２引数に 0 を指定していることに注意します。

7 出力

以下のグラフが出力されます。-(x - 2) * (x + 2) が正になる範囲、すなわち -2 < x < 2 の範囲で 1 に、それ以外の範囲では 0 になることがわかります。

図3. 矩形関数

8 まとめ

sympy ライブラリの関数 Heaviside を用いて、 矩形関数を実装しました。上に凸で正の値をとる２次関数の値をヘビサイド関数に与えると、矩形関数になることがわかりました。

以上です。

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

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

Table of Contents

• 1. 参考文献
• 2. 環境構築
  • 2.1. 前提条件
  • 2.2. "plaidml-keras" のインストール
  • 2.3. 初期設定コマンド
• 3. Keras で GPU を使用するために追加するコード
• 4. MNIST の学習
• 5. まとめ

1 参考文献

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

• Tensorflowのdockerを使ってみる（macOS）

qiita.com

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

PlaidML Kerasでやっていく #TokyoR 73 from Akifumi Eguchi

www.slideshare.net

• Plaid ML Installation Instructions macOS

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 Pro の Radeon GPU を使用することができました。

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

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

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

Table of Contents

• 1. ボロノイ図とは
• 2. ボロノイ図を書くプログラム実装する
  • 2.1. コード
  • 2.2. コードの説明
• 3. 出力結果
• 4. マウスクリックしたときに一番近いデータポイントを表示するプログラムを実装する
  • 4.1. コード
• 5. まとめ

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'])

# データをブロードキャスト(データを３階のテンソルへ変換)
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))

# ２乗距離が最小の平面のインデックスを求める
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 ３次元の r 配列の要素である行列はデータポイントとの距離を表しています。
29	z = r.argmin(axis=0) ３次元の r 配列の0軸方向で一番値の小さいインデックスを求めています。これは、xy座標上のある領域において、その領域では、どのデータポイントが一番近いか求めているということです。
33	ax.plot(df.x, df.y, lw=0, marker='o') plot 関数により、データポイントの描画を行います。 scatter でも良いでしょう。
34	ax.pcolormesh(mx, my, z) pcolormesh 関数により、領域に色を付けます

3 出力結果

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

図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)

# データをブロードキャスト(データを３階のテンソルへ変換)
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))

# ２乗距離が最小の平面のインデックスを求める
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)

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

5 まとめ

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

以上です。

Author: MINAGAWA KEI

Created: 2019-12-15 日 20:24

Validate

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

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

Table of Contents

• 1. 今回、作成したい配列
• 2. 作成する方法 その1
• 3. 作成する方法 その2
• 4. コード
• 5. 出力結果
• 6. コードの説明
• 7. まとめ

1 今回、作成したい配列

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

図1. 作成したい配列

2 作成する方法 その1

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

1. １次元配列の一つの要素を取り出して、その取り出した値を要素とする行列を作成する
2. 作成した行列をリストにappend
3. １次元配列の全ての要素について 上記(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 を配列を要素数をそのままに次元数を３次元にを増やした配列を表しています。np.newaxis は None で定義されていますが、このように使用すると、配列に新しい軸を与え、次元を増やすことができます。第２引数は、変換したい配列の shape を与えています。

7 まとめ

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

以上です。

Author: MINAGAWA KEI

Created: 2019-12-15 日 13:43

Validate

keimina 2019-12-15 14:10

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

Tweet

• もっと読む

コメントを書く