2019-04-27

SciPy Japan 2019 に参加しました

前に申し込んでいた SciPy Japan 2019 に参加してきました。

予定通り4/23, 4/24 に開催されました。無事に開催されて良かったです。開催の内容は、午前中はチュートリアル、午後は Enthoughtさんによる公演や Science に関する問題を(一部) Python を使って解決したというようなトークのセッションでした。私は個人で参加したのですが、企業の一員として来られている方が多い印象でした。あっさりと書いていますが、トークの一部は私が求めていたような内容そのもので(トークの選別した人グッジョブです)感動しました、発表者の方が非常にわかりやすく説明してくれましたし、ワクワク感を味わえました。どれくらいワクワクしたかというと、私が小学生だったらこの発表を聞いたらサイエンスの分野に足を踏み入れたくなるくらいです。

SciPy のコミュニティを推進している人たち(主にEnthoughtさん)の発表を聞いて思ったのは、 SciPy の根底にあるものは Science だということです。Science に主眼を置いて、 Python やそれに付随する各種ライブラリやフレームワークは、あくまで目的達成のための手段に過ぎないということです。それなら SciPy ではなく SciCode でいいじゃないかと思うのですが、コードを書くには Python の方が都合が良かったため SciPy と名付けたということです(と私は解釈した)。実際その決断は間違ってなさそうで、例えば大学の講師がデータ処理にPythonを主要言語に選択したなどのような話は、その決断の正しさに一票を追加するものでしょう。今はPython言語の勢いがありますが、SciPy コミュニティは、この勢いに負けずに、 Science を根底とするというコンセプトは変えないで欲しいと思います。

f:id:keimina:20190424184727j:plain

今の所、SciPy は 2020年も日本に来てくれるようです、個人的に、それまでに以下を行いたいと思っています。

1. 英語力(聞く、話す)を身に付ける
2. SciPy っぽいコミュニティを作る or 参加する

以上

2019-03-25

sympy でニューラルネットワークの重みの更新に使用する式を計算グラフにしてみる

python sympy 計算グラフ自動微分 neural network deep learning 勾配降下法最急降下法

1. はじめに

sympy(https://www.sympy.org/en/index.html) を使用し、ニューラルネットワークの重みの更新に使用する式を計算グラフにしてみます。具体的には、 sympy を使って、損失関数の勾配の計算を計算グラフとして表し、それを graphviz で画像に出力します。sympy と graphviz を使用すると、数式を簡単に計算グラフに変換できます。ちなみに計算グラフとは、以下のようなものです。(実際にこれらを使用し作成したものです)

f:id:keimina:20190325203008p:plain — 図1. c = a + b の計算グラフ

このように計算式をグラフとして表現したものを計算グラフと呼びます。

2. 使用するモジュールと、そのインストール

本記事では、以下の二つのモジュールを使用します。グラフの描画に graphviz というモジュールを使用します。sympy は計算グラフをテキストで出力できますが、単体では画像で出力することができないためグラフ描画アプリ(graphviz)を使用します。

graphviz は以下のコマンドでインストールしました。

pip install graphviz

コマンドの出力結果

Collecting graphviz
  Downloading https://files.pythonhosted.org/packages/1f/e2/ef2581b5b86625657afd32030f90cf2717456c1d2b711ba074bf007c0f1a/graphviz-0.10.1-py2.py3-none-any.whl
distributed 1.21.8 requires msgpack, which is not installed.
Installing collected packages: graphviz
Successfully installed graphviz-0.10.1
You are using pip version 10.0.1, however version 19.0.3 is available.
You should consider upgrading via the 'pip install --upgrade pip' command.

コマンド実行後、"requires msgpack" と出力されていますが、最終的に"Successfully installed" と出力されているため、問題なしとみなし気にせず先に進みました。

今回は、計算グラフを画像で出力したいため、さらに以下の手順行う必要がありました。

MacPorts のインストール(インストール手順は以下のサイトを参照)
- https://www.macports.org
以下のコマンドの実行しアプリをインストール

sudo port install graphviz

sympy のインストール手順について説明は省略します。

3. ニューラルネットワークのモデルと、順伝播の式

本記事で扱うニューラルネットワークのモデルと、順伝播の式は以下の通りです。

f:id:keimina:20190325203628p:plain — 図2. ニューラルネットワークのモデルと、順伝播の式

4. ソースコード

3. のニューラルネットワークの損失関数の勾配を計算するコードにすると以下のようになります。
※本記事は、重みの更新を行う計算を、計算グラフに出力することを目的としているため、説明をわかりやすくするため損失関数 l は o1 としています。本来であれば、l = (o1 - y)**2 のような定義をすべきですが、パラメータが増えると、計算グラフが大きくなり見づらくなるため、簡略化しています。ご了承ください。
※ソースコードで numpy を使用して行列計算していますが、sympy での行列計算の仕方がわからなかったためです。あまり真似しないように。。。

from sympy import symbols, Function, tanh
from sympy.printing import dot
import graphviz
import numpy as np

# sympy で使用する変数(シンボル)を定義する
# 入力値、中間層の値
i1, i2 = symbols("i1, i2")
h11, h12 = symbols("h11, h12")

# 重み
w111, w112, w121, w122 = symbols("w111, w112, w121, w122")
w211, w212 = symbols("w211, w212")

# バイアス
b11, b12 = symbols("b11, b12")
b21 = symbols("b21")

# numpy を利用して行列形式にする
I = np.array([[i1], [i2]])
W1 = np.array([[w111, w112], [w121, w122]])
W2 = np.array([[w211, w212]])
B1 = np.array([[b11], [b12]])
B2 = np.array([[b21]])

# 活性化関数
TANH = np.vectorize(tanh)

# 順伝播の数式を作成する
H1 = TANH(W1.dot(I) + B1)
O = TANH(W2.dot(H1) + B2)

# 行列からスカラを取り出す
o1 = O[0, 0]

# 最急降下法 によって重み w1 を更新する式(学習係数は 1 とする)
# 損失関数を定義する、重み更新を行う計算グラフを出したいだけなので l = o1  としてしまう
l = o1
dl_dw111 = l.diff(w111)

expr_forward = o1
src_forward = graphviz.Source(dot.dotprint(expr_forward))

expr_backward = dl_dw111
src_backward = graphviz.Source(dot.dotprint(expr_backward))

# 私の環境ではExecutableNotFound エラーが発生したので、以下のコードで回避しました
app_path = '/opt/local/bin/dot'
graphviz.ENGINES.add(app_path)
src_forward._engine = app_path
src_backward._engine = app_path

# 画像ファイルを作成する
filename_forward = "calc_forward_graph"
src_forward.render(filename_forward)
filename_backward = "calc_backward_graph"
src_backward.render(filename_backward)

順伝播の式が o1 です。勾配降下法によって損失関数 l を最小にするため重みパラメータ w111 を更新するときに行われる計算が dl_dw111 です。コードの最後の方で、sympy と graphviz を使用して、o1 と dl_dw111 の計算を行う計算グラフを出力しています。sympy.printing.dot で計算グラフを graphviz で使用可能な dot フォーマット形式のテキストに変換しています。そして、 dot フォーマット形式のテキストを graphviz.Source で graphviz に取り込み、 render メソッドを呼びだし pdf ファイルを作成しています。

5. 出力結果

出力結果は以下になります。図は上から順に、順伝播の計算を行うための計算グラフ(o1 の計算を行うグラフ)、損失関数lを最小化すべく w111 を勾配降下法により更新するための計算グラフ(dl_dw111 の計算を行うグラフ)、となっています

f:id:keimina:20190325204132p:plain — 図3. 順伝播の計算グラフ("calc_forward_graph")

f:id:keimina:20190325204331p:plain — 図4. w111 の重みを更新する数式の計算グラフ("calc_backward_graph")

楕円の中に書かれている Mul 、 Add 、 Pow 、 tanh はどのような計算を行うかを表していて、それぞれ、乗算、加算、累乗、ハイパボリックタンジェントを表しています。ニューラルネットの学習時には、図中の w111、w112 、w121、 b11 、、、などの変数のように見えるものは全て値が割り振られているため、この計算グラフの計算を行うと、何かしらの値(スカラ値)が求まります。そしてこの値を使用して重み w111 を更新します。

6. まとめ

sympy と graphviz を使用して、ニューラルネットワークの重みの更新に使用する式を計算グラフにしました。このように、ニューラルネットワークのモデルが決まると、重みを更新するための式、すなわち計算グラフも決まります。

7. 参考文献

ゼロから作るDeep Learning――Pythonで学ぶディープラーニングの理論と実装
- https://www.oreilly.co.jp/books/9784873117584/

2019-03-23

Numpy 入門 (機械学習勉強会 in 新潟 2019/03/23)

Numpy 入門

Numpy 入門
(機械学習勉強会 in 新潟 2019/03/23)

1. はじめに
2. numpy とは何か
3. numpy モジュールのインポート方法について
4. 数値を作成する
- 4.1. np.int32 とは
5. 「１次元配列」の数値計算をする
6. 「２次元配列」の数値計算をする
7. 配列を作成する関数
8. 配列の要素へのアクセス
- 8.1. 配列の要素の取得
- 8.2. 配列の要素の書き換え
9. 真偽値(True/False)について
10. 配列の制約
11. 配列の長さや、形状が異なる場合に行われる計算について
- 11.1. エラーになるか計算結果が求まるかの判定方法
- 11.2. Broadcasting rule の適用例
12. 速度比較
13. まとめ
14. 最後に

1 はじめに

本記事は、機械学習勉強会 in 新潟 2019/03/23 の発表資料です。 numpy の入門者向けの記事になります。 numpy とは何か、何ができるかなどを説明します。numpy で簡単な数値計算を行えるようになることを目的とします。

※時間の都合上、 numpy の機能のほんの一部を紹介しています。本記事は以下の参考文献を参考に書いています。 numpy についての詳細は以下の参考文献を参照してください。

参考文献
1. NumPy v1.16 Manual
  - https://docs.scipy.org/doc/numpy/
2. 現場で使える! NumPyデータ処理入門機械学習・データサイエンスで役立つ高速処理手法
  - https://www.shoeisha.co.jp/book/detail/9784798155913

2 numpy とは何か

Python の数値計算のためのライブラリです。 numpy によってできることとメリットは以下になります

No.	ライブラリ	説明	できること、メリット
1	numpy	数値計算ライブラリ	1. 行列の計算を行うコードを短く書ける
			2. 行列の計算、ベクトルの計算の高速化( 速度比較参照)

3 numpy モジュールのインポート方法について

numpy を使用する場合、numpy モジュールのインポートは以下のように行います。

コード1. 「モジュールのインポート」:

import numpy as np

これにより、'np' の２文字で numpy を使用することができます。

コード2. 「モジュールの使用方法」:

import numpy as np
arr = np.array([1,2,3])
arr

Out: array([1, 2, 3])

使用する名前は numpy など np 以外の名前にもできますが、 np と書くのが一般的になっているためこちらを使用しましょう。

4 数値を作成する

はじめに、numpy で数値を作成しましょう。数値は、スカラ(Scalar)と呼ばれることがあります。これらは、後述する配列の要素( 5.1.1 参照)となりますので、しっかり覚えましょう。 np.int32 を使用すると numpy の整数を作成できます。

コード3. 「numpy で数値の作成」:

import numpy as np
answer = np.int32(100)
answer

Out: 100

4.1 np.int32 とは

np.int32 はデータタイプ(dtype)として定義されています。 np.int32(100) を実行すると np.int32 データタイプのオブジェクトが作成されます。 np.int32 以外にも数値を表わすデータタイプが用意されています。

データタイプ	数値の種類	用途
np.uint8	8bit 符号無し整数	グレースケール画像
np.int32	32bit 整数	画像処理
np.float32	32bit 浮動少数	画像処理、平均値、確率
np.float64	64bit 浮動少数	数値計算(精度が求められるもの)
…

真偽値(True/False)を現す場合は np.bool データタイプが用意されています。

5 「１次元配列」の数値計算をする

配列の基本となる１次元配について学習しましょう。

5.1 「１次元配列」を作成する

１次元の配列を作成しましょう。１次元配列とは以下のようなものです。配列は'要素'と呼ばれるもので構成されています。

値1

値2

値3

,,,

5.1.1 配列の'要素'とは

一般的に、上記 5.1 のような配列の中にある値1や値2などのことを配列の”要素”(Element)と呼びます。１次元配列だけでなく、２次元配列や３次元配列も同様です。配列は要素によって構成されています。

※ numpy では配列を作成する関数の引数に dtype を指定することで、その配列の要素のデータタイプを指定できます。引数 dtype を指定しない場合は、自動的に配列のデータタイプが決定されます。 numpy では配列は同一のデータタイプの要素から構成されます。

5.1.2 np.array 関数

np.array 関数を使用すると、標準 Python のリストを配列に変換することができます.

コード4. 「標準 Python のリストを配列に変換する」:

import numpy as np
arr = np.array([1, 2, 3])
arr

Out: array([1, 2, 3])

5.1.3 np.arange 関数

np.arange 関数を使用すると、数値の１次元配列を作成することができます。

コード5. 「np.arange 関数を使用し、数値の１次元配列を作成」:

import numpy as np
arr = np.arange(3)
arr

Out: array([0, 1, 2])

5.2 「１次元配列」の形状 shape について

配列には shape と呼ばれる、配列の形状を表す情報があります。配列の形状 shape を調べるには以下のようにします。この例の場合、配列は長さ３の１次元配列であることがわかります。

コード6. 「配列の形状 shape を調べる」:

import numpy as np
arr = np.arange(3)
answer = arr.shape
answer

Out: (3,)

5.3 「１次元配列」 + 「１次元配列」

１次元配列と１次元配列の四則演算の基本を学習します。ここでは例として、１次元配列と１次元配列の足し算を行います。以下のように、１次元配列を二つ作成し、足し算記号(+) を使用して求めます。計算は、要素ごと(element wise)に行われます。

コード7. 「１次元配列と１次元配列の足し算」:

import numpy as np
arr_1 = np.array([1, 2, 3])
arr_2 = np.array([0, 0, 1])
answer = arr_1 + arr_2
answer

Out: array([1, 2, 4])

その他の演算子(-*/など)についても、同じように、配列の要素ごとに計算されます。形状の異なる配列同士の計算は、エラーになるか、ならない場合があります。形状の異なる配列同士の計算については 11 を参照してください。

5.4 配列の要素に数学関数を適用する

１次元配列の要素ごと(element wise)に数学関数を適用してみましょう。ここでは例として平方根を適用してみましょう。

コード8. 「配列の各要素に平方根を適用」:

import numpy as np
arr = np.array([1, 2, 3])
answer = np.sqrt(arr)
answer

Out: array([1.        , 1.41421356, 1.73205081])

np.sqrt 関数以外にも様々な数学関数が用意されています。それらの関数も、上記のように使用することで、要素ごとに適用することができます。以下に数学関数の一部を記載します。

No.	関数	説明
1	np.exp	指数関数
2	np.log	対数関数
3	np.sin	sin関数
4	np.arcsin	arcsin関数

5.5 １次元配列の要素の合計を求める

合計を求める場合は、 np.sum 関数を使用します。

コード9. 「np.sum 関数を使用して１次元配列の要素の合計を求める」:

import numpy as np
arr = np.arange(11)
answer = np.sum(arr)
answer

Out: 55

6 「２次元配列」の数値計算をする

２次元配列について学習しましょう。

6.1 「２次元配列」を作成する

２次元の配列を作成しましょう。２次元配列とは以下のようなものです。

値11	値12	値13	,,,
値21	値22	値23	,,,
値31	値32	値33	,,,
,,,

6.1.1 np.array 関数

np.array 関数を使用すると、２次元配列を作成することができます。(１次元配列の時と同様)

コード10. 「np.array 関数を使用して２次元配列を作成する」:

import numpy as np
arr = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
arr

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

6.1.2 reshape メソッド

配列の reshape メソッドを使用すると配列の形状を変えることができます。ここでは、 1次元配列の形状を変えて、２次元配列を作成してみましょう。

コード11. 「reshape メソッドを使用して1次元配列の形状を変えて、２次元配列を作成する」:

import numpy as np
arr = np.arange(9)
answer = arr.reshape((3, 3))
answer

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

6.2 「２次元配列」 + 「２次元配列」

２次元配列同士の四則演算、関数を適用した時の動作は、１次元配列の時と同様、要素ごとに計算が行われます。 (そのため２次元配列の四則演算については説明を省略します。)２次元配列を使用すると、行列の掛け算を行うことができます。

行列の掛け算とは、具体的には以下のような計算のことです。行列 A と行列 B の掛け算(積)は A・B と書き、A と B と A・B は以下のような関係になります。 ※説明の簡略化のため行列の大きさは２x２としています。

行列 A

a b

c d
行列 B

e f

g h
行列 A・B

a*e+b*g a*f+b*h

c*e+d*g c*f+d*h

行列の掛け算を行うには np.dot 関数を使用します。

コード12. 「np.dot 関数を使用して行列の掛け算を行う」:

import numpy as np
arr_1 = np.arange(9).reshape((3, 3))
arr_2 = np.arange(9).reshape((3, 3))
answer = np.dot(arr_1, arr_2)
answer

Out: 
array([[ 15,  18,  21],
       [ 42,  54,  66],
       [ 69,  90, 111]])

配列オブジェクトの dot メソッドを使用して以下のように書くこともできます。

コード13. 「dot メソッドを使用して行列の掛け算を行う」:

import numpy as np
arr_1 = np.arange(9).reshape((3, 3))
arr_2 = np.arange(9).reshape((3, 3))
answer = arr_1.dot(arr_2)
answer

Out: 
array([[ 15,  18,  21],
       [ 42,  54,  66],
       [ 69,  90, 111]])

注意: np.dot は '∗'による掛け算と同じではありません。 '∗' による掛け算は要素ごと(element wise)に掛け算が行わます。一方 np.dot は行列の掛け算ですので、積和演算(掛け算と足し算)が行われます。

6.3 「２次元配列」の要素の合計を求める

'すべての要素の合計'を求める場合は、１次元配列の時と変わりませんので省略しますが、２次元配列では、さらに、'列の合計'と'行の合計'を求めることができます。

6.3.1 列の合計

np.sum 関数の引数で axis=0 とすることで列の合計を求めることができます。

コード14. 「列の合計を求める」:

import numpy as np
arr = np.arange(9).reshape((3, 3))
answer = np.sum(arr, axis=0)
answer

Out: array([ 9, 12, 15])

6.3.2 行の合計

np.sum 関数の引数で axis=1 とすることで行の合計を求めることができます。

コード15. 「行の合計を求める」:

import numpy as np
arr = np.arange(9).reshape((3, 3))
answer = np.sum(arr, axis=1)
answer

Out: array([ 3, 12, 21])

6.3.3 配列の次元数を保つには

6.3.1 、 6.3.2 で np.sum を使用して求めた配列を見ると、２次元配列の入力に対して出力が１次元配列となっており、次元数が変わってしまうことに気付くかと思います。 np.sum 関数の引数で keepdims=True とすることで、配列の次元数を保つことができます。

コード16. 「np.sum 関数の引数で keepdims=True とし配列の構造をなるべく保つ」:

import numpy as np
arr = np.arange(9).reshape((3, 3))
answer = np.sum(arr, axis=1, keepdims=True)
answer

Out: 
array([[ 3],
       [12],
       [21]])

6.4 行列計算用の関数

２次元配列には行列計算用の関数を使用することができます。行列計算用の関数は np.linalg にあります。その一部を紹介します。

関数	説明
np.linalg.inv	逆行列を求めるための関数
np.linalg.norm	行列のノルムを求めるための関数

7 配列を作成する関数

配列を作成する方法は、5.1 と 6.1 で紹介した方法が全てではありません。以下に、配列を作成するための関数の一部を紹介します。

7.1 固定された値を要素にもつ配列を作成する関数

No.	関数	説明
1	np.zeros	全ての要素が 0 の配列を作成します
2	np.ones	全ての要素が 1 の配列を作成します
3	np.full	全ての要素が x の配列を作成します(x は引数で指定)

コード17. 「np.zeros 関数を使用し全ての要素が 0 の配列を作成」:

import numpy as np
arr = np.zeros((3, 3))
arr

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

7.2 単位行列を作成する関数

No.	関数	説明
1	np.eye	単位行列を表す配列を作成します

コード18. 「np.eye 関数を使用し単位行列を表す配列を作成」:

import numpy as np
arr = np.eye(3)
arr

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

7.3 ランダムな値を要素にもつ配列を作成する関数

No.	関数	説明
1	np.random.randint	ランダムな整数の配列を作成します
2	np.random.random	ランダムな少数の配列を作成します

コード19. 「np.random.randint 関数を使用しランダムな整数の配列を作成」:

import numpy as np
arr = np.random.randint(0, 10, (3, 3))
arr

Out:
array([[4, 0, 7],
       [2, 3, 4],
       [3, 0, 5]])

コード20. 「np.random.random 関数を使用しランダムな少数の配列を作成」:

import numpy as np
arr = np.random.random((3, 3))
arr

Out:
array([[0.23963676, 0.14041946, 0.67624721],
       [0.75209081, 0.3182057 , 0.18579186],
       [0.0888395 , 0.87384273, 0.24740674]])

7.4 似たような配列を作成する関数

配列の計算をプログラミングしているときに、ある配列と同じ形状をもつ配列だが、値は違うという配列を作成したい時があります。そのような場合は、以下の関数を使用します。

No.	関数	説明
1	np.zeros_like	引数で指定した配列と同じ形状
		＆
		全ての要素が 0 の配列を作成します

2	np.ones_like	引数で指定した配列と同じ形状
		＆
		全ての要素が 1 の配列を作成します

コード21. 「np.zeros_like 関数を使用し形状が同じ配列を作成」:

import numpy as np
arr = np.random.randint(0, 10, (3, 3))
arr_2 = np.zeros_like(arr)
arr_2

Out[4]: 
array([[0, 0, 0],
       [0, 0, 0],
       [0, 0, 0]])

7.5 ある範囲の数値を要素にもつ配列を作成する関数

No.	関数	説明
1	np.linspace	数直線の範囲を等幅で分割した１次元配列を作成します
		１変数関数を調べる時などに使用できます
2	np.meshgrid	２次元座標などの座標を配列を使い表現するときに使用します。
		２変数関数を調べる時などに使用できます

コード22. 「np.linspace 関数の使用例」:

import numpy as np
arr = np.linspace(0.0, 0.9, 10)
arr

Out: array([0. , 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9])

コード23. 「np.meshgrid 関数の使用例」:

import numpy as np
x = np.linspace(0, 0.4, 5)
y = np.linspace(0, 0.4, 5)
x_mesh, y_mesh = np.meshgrid(x, y)
answer = x_mesh, y_mesh
answer

Out:
(array([[0. , 0.1, 0.2, 0.3, 0.4],
        [0. , 0.1, 0.2, 0.3, 0.4],
        [0. , 0.1, 0.2, 0.3, 0.4],
        [0. , 0.1, 0.2, 0.3, 0.4],
        [0. , 0.1, 0.2, 0.3, 0.4]]), array([[0. , 0. , 0. , 0. , 0. ],
        [0.1, 0.1, 0.1, 0.1, 0.1],
        [0.2, 0.2, 0.2, 0.2, 0.2],
        [0.3, 0.3, 0.3, 0.3, 0.3],
        [0.4, 0.4, 0.4, 0.4, 0.4]]))

8 配列の要素へのアクセス

8.1 配列の要素の取得

添字(インデックス)を使用することで、配列の要素を取得することができます。

8.1.1 要素の取得の基本

以下のようにすると、配列の要素の値を取得することができます。

コード24. 「１次元配列の配列の要素の取得の基本」:

import numpy as np
arr = np.array([1, 3, 5])
arr[2]

Out: 5

コード25. 「２次元配列の配列の要素の取得の基本(1)」:

import numpy as np
arr = np.array([[1, 3, 5], [7, 11, 13], [17, 21, 23]])
arr[2][2]

Out: 23

上記コードは以下のようにも書くことができます。

コード26. 「２次元配列の配列の要素の取得の基本(2)」:

import numpy as np
arr = np.array([[1, 3, 5], [7, 11, 13], [17, 21, 23]])
arr[2, 2]

Out: 23

8.1.2 要素を範囲で指定して取得

以下のようにすると、指定した範囲にある配列の要素を取得することができます。

コード27. 「１次元配列の要素を範囲で指定して取得」:

import numpy as np
arr = np.array([1, 3, 5])
arr[1:]

Out: array([3, 5])

コード28. 「２次元配列の要素を範囲で指定して取得」:

import numpy as np
arr = np.array([[1, 3, 5], [7, 11, 13], [17, 21, 23]])
arr[1:, 1:]

Out: 
array([[11, 13],
       [21, 23]])

8.1.3 逆順の配列の取得

以下のようにすると、逆順の配列を取得できます。

コード29. 「逆順の配列の取得」:

import numpy as np
arr = np.arange(10)
arr[::-1]

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

※8.1.2 や 8.1.3 で取得した arr[1:, 1:] や arr[::-1] は、あたかも配列のコピーのように思われるかもしれませんが、実際は配列のコピーではありません。これらはビューと呼ばれるものであり、実態は arr のデータを参照しているものです。扱いには注意しましょう。 8.2 でも説明します。

8.2 配列の要素の書き換え

8.2.1 要素の書き換えの基本

配列の要素を添字(インデックス)で指定して別の値に書き換えることができます。以下のようにすると、配列の要素の値を書き換えることができます。

コード30. 「１次元配列の配列の要素の書き換えの基本」:

import numpy as np
arr = np.array([1, 3, 5])
arr[2] = 100
arr

Out: array([  1,   3, 100])

コード31. 「２次元配列の配列の要素の書き換えの基本(1)」:

import numpy as np
arr = np.array([[1, 3, 5], [7, 11, 13], [17, 21, 23]])
arr[2][2] = 100
arr

Out:
array([[  1,   3,   5],
    [  7,  11,  13],
    [ 17,  21, 100]])

上記コードは以下のようにも書くことができます。

コード32. 「２次元配列の配列の要素の書き換えの基本(2)」:

import numpy as np
arr = np.array([[1, 3, 5], [7, 11, 13], [17, 21, 23]])
arr[2, 2] = 100
arr

Out:
array([[  1,   3,   5],
       [  7,  11,  13],
       [ 17,  21, 100]])

8.2.2 要素を範囲で指定して書き換え

以下のようにすると、指定した範囲にある配列の要素を書き換えることができます。

コード33. 「１次元配列の要素を範囲で指定して書き換え」:

import numpy as np
arr = np.array([1, 3, 5])
arr[1:] = 100
arr

Out: array([  1, 100, 100])

コード34. 「２次元配列の要素を範囲で指定して書き換え」:

import numpy as np
arr = np.array([[1, 3, 5], [7, 11, 13], [17, 21, 23]])
arr[1:, 1:] = 0
arr

Out: 
array([[ 1,  3,  5],
       [ 7,  0,  0],
       [17,  0,  0]])

上記コードの arr[1:] や arr[1:, 1:] は配列のコピーではないことに注意しましょう。このようなものを numpy ではビュー(View)と呼びます。実態は arr のデータを参照しているものです。このようなビュー arr[1:] や arr[1:, 1:] に対する値の代入は、ビューの参照先である arr に対する値の代入であることを意味します。

よくあるのが、値を保持しておきたい配列のビューに代入を行なってしまいその配列を破壊してしまうことです。値を代入したい場合は、以下のように、配列のコピーを作成して、そのコピーに対して代入しましょう。

コード35. 「配列のコピーを作成して、そのコピーに対して代入する」:

import numpy as np
arr_orig = np.array([[1, 3, 5], [7, 11, 13], [17, 21, 23]])
arr = arr_orig.copy()
arr[1:, 1:] = 0
arr_orig

Out: 
array([[ 1,  3,  5],
       [ 7, 11, 13],
       [17, 21, 23]])

9 真偽値(True/False)について

'==' や '<' などの比較演算子を使用すると配列の要素ごとに比較が行われます。

コード36. 「比較演算子による配列の比較」:

import numpy as np
arr_1 = np.arange(9).reshape((3, 3))**2
answer = arr_1 >= 10
answer

Out: 
array([[False, False, False],
       [False,  True,  True],
       [ True,  True,  True]])

この真偽値を配列のインデックスに指定すると、真偽値のTrue の箇所の要素だけ集めた１次元配列を求めることができます。

コード37. 「真偽値の配列による要素の収集」:

import numpy as np
arr_1 = np.arange(9).reshape((3, 3))**2
answer = arr_1[arr_1 >= 10]
answer

Out: array([16, 25, 36, 49, 64])

以下のように、 np.where を使用すると値が True のインデックスを求めることができます。 np.where の戻り値は、「インデックスを要素にもつ配列」を要素にもつタプル(tuple)です。

コード38. 「np.where による条件を満たす配列のインデックスの取得」:

import numpy as np
arr_1 = np.arange(9).reshape((3, 3))**2
answer = np.where(arr_1 >= 10)
answer

Out: (array([1, 1, 2, 2, 2]), array([1, 2, 0, 1, 2]))

10 配列の制約

ここまで１次元配列、２次元配列の計算方法について説明しました。 numpy の配列は、標準Python のリストと違い、一度作成した配列に要素を追加したり、削除したりすることはできません。一度作成した配列のデータタイプ(dtype)を変更することもできません。これらの操作を行いたい場合は、新たに配列を作成してください。これらの制約を設けることにより高速(C言語並)な計算が可能になっています。

11 配列の長さや、形状が異なる場合に行われる計算について

形状の異なる配列同士の演算を行う場合、それ専用の規則(broadcasting rule)に従い処理が行われます。その規則に従いエラーになることもあれば計算結果が求まることもあります。

11.1 エラーになるか計算結果が求まるかの判定方法

配列 arr₁ と arr₂ のそれぞれの配列の shape を一番右から左にかけて以下条件を満たす時、計算結果が求まります。条件を満たさない場合はエラーとなります。また、エラーにならず、形状が異なる場合は、お互いに異なる形状を補うよう値が補完され、計算されます(詳細は割愛します)。

(以下は「Broadcasting」(https://docs.scipy.org/doc/numpy/user/basics.broadcasting.html) より引用)

両方の値が同じ、または、
片方の数値が1

11.2 Broadcasting rule の適用例

OK

配列 shape

arr₁ 10

arr₂ 1

結果 10
OK

配列 shape

arr₁ 10 1

arr₂ 1 10

結果 10 10
OK

配列 shape

arr₁ 10 10 1

arr₂ 1 1 3

結果 10 10 3
NG

配列 shape

arr₁ 10 10 3

arr₂ 1 1 2

結果エラー

配列	shape
arr₁	10
arr₂	1
結果	10

配列	shape
arr₁	10	1
arr₂	1	10
結果	10	10

配列	shape
arr₁	10	10	1
arr₂	1	1	3
結果	10	10	3

配列	shape
arr₁	10	10	3
arr₂	1	1	2
結果	エラー

12 速度比較

numpy を使用すると、どれくらい速くなるのでしょうか。 100x100の要素からなる２次元配列の各要素を +1 する処理を numpy を使用する前と後で処理時間を比較してみました。処理時間の計測には、マジックコマンド %timeit を使用しています。その結果、numpy の使用により処理速度は 100倍になりました。また可読性も向上したと思います。

コード39. 「速度比較 numpy 使用前」:

def func(x):
   y = x.copy()
   for i in range(100):
       for j in range(100):
           y[i][j] += 1
   return y

# リストの初期化
x = [[0 for i in range(100)] for j in range(100)]
# 処理時間計測
%timeit func(x)

785 µs ± 24.2 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

コード40. 「速度比較 numpy 使用後」:

import numpy as np

def func(x):
  y = x.copy()
  y += 1
  return y

# 配列の初期化
x = np.zeros((100, 100))
# 処理時間計測
%timeit func(x)

7.11 µs ± 284 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

13 まとめ

配列の作成方法、配列同士の計算や、配列への関数の適用方法について学習しました。

numpy を使って得られる恩恵は以下の２つです。

行列計算のコード量が少なくなる
処理速度が早くなる

14 最後に

既存の行列の計算のソースコードの可読性の向上や、処理速度向上などを目的に、 numpy を使用したいと思われるかもしれません。まず計算の対象が 10 の制約の範囲内で処理が可能であることを確認しましょう。次に、以下の表のように可読性と、処理速度を天秤にかけ、 numpy で実装するか標準 Python で実装するか考えると良いのではないかと思います。 (numpy を使ったからといって常に「可読性が高い＆処理速度が早い」コードになる訳ではありません。)

可読性	処理速度	numpy を使用する？
高くなる	高くなる	YES
高くなる	低くなる	YES or NO
低くなる	高くなる	YES or NO
低くなる	低くなる	NO

numpy を使用してコーディングするとテクニカルで難解なコードになってしまうことがあります。可読性重視であればテクニカルなコーディングはしないほうが良いでしょう。 numpy を使用してコーディングする際は、numpy を使用する当初の目的を見失わないようにしましょう。

それでは、よい numpy ライフを!

著者: Kei Minagawa (http://keimina.hatenablog.jp)

Created: 2019-03-23 Sat 11:25

Validate

2019-03-14

SciPy Japan 2019 に参加します

SciPy Japan 2019 に参加します。

www.scipyjapan2019.scipy.org

開催は4/23, 4/24の二日間、場所は東京。チュートリアル、トーク、ポスター発表が行われるようです。チュートリアルも二日間に渡って行われるようです。

英語のリスニングができずに、ほとんど理解できないかもしれず、無謀なのですが、せっかくの機会ですので、参加してみます。果たしてどうなることか、、、

2019-02-25

Wikipediaの文書を使って文の長さの平均を求める

python 統計ビッグデータパラメタ推定 scipy 対数正規分布データサイエンス分析

Wikipediaの文書を使って文の長さの平均値と最頻値を求めました。機械学習勉強会 in 新潟で発表した内容です。Wikipedia の約60億文字以上のXMLファイルから、独自に定義した文を抜き出し、その文に対して統計量を求めました。使用言語はいつも通りPythonです。

1. 動作環境について

コードは以下の環境で動作確認を行っています。

CPU:6コア
メモリ:256GB
HDD:50GB

※環境はGCE(Google Compute Engine)で行っています。実行する場合はメモリ量に注意してください。

2. Wikipedia 文書のアーカイブについて

Wikipedia の文書は以下からダウンロードできます。 "jawiki-latest-pages-articles.xml.bz2"の中にあるXMLが今回解析に使用したファイルになります。

dumps.wikimedia.org

3. 発表スライド

発表スライドは、一番下に置いておきましたので、そちらをご覧ください。

4. やってみた感想

まず、このような分析は、自然言語処理の分野では頻度分析と呼ばれるらしいです。統計ですね。最近は統計のビジネス書をコツコツ読んで勉強していたたため、おかげでサンプルデータの分布を見た瞬間パラメタ推定できると気づくことができました。
パラメタ推定は本で読んでいたけど実際やってみると、ニューラルネットワークなど必要とせず非常にシンプルなのですが、強力でした。今回は、分布が対数正規分布にフィットしたのですがその理由を考察するときに、仮定した関数の性質が役に立ちました。結果、個人的に納得感のある考察を得ることができました（対数正規分布に従うのはなぜか？の考察です）。
最後ですが、Wikipediaを日常的に使用していたため、その中に法則性が存在するという事実に驚きました。おそらく解析しないと一生気付かなかったと思います。個人的に、その法則性の存在を客観的に示せたので非常に有意義だったと思います。また、求めた最頻値の値自体にも、価値があると思います。この値はWikipediaで遭遇する可能性が一番高い文の文字列の長さなので、それは、もしかすると我々が最も親しんでいる文字列長かもしれません。そのため24文字を目指して文を書くとわかりやすさの改善に繋がるつかもしれません。ビッグデータの時代ですが、手付かずのデータが山ほどあるとのことなので、まだ誰も知らない事実がデータの中に大量に眠っていると推測されます。今後それらのデータの中から、想像もしなかった事実が出てくることを期待しています。

5. 発表会で上がった質問や意見

発表会で上がった質問をここにメモしておきます。これらの質問や意見については、今後扱うかもしれませんし、扱わないかもしれません。個人的には、最頻値を求めて十分満足してしまったため、あとやったら面白そうなことは、コーパスで教師なし学習（クラスタリングなど）をさせてテキストマイニングするくらいしか思いつかないです。あと人工無脳。

1. 他の分布で推定してみたらどうなるのか？例えばポアソン分布でやってみたらどうなるか？
→ ググったら一番最初に出てきたキーワードが対数正規分布だったのでそれでやってみた、ポアソン分布ではパラメタ推定していない。（ただし、対数正規分布でパラメタ推定を行ったきっかけはどうであれ、実際は、考察に記載したように対数正規分布を仮定するもっともらしい理由があるため、ポアソン分布や他の分布でパラメタ推定は行う必要性を感じなかった）
2. 集めた文(コーパスとも言える)を今後どういったことに使えるのか？
→ オートエンコーダによる教師なし学習(クラスタリング)を行いテキストマイニングを行うことが考えられる。
3. Wikipediaではなくニュースサイトなどの文の場合は違う分布になるのか？
→ 本などをみるとコーパスの分野が変わると違う分布になったり、あるコーパスで学習した機械学習のモデルが、別の分野のコーパスを入力して予測すると精度が落ちるなどの事例がある様子。そのため、ニュースサイトの分布とは異なるものになる可能性が高い。
4. XMLの読み込みにかかる時間はどれくらいか
→ 約15分くらい
5. 正規表現の適用にかかる時間はどれくらいか
→ ６コアCPUで分散処理しているため早い、約５分くらい。CPUを増やせばほぼリニアに早くなる（←発表では１０分くらいと言ったかもしれませんが訂正します）

発表スライド

発表したスライドになります。正規表現が冗長になっており(スクエアブラケットの中のバックスラッシュのこと)、もっとコンパクトにできますが、そこは認識しています。時間がないため修正せずにあげています。ご指摘ありましたらコメントをお願いいたします。
以上。

f:id:keimina:20190225193006j:plain — スライド1

f:id:keimina:20190225193114j:plain — スライド2

f:id:keimina:20190225193134j:plain — スライド3

f:id:keimina:20190225193152j:plain — スライド4

f:id:keimina:20190225193210j:plain — スライド5

f:id:keimina:20190225193306j:plain — スライド6

f:id:keimina:20190225193403j:plain — スライド7

f:id:keimina:20190225193417j:plain — スライド8

f:id:keimina:20190225193434j:plain — スライド9

f:id:keimina:20190225193448j:plain — スライド10

f:id:keimina:20190225193501j:plain — スライド11

f:id:keimina:20190225193520j:plain — スライド12

f:id:keimina:20190225193609j:plain — スライド13

2019-01-26

Pythonで２次元データ作成ーポリゴンに収まる点の集合

前回の続きです。今回は、前回に引き続き、ポリゴンに収まる点の集合を作成します。アルゴリズムは簡単で、ランダムに点を打ち、それがポリゴンの内側であれば点を残す、ポリゴンの外側であれば点を消すを繰り返すだけです。「ポリゴンの内側か？」を判定する関数を作成し、その関数に点の位置を次々と問い合わせればよいでしょう。
※ここでは「ポリゴン」は多角形の意味で使っています。

前回の記事:
keimina.hatenablog.jp

前々回の記事：
keimina.hatenablog.jp

N角形(N>=3)の領域内の点の集合を得ることができますので、前回や前々回の記事で作成した点の集合も、この記事の方法だけで作成することができるようになります。

目次
1. 点がポリゴンの領域の内側にある点か判定する方法
2. 角度の求め方
3. 角度を求めるための関数を作成する
4. 時計回りか、反時計回りかを判定する
5. 角度を求めるための関数2(反時計回り、時計回りの情報を追加)
6. 点がポリゴンの領域の内側にある点かどうかを判定する関数を作成しテストする
7. 出力結果
8. まとめ

1. 点がポリゴンの領域の内側にある点か判定する方法

点がポリゴンの領域の内側にあるか判定する方法は以下に色々書かれていました。

stackoverflow.com

上記の回答の一つを見たところ、点がポリゴンの内側にあるとき、その点とポリゴンを構成する点から算出される角度の総和は必ず360度になるらしいです。図を書いてイメージをつかみます。以下が、点A,点B,点C,,,点Jを結んだポリゴンの内部に点Pがあった時の図になります。点Bから点Cあたりでガタガタと移動する箇所があります。点Bから点Cに移動する時のように角度θiがマイナスになるときもあるのですが、その後のガタガタで相殺されます。θiの総和が360°のとき左回りに１回転できたということになります。※かなりざっくりした説明で厳密ではありません。

f:id:keimina:20190126202512p:plain — 図1. 点がポリゴンの領域の内側にあるとき

一方、点がポリゴンの外側にあるとき、360度にはならないと書かれています。これも図を書いてイメージをつかみます。θiの箇所についてみて見ると前半は右回りで、後半は左回りとなっています。右回りと、左回りで相殺されますので、θiの総和は0°になることがわかります。

f:id:keimina:20190126202605p:plain — 図2. 点がポリゴンの領域の外側にあるとき

ポリゴンのような多角形の図形は、円の一つだと考えられます。そのため多角形を円とみなして考えるとわかりやすいです。円形の壁の中に、自分がいたとして、ぐるっと見回して壁に囲まれていたら、円の中だとわかります。一方、ぐるっと見回して左には壁があったけど右には壁がなかったなどの場合は壁の中でないことがわかります。算出される角度の総和が360°だったというのは、壁をたどると左周りに１回転できたということと同じです。１回転できなかった場合は、壁を辿っても１回転できなかったということですので、壁がつながっていない、すなわち壁の内側ではないということがわかります。
今回は本アルゴリズムを使用します。(アルゴリズムがちゃんと機能するのかも確認します)

2. 角度の求め方

角度を求めるには以下の数学公式を用います。

cosθ = (PA・PB)/(|PA|*|PB|)
θ = arccos(cosθ)

※「・」は内積を、arccosはアークコサインです。
詳細は以下を参照してください。

ja.wikipedia.org

3. 角度を求めるための関数を作成する

numpyで角度を求める関数を実装すると以下になります。

コード：

import numpy as np
from numpy import inner
from numpy.linalg import norm

def angle(P, A, B):
    '''
    args
    ----
        P: np.array(x, y))
        A: np.array(x, y))
        B: np.array(x, y))
    returns
    -------
        theta: x (degree)
    '''
    # 全体をシフトさせる
    PA = A - P
    PB = B - P
    # cosθを求める
    costheta = inner(PA, PB)/norm(PA)/norm(PB)
    # 角度を求める(度)
    theta = np.arccos(costheta)/np.pi * 180
    return theta

# テスト
P = np.array((0,0))
A = np.array((1,1))
B = np.array((1,0))
print(angle(P, A, B))

出力：

45.00000000000001

4. 時計回りか、反時計回りかを判定する

3. によって角度の絶対値はもとめることができました。しかし、反時計回りか時計回りか区別できません。時計回りか、反時計回りか区別できないと、角度の総和が常に360°になってしまうため、区別できるようにする必要があります。そこで、点A(a, b)と点B(c, d)と点P(0, 0)によってできるベクトルPAとPBのなす角θについて、時計回りと反時計回りを以下のように定義します。

(i) d > (b/a) * c
ならば、反時計回り方向の角度

(ii) d < (b/a) * c
ならば時計回り方向の角度

※説明を簡単にするため点P(0, 0)としています。点Pが(0, 0)以外のときでも、全体をシフトさせて点Pが(0, 0)になるように点Aと点Bを移動させても、角度に関する情報は失われないので、P(0, 0)でも問題ありません。

また(i)(ii)の不等号の式の両辺にaをかけて、数式を変形すると以下のように書き直せます。

(i)' a*d - b*c > 0
(ii)' a*d - b*c < 0

図で表すと以下になります。点Bが緑色の領域にあればPABの角度は反時計回り方向。緑色の領域にあれば、時計回り方向というふうに定義しています。
なぜ、この定義をするのかは、説明しづらいですが、図1、図2を実際に手でなぞって見ると、まず回転方向が必要なことがわかります。また0°<θ<180°のときは足し算で-180°<θ<0°のときは引き算のように計算方法を変えなければいけないこともわかると思います。そのため、このような定義が必要であるということがわかります。

f:id:keimina:20190126202901p:plain — 図3. この記事における、反時計回り、時計回りの定義

5. 角度を求めるための関数2(反時計回り、時計回りの情報を追加)

3.で作成したangle関数を、4.の定義を使用して、反時計周りならプラスを、時計回りならマイナスの角度を返すように修正します。

コード：

import numpy as np
from numpy import inner
from numpy.linalg import norm

def angle2(P, A, B):
    '''
    args
    ----
        P: np.array(x, y))
        A: np.array(x, y))
        B: np.array(x, y))
    returns
    -------
        theta: x (degree)
    '''
    # 全体をシフトさせる
    PA = A - P
    PB = B - P
    # cosθを求める
    costheta = inner(PA, PB)/norm(PA)/norm(PB)
    # 角度を求める(ラジアン)
    theta = np.arccos(costheta)
    # 角度を求める(度)
    theta = np.arccos(costheta)/np.pi * 180
    # 
    if PA[0]*PB[1] - PA[1]*PB[0] > 0:
        # 条件(i)'に該当するときはなにもしない(ここではthetaはプラス)
        pass
    else:
        # 条件(ii)'に該当するときはthetaをマイナスにする
        theta *= -1.0
    return theta

# テスト
P = np.array((0,0))
A = np.array((1,1))
B = np.array((1,0))
print(angle2(P, A, B))

出力：

-45.00000000000001

6. 点がポリゴンの領域の内側にある点かどうかを判定する関数を作成しテストする

準備が整いましたので、点がポリゴンの領域の内側にある点かどうかを判定する関数を作成し、機能するかテストします。

点(x, y)が1.のアルゴリズムを使用しポリゴンの領域内の点であるかどうかを判定するis_inside_polygon関数を実装します。以下のようになります。
この関数は[(x1, y1), (x2, y2),,,]のような点(x, y)の集合pointsと点p(x, y)を受け取り、pointsからなるポリゴンの領域内に点pが含まれるか判定します。

コード：

import random
import matplotlib.pyplot as plt

import numpy as np
from numpy import inner
from numpy.linalg import norm

def angle2(P, A, B):
    '''
    args
    ----
        P: np.array((x, y))
        A: np.array((x, y))
        B: np.array((x, y))
    returns
    -------
        theta: float (dgree)
    '''
    # 全体をシフトさせる
    PA = A - P
    PB = B - P
    # cosθを求める
    costheta = inner(PA, PB)/norm(PA)/norm(PB)
    # 角度を求める(度)
    theta = np.arccos(costheta)/np.pi * 180
    # 時計回りか、反時計回りか判定し、thetaをプラスかマイナスにして返す
    if PA[0]*PB[1] - PA[1]*PB[0] > 0:
        # 条件(i)'に該当するときはなにもしない(ここではthetaはプラス)
        pass
    else:
        # 条件(ii)'に該当するときはthetaをマイナスにする
        theta *= -1.0
    return theta

def is_inside_polygon(points, P):
    '''
    args
    ----
        points: list of np.array((x, y))
        P: np.array((x, y))
    returns
    -------
        theta: x (degree)
    '''
    points = [np.array(p) for p in points]
    P = np.array(P)
    angle_total = 0.0
    for A, B in zip(points, points[1:] + points[0:1]): # Notice! "points[1:] + points[0:1]" means something like "points[1:].extend(points[0:1])"
        angle_total += angle2(P, A, B)
    if np.isclose(abs(angle_total), 360):
        # 360°の場合はポリゴンの内側なのでTrueを返す
        return True
    else:
        # ポリゴンの外側
        return False

# 三角形
#points = ((0,0), (1,0), (0, 1))

# リボン?
#points = ((1.0, 0.0), (0.5, 0.5), (1.0, 1.0), (-1.0, 1.0), (-0.5, 0.5), (-1.0, 0.0))

# 十字
points = ((1.0, -3.0), (1.0, -1.0), (3.0, -1.0), (3.0, 1.0), (1.0, 1.0), (1.0, 3.0), (-1.0, 3.0), (-1.0, 1.0), (-3.0, 1.0), (-3.0, -1.0), (-1.0, -1.0), (-1.0, -3.0))
points = [(i[0]/3, i[1]/3) for i in points]

N = 5000
for _ in range(N):
    # 点を打つ場所をランダムに決める(ただし、-1<x<1, -1<y<1とする)
    x = (random.random() - 0.5) * 2
    y = (random.random() - 0.5) * 2
    p = (x, y)
    # 決めた点(x, y)が三角形の内側にあればグラフに点を打つ
    if is_inside_polygon(points, p):
        plt.plot(x, y, color='k', marker='.')

# グラフを表示する
plt.gca().set_xlim(-2,2)
plt.gca().set_ylim(-2,2)
plt.gca().set_aspect('equal')
plt.show()

7. 出力結果

上から順に三角形、リボン、十字形になります。ランダムに点を打っているのですが、pointsで定義された多角形の内側にある点だけ残ることがわかります。
※形を変えるには、多角形の定義(#points... の箇所)を変えてください。

f:id:keimina:20190126203803p:plain — 図4. 三角形に収まる点の集合

f:id:keimina:20190126204014p:plain — 図5. リボン形に収まる点の集合

f:id:keimina:20190126204048p:plain — 図6. 十字形に収まる点の集合

8. まとめ

ポリゴンに収まる点の集合を作成できました。

以上です。

2019-01-20

Pythonで２次元データ作成ー三角形に収まる点の集合

前回の続きです。今回は、前回に引き続き、Pythonで三角形に収まる点の集合を作成します。アルゴリズムは簡単で、ランダムに点を打ち、それが三角形の内側であれば点を残す、三角形の外側であれば点を消すを繰り返すだけです。「三角形の内側か？」を判定する関数を作成し、その関数に点の位置を次々と問い合わせればよいでしょう。ここでは三角形は直角二等辺三角形ということにします。

前回の記事:
keimina.hatenablog.jp

1. 三角形の領域の定義

実装する前に使用する三角形の形を明確にします。以下のように使用する三角形の領域を定義します。

三角形の領域の定義：
　x-y座標上の点(x, y)で「x > 0 かつ y > 0 かつ y < -x + 1」を満たす点の集合

この条件を満たす点(x, y)は３つの点(0, 0)と(1, 0)と(0, 1)を結ぶとできる三角形の領域を構成します。なぜそうなるかについては、以下の①②③を参考にしてください。

①「x > 0 かつ y > 0」の領域は x-y座標の右上の領域すべてを表します。
②「y < -x + 1』の領域は直線 y = -x + 1 より下側の領域をすべてを表します。
③「x > 0 かつ y > 0 かつ y < -x + 1」の領域は上記①②の領域が重なった領域であり、三角形になります。

①②③を図で表すと以下になります。

f:id:keimina:20190120142434j:plain — 図1. 三角形の領域の定義

2. 三角形の領域内の点であるかどうかを判定する関数を作成する

点(x, y)が上記1.で定義したの三角形の領域内の点であるかどうかを判定する関数を実装します。以下のようになります。

def is_inside_triangle(x, y):
    if x >= 0 and \
       y >= 0 and \
       y <= -x + 1:
        # 「xが0以上の領域」かつ「yが0以上の領域」かつ「y=-x+1の直線より下の領域」
        # すなわち点(x,y)が三角形の内側ならTrue
        return True
    else:
        # 三角形の外側ならFalse
        return False

is_inside_triangle関数は引数で与えられた点(x, y)の情報が1.で定義した三角形の内側にあるかどうかを判定します。三角形の内側であればTrue、それ以外はFalseを戻します。

3. コード1(ライブラリを極力使わない場合)

2.で作成した関数を使用することで、三角形に収まる点の集合を求めることができます。ライブラリを極力つかわないで全体のコードを書くとこんな感じになるかと思います。前回の記事のコードにある、is_inside_square関数を今回作成したis_inside_triangle関数に置き換えただけです。

※matplotlibは可視化のため使用

import random
import matplotlib.pyplot as plt

def is_inside_triangle(x, y):
    if x >= 0 and \
       y >= 0 and \
       y <= -x + 1:
        # 「xが0以上の領域」かつ「yが0以上の領域」かつ「y=-x+1の直線より下の領域」
        # すなわち点(x,y)が三角形の内側ならTrue
        return True
    else:
        # 三角形の外側ならFalse
        return False

N = 1000
for _ in range(N):
    # 点を打つ場所をランダムに決める
    x = random.random()
    y = random.random()
    # 決めた点(x, y)が三角形の内側にあればグラフに点を打つ
    if is_inside_triangle(x, y):
        plt.plot(x, y, color='k', marker='.')

# グラフを表示する
plt.gca().set_xlim(-2,2)
plt.gca().set_ylim(-2,2)
plt.gca().set_aspect('equal')
plt.show()

4. コード2(Numupyを使用して書いた場合)

Numpyを使用して書くと以下のようになると思います。これも、前回の記事のコードにある、is_inside_square関数を今回作成したis_inside_triangle関数に置き換えただけです。

import numpy as np
import matplotlib.pyplot as plt

def is_inside_triangle(x, y):
    if x >= 0 and \
       y >= 0 and \
       y <= -x + 1:
        # 「xが0以上の領域」かつ「yが0以上の領域」かつ「y=-x+1の直線より下の領域」
        # すなわち点(x,y)が三角形の内側ならTrue
        return True
    else:
        # 三角形の外側ならFalse
        return False

# 関数をnumpy化する
is_inside_triangle = np.vectorize(is_inside_triangle)

N = 1000

# ランダムに点を打つ場所を決める
x = np.random.random(N)
y = np.random.random(N)

# 三角形の内側にある点だけを抜き出す
conditions = is_inside_triangle(x, y)
x = x[conditions]
y = y[conditions]

# グラフに点を打つ
plt.plot(x, y, color='k', marker='.')

# グラフを表示する
plt.gca().set_xlim(-2,2)
plt.gca().set_ylim(-2,2)
plt.gca().set_aspect('equal')
plt.show()

5. 出力

以下のようなグラフが出力されると思います。Pythonで三角形に収まる点の集合を作成できました。

f:id:keimina:20190120145112p:plain — 図2. 出力

以上です。

ae+bg	af+bh
ce+dg	cf+dh

1. はじめに

2. 使用するモジュールと、そのインストール

3. ニューラルネットワークのモデルと、順伝播の式

4. ソースコード

5. 出力結果

6. まとめ

7. 参考文献

Numpy 入門 (機械学習勉強会 in 新潟 2019/03/23)

目次

1 はじめに

2 numpy とは何か

3 numpy モジュールのインポート方法について

4 数値を作成する

4.1 np.int32 とは

5 「１次元配列」の数値計算をする

5.1 「１次元配列」を作成する

5.1.1 配列の'要素'とは

5.1.2 np.array 関数

5.1.3 np.arange 関数

5.2 「１次元配列」の形状 shape について

5.3 「１次元配列」 + 「１次元配列」

5.4 配列の要素に数学関数を適用する

5.5 １次元配列の要素の合計を求める

6 「２次元配列」の数値計算をする

6.1 「２次元配列」を作成する

6.1.1 np.array 関数

6.1.2 reshape メソッド

6.2 「２次元配列」 + 「２次元配列」

6.3 「２次元配列」の要素の合計を求める

6.3.1 列の合計

6.3.2 行の合計

6.3.3 配列の次元数を保つには

6.4 行列計算用の関数

7 配列を作成する関数

7.1 固定された値を要素にもつ配列を作成する関数

7.2 単位行列を作成する関数

7.3 ランダムな値を要素にもつ配列を作成する関数

7.4 似たような配列を作成する関数

7.5 ある範囲の数値を要素にもつ配列を作成する関数

8 配列の要素へのアクセス

8.1 配列の要素の取得

8.1.1 要素の取得の基本

8.1.2 要素を範囲で指定して取得

8.1.3 逆順の配列の取得

8.2 配列の要素の書き換え

8.2.1 要素の書き換えの基本

8.2.2 要素を範囲で指定して書き換え

9 真偽値(True/False)について

10 配列の制約

11 配列の長さや、形状が異なる場合に行われる計算について

11.1 エラーになるか計算結果が求まるかの判定方法

11.2 Broadcasting rule の適用例

12 速度比較

13 まとめ

14 最後に

目次

1. 動作環境について

2. Wikipedia 文書のアーカイブについて

3. 発表スライド

4. やってみた感想

5. 発表会で上がった質問や意見

発表スライド

目次

1. 点がポリゴンの領域の内側にある点か判定する方法

2. 角度の求め方

3. 角度を求めるための関数を作成する

4. 時計回りか、反時計回りかを判定する

5. 角度を求めるための関数2(反時計回り、時計回りの情報を追加)

6. 点がポリゴンの領域の内側にある点かどうかを判定する関数を作成しテストする

7. 出力結果

8. まとめ

目次

1. 三角形の領域の定義

2. 三角形の領域内の点であるかどうかを判定する関数を作成する

3. コード1(ライブラリを極力使わない場合)

4. コード2(Numupyを使用して書いた場合)

5. 出力

Numpy 入門
(機械学習勉強会 in 新潟 2019/03/23)