三面体のサイコロを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 の範囲となると推測しました。 この範囲外にあった場合、それは割とよく発生するものでない(異常な状態)ということがわかります。ただし、この推測が役に立つのは、あくまで仮定が成立する時の話で、この仮定とデータが何ら紐づかない場合は実際あまり意味のない議論になります。しかし、データと仮定が紐づかない場合でも、このような条件が議論する上で重要であるということがわかっているのであれば、何の指標も方向性も定まっていない状態より良いのではないかと思います。

以上です。

Pandas の DataFrame をエクセルで開く自作関数

表データを可視化したい時、まずエクセルで生のデータを表示したい時があります。これは具体的には「Pandas の DataFrame を Excel で保存し、それを開く」という作業です。この作業を自動化したら、便利だったので、ここに紹介します。

1. 動作環境

OS MacOS Mojave
Python Anaconda 1.6.14 (Python 3.6)

※少なくとも pandas と xlwings が必要です。

2. 作成した関数

def show_excel(df, out_filename='tmp.xlsx'):
    '''
    Parameters
    ----------
        df:
            DataFrame to show

        out_filename:
            Filename of the DataFrame(optional) default 'tmp.xlsx'
    
    Returns:
        None
    '''
    df.to_excel(out_filename)
    app = xw.App(visible=None, add_book=True)
    wb = xw.books.open(out_filename)
    wb.activate(steal_focus=True)

3. 作成した関数の説明

以下を行う関数です。

1. 引数で指定した DataFrame をカレントディレクトリに tmp.xlsx として保存

2. tmp.xlsx を開いて、最前面に持ってくる

4. 実行方法

import pandas as pd
import numpy as np
import xlwings as xw

df = pd.DataFrame(np.arange(15).reshape(5, 3), columns=['x', 'y', 'z'])
show_excel(df)

# エクセルのプロセスを全て終了(kill) する
# [app.quit() for app in xw.apps]

5. 実行結果

実行すると、 エクセルが起動し df の内容が表示されます。

f:id:keimina:20190922105956p:plain
実行結果

6. 注意事項

  • デフォルトでカレントディレクトリの 'tmp.xlsx' が上書き保存されますので注意してください

  • エクセルの新しいプロセスが起動します
  • 起動したエクセルを終了する場合は、プロセスを終了(kill)してください(実行方法は使用方法の一番下のコメントに記載)


以上です。

(勉強会用資料)Python Pandas もくもく勉強会 in 新潟 #3

1. はじめに

新潟で Python Pandas もくもく勉強会 というものを主催しています。本記事はその 2019/08/21 に行われる勉強会で私が行う 5分LT の資料です。1週間くらい Pandas を使ってみて、記憶に残っているものをまとめました。

connpass.com

2019/08/21 の Pandas 勉強会用
updated: <2019-08-20 Tue>

Pandas 覚えたことまとめ

2. import 文

import pandas as pd
import numpy as np

3. 出現頻度のカウント

data = ['YES', 'YES', 'NO', 'NO', '-']
s = pd.Series(data)
# series の value_counts メソッドを使用する
s.value_counts()

# NO     2
# YES    2
# -      1
# dtype: int64

# DataFrame の場合
data = ['YES', 'YES', 'NO', 'NO', '-']
data2 = ['A', 'A', 'A', 'B', 'C']
s2 = pd.Series(data2)
df = pd.concat([s, s2], axis=1)
# df に value_counts メソッドは存在しない
# df の列全体に value_counts を適用するときは apply を使用するとできる
df.apply(lambda s: s.value_counts())

#        0    1
# -    1.0  NaN
# A    NaN  3.0
# B    NaN  1.0
# C    NaN  1.0
# NO   2.0  NaN
# YES  2.0  NaN

4. Excel あるある

4.1. 作者の主観により特定の行がデータの区切りとなっていると思われるエクセルをどうにかする

data = ['・赤い果物', 'いちご', 'りんご', '・緑色の果物', 'キウイ', 'メロン', '・紫色の果物', 'ぶどう']
s = pd.Series(data)
s

# 0     ・赤い果物
# 1       いちご
# 2       りんご
# 3    ・緑色の果物
# 4       キウイ
# 5       メロン
# 6    ・紫色の果物
# 7       ぶどう
# dtype: object

# DataFrame に変換(あとで pd.merged で結合できるようにするため)
df = s.to_frame()

# '・赤い果物', '・緑色の果物', '・紫色の果物' に共通の'・'を含む行を抽出、name属性を設定しておく
s2 = s[s.map(lambda x: '・' in x)]
s2.name = '種類'
s2

# 0     ・赤い果物
# 3    ・緑色の果物
# 6    ・紫色の果物
# Name: 種類, dtype: object

# DataFrame に変換(あとで pd.merged で結合できるようにするため)
df2 = s2.to_frame()

# インデックスをキーとして結合
df3 = pd.merge(df2, df, left_index=True, right_index=True, how='outer')
df3

#        種類       0
# 0   ・赤い果物   ・赤い果物
# 1     NaN     いちご
# 2     NaN     りんご
# 3  ・緑色の果物  ・緑色の果物
# 4     NaN     キウイ
# 5     NaN     メロン
# 6  ・紫色の果物  ・紫色の果物
# 7     NaN     ぶどう

# 「種類」列の NaN を 埋める
df3['種類'].fillna(method='ffill', inplace=True)
df3

#        種類       0
# 0   ・赤い果物   ・赤い果物
# 1   ・赤い果物     いちご
# 2   ・赤い果物     りんご
# 3  ・緑色の果物  ・緑色の果物
# 4  ・緑色の果物     キウイ
# 5  ・緑色の果物     メロン
# 6  ・紫色の果物  ・紫色の果物
# 7  ・紫色の果物     ぶどう

# df3 から s2 のインデックスの行を削除する
df3.drop(index=s2.index, inplace=True)
df3

#        種類    0
# 1   ・赤い果物  いちご
# 2   ・赤い果物  りんご
# 4  ・緑色の果物  キウイ
# 5  ・緑色の果物  メロン
# 7  ・紫色の果物  ぶどう

(勉強会用資料)Python Pandas もくもく勉強会 in 新潟 #2

1. はじめに

新潟で Python Pandas もくもく勉強会 というものを主催しています。本記事はその 2019/08/04 に行われる勉強会で私が行う 5分LT の資料です。(資料置き場がないのでブログを活用しています)。1週間くらい Pandas を使ってみて、記憶に残っているものをまとめました。

connpass.com

2. import 文

import pandas as pd

3. Series または DataFrame の index を1始まりにする方法

3.1. Series の場合

3.1.1. 使用するデータ
data = ['a', 'b', 'c']
s = pd.Series(data)
s
#################
# In [392]: s   #
# Out[392]:     #
# 0    a        #
# 1    b        #
# 2    c        #
# dtype: object #
#################
3.1.2. 方法
'''index を指定しない場合、自動で 0 始まりのインデックスとなるので単に +1 すればよい'''
s.index += 1
s
#################
# In [89]: s    #
# Out[89]:      #
# 1    a        #
# 2    b        #
# 3    c        #
# dtype: object #
#################

3.2. DataFrame の場合

3.2.1. 使用するデータ
data = [['a', 'b', 'c'],
        ['d', 'e', 'f'],
        ['g', 'h', 'i']]
df = pd.DataFrame(data)
df
################
# In [397]: df #
# Out[398]:    #
#    0  1  2   #
# 0  a  b  c   #
# 1  d  e  f   #
# 2  g  h  i   #
################
3.2.2. 方法
'''(Series と同じ)index を指定しない場合、自動で 0 始まりのインデックスとなるので単に +1 すればよい'''
df.index += 1
df
################
# In [189]: df #
# Out[189]:    #
#    0  1  2   #
# 1  a  b  c   #
# 2  d  e  f   #
# 3  g  h  i   #
################

4. Series または DataFrame を比較して一致する行を探す

4.1. Series の場合

4.1.1. 使用するデータ
data = ['a', 'b', 'c']
data2 = ['a', 'x', 'c']
s = pd.Series(data)
s2 = pd.Series(data2)
s
s2
#################
# In [128]: s   #
# Out[128]:     #
# 0    a        #
# 1    b        #
# 2    c        #
# dtype: object #
#               #
# In [129]: s2  #
# Out[129]:     #
# 0    a        #
# 1    x        #
# 2    c        #
# dtype: object #
#################
4.1.2. 比較を行う
s3 = s==s2
s3
#################
# In [130]: s3  #
# Out[130]:     #
# 0     True    #
# 1    False    #
# 2     True    #
# dtype: bool   #
#################
4.1.3. フィルタリング
s4 = s[s3]
s4
################
# In [380]: s4 #
# Out[380]:    #
# 0    a       #
# 2    c       #
################

4.2. DataFrame の場合

4.2.1. 使用するデータ
data = [['a', 'b', 'c'],
        ['d', 'e', 'f'],
        ['g', 'h', 'i']]

data2 = [['a', 'b', 'c'],
         ['d', 'x', 'f'],
         ['g', 'h', 'i']]

df = pd.DataFrame(data)
df2 = pd.DataFrame(data2)
df
df2
#################
# In [227]: df  #
# Out[227]:     #
#    0  1  2    #
# 0  a  b  c    #
# 1  d  e  f    #
# 2  g  h  i    #
#               #
# In [228]: df2 #
# Out[228]:     #
#    0  1  2    #
# 0  a  b  c    #
# 1  d  x  f    #
# 2  g  h  i    #
#################
4.2.2. 一致する行を見つける方法
4.2.2.1. 比較
df3 = df==df2
df3
########################
# In [244]: df3        #
# Out[244]:            #
#       0      1     2 #
# 0  True   True  True #
# 1  True  False  True #
# 2  True   True  True #
########################
4.2.2.2.フィルタリング
# df と df2 で一致する行を抽出したい。
# df3 で全ての行が True の行のことなので、その行を抽出する
# all を使用するとできる
# all で reduce する方向は axis で決める、列方向に reduce するには axis=1 とする
# all は全ての行(観測点)に適用される
df4 = df3.all(axis=1)
df4
#################
# In [253]: df4 #
# Out[255]:     #
# 0     True    #
# 1    False    #
# 2     True    #
# dtype: bool   #
#################

# 以下のようにすると df と df2 の行において一致するの行だけ抽出できる
df5 = df[df4]
df5
#################
# In [280]: df5 #
# Out[280]:     #
#    0  1  2    #
# 0  a  b  c    #
# 2  g  h  i    #
#################
4.2.3. 一致しない行を見つける方法
4.2.3.1. 比較
df6 = df!=df2
df6
##########################
# In [288]: df6          #
# Out[288]:              #
#        0      1      2 #
# 0  False  False  False #
# 1  False   True  False #
# 2  False  False  False #
##########################
4.2.3.2. フィルタリング
# 注意!!: ひとつでも値が異なる行を見つけたいので、 all ではなく any を使用する
# (個人的にはまったのでメモ)
df7 = df6.any(axis=1)
df7
#################
# In [319]: df7 #
# Out[319]:     #
# 0    False    #
# 1     True    #
# 2    False    #
# dtype: bool   #
#################

df8 = df[df7]
df8
#################
# In [333]: df8 #
# Out[333]:     #
#    0  1  2    #
# 1  d  e  f    #
#################

# (Tips) 上記以外の書き方で、等価となる式
# 「"一致する行"でない行」を抽出
df9 = df[~df4]
df9
#################
# In [421]: df9 #
# Out[421]:     #
#    0  1  2    #
# 1  d  e  f    #
#################

5. np.full みたいに要素全てが同じ値からなる Series を作る方法

data = ['a']*5
s = pd.Series(data)
s
#################
# In [387]: s   #
# Out[387]:     #
# 0    a        #
# 1    a        #
# 2    a        #
# 3    a        #
# 4    a        #
# dtype: object #
#################

(エクセルで作成した csv ファイルについて)一つのセルに書かれている複数の行を、テキストエディタ Emacs で一つの行で見れるようにする方法

1. はじめに

エクセルで作成した csv ファイルで一つのセルの中に複数の行の記載があるファイルは、よく見かけると思います。今回の記事では、この csv ファイルをテキストエディタ Emacs で開き、一つのセルに書かれている複数の行を一つの行で見れるように変換する方法を書きます。さらに、変換したあと逆変換し元に戻す方法も書きます。使用するエディタは Emacs です。

2. 注意事項、制限

以下の csv では正しく動作しません。
- CRLF(キャリッジリターン + ラインフィード)以外で、 CR(キャリッジリターン)を使用している csv ファイル

※変換、逆変換するときに CR を区切り文字として使用するため、 CR を単独で使用している csv ファイルの場合、正しく動作しません。

3. 入力ファイル

3.1. エクセル

以下の様なデータをエクセルで作成します。2行目のA列のセルと2行目のB列のセルが複数行になっています。これを csv で保存したものが入力ファイルです。

f:id:keimina:20190731224004p:plain

3.2. 入力ファイルの初期状態

入力ファイルを Emacs で開くと以下の様に表示されます。^M と表示されているのは CR(キャリッジリターン)を表す文字です。LF(ラインフィード)は改行として表示されています。

f:id:keimina:20190731224102p:plain

4. 方法

複数行から1行への変換、1行から複数行へ変換する具体的な方法は以下の通りです。Emacscsv ファイルを開き以下のような手順を実行すると複数行、CR(キャリッジリターン)区切りの1行にできます。理論については「5. 背景、理論」に記載しました。

4.1. 1行への変換

複数行から1行へ変換するには以下を行います。

4.1.1. 検索する正規表現

検索する正規表現を入力します。

f:id:keimina:20190731224137p:plain

4.1.2. 置換のための emacs lisp

置換するための文字列として emacs lisp 式を入力します。

f:id:keimina:20190731224208p:plain

4.1.3. 置換実行後

置換実行後の状態です。複数行になっていた箇所が、CR区切りでスッキリと1行にまとまりました。

f:id:keimina:20190731224234p:plain

4.2. 複数行への変換

1行から複数行へ変換(逆変換)するには以下を行います。

4.2.1. 検索する正規表現

検索する正規表現を入力します。

f:id:keimina:20190731224349p:plain

4.2.2. 置換のための emacs lisp

置換するための文字列として emacs lisp 式を入力します。

f:id:keimina:20190731224419p:plain

4.2.3. 置換実行後

置換実行後の状態です。入力ファイルの初期状態に戻りました。

f:id:keimina:20190731224445p:plain

5. 背景、理論

エクセルで作成した csv ファイルは、セルの中の改行は LF で、そうでない改行(個々のサンプルを隔てるための改行)は CRLF で表されるようです。今回は LF で記載された複数の行を一つにすることが目的なので、Emacs に表示されている CRLF を普通の改行に、普通の改行を CR に変換すればよさそうだと考えました。すなわち CRLF を LF に LF を CR に同時に変換するということです。そして、そうやってできたデータを元に戻す、逆の変換には、LF を CRLF に CR を LF に変換すればよさそうです。このような変換は、一般的にはプログラミング言語を使うと正規表現のグルーピングを用いた置換で実現できると思います。テキストエディタ Emacs では、正規表現で指定した文字列を置換するための文字列として emacs lisp 式を使用することができます(「バックスラッシュ+カンマ+式」の形式で記載します)ので、こちらを使用して置換しました。

6. まとめ

テキストエディタ Emacs を使用して、csv で一つのセルに書かれている複数の行を CR 区切りで一つの行として見れるように変換しました。さらに、それを逆変換し元に戻しました。

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

Matplotlib でお絵かきアプリっぽいものを作る

Matplotlib でお絵かきアプリっぽいものを作ります。

1. matplotlib のバージョン

GUIAPIの仕様は変わる可能性が高いため、matplotlib のバージョンを記載しておきます。
バージョン 2.2.2 で動作確認済みです。以下のようにするとバージョンの確認ができます。

import matplotlib
print(matplotlib.__version__)
# Out: '2.2.2'

2. 仕様

アプリの仕様は以下とします。

No 状態 動作
1 マウス左ボタン プレス ペンを下ろす
2 マウス左ボタン リリース ペンを上げて、軌跡を消去
3 マウス左ボタンをプレス中にマウスを動かす ペンで軌跡を描画
4 マウス左ボタンをリリース中にマウスを動かす なにもしない

簡単にいってしまうと、マウス左ボタンを押している間だけマウスの軌跡を描画するアプリです。ボタンを離すと軌跡は消えます。一筆書きをしている間は軌跡をみることができます。

3. コード

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Circle, Arrow

fig, ax = plt.subplots()

being_pressed = False
pre_x = None
pre_y = None

def on_motion(event):
    global pre_x, pre_y
    if None in (event.xdata, event.ydata): # マウスが画面の外にいるときは何もしない
        #print('x=%d, y=%d, xdata=%f, ydata=%f' %(event.x, event.y, event.xdata, event.ydata))
        return
    if being_pressed:
        ax.add_patch(Circle((event.xdata, event.ydata), 0.01, color='r'))
        if None not in (pre_x, pre_y): # pre_x と pre_y が存在するとき矢印を描画
            ax.add_patch(
                Arrow(pre_x, pre_y, event.xdata-pre_x, event.ydata-pre_y,
                      linestyle='solid', color='blue', width=0.02))
        pre_x = event.xdata
        pre_y = event.ydata
        fig.canvas.draw() # 必須!

def on_button_press(event):
    print('button pressed')
    global being_pressed
    being_pressed = True

def on_button_release(event):
    print('button released')
    global being_pressed, pre_x, pre_y
    being_pressed = False
    pre_x = None
    pre_y = None
    # reset patches
    for patch in reversed(ax.patches):
        patch.remove()
    fig.canvas.draw() # 必須!


cid_1 = fig.canvas.mpl_connect('motion_notify_event', on_motion)
cid_2 = fig.canvas.mpl_connect('button_press_event', on_button_press)
cid_3 = fig.canvas.mpl_connect('button_release_event', on_button_release)

plt.show()

4. コードの説明

マウスの左ボタンを押しているかどうかを being_pressed 変数で管理することとしています。pre_x, pre_y はボタンが押されているときの、マウスを動かす前のマウスの座標を保存するために定義しています。(pre_x, pre_y はボタンが押されてないときは存在しないので None で初期化します)

ボタンが押されたときに独自関数(on_button_press関数)を呼んでほしいため

fig.canvas.mpl_connect('button_press_event', on_button_press)

としています。これにより canvas で button_press_event が発生すると on_button_press 関数が呼ばれるようになります。同様に、 ボタンが離されたとき(button_release_event イベントが発生したとき)には on_button_release、 マウスが動いたとき(motion_notify_event イベントが発生したとき)には on_motion 関数が呼ばれるように設定します。

話はそれますが、matplotlib で他にどのようなイベントが提供されているか興味のある方は、以下を参照してください。
https://matplotlib.org/users/event_handling.html より抜粋

Event name Class and description
'button_press_event' MouseEvent - mouse button is pressed
'button_release_event' MouseEvent - mouse button is released
'draw_event' DrawEvent - canvas draw (but before screen update)
'key_press_event' KeyEvent - key is pressed
'key_release_event' KeyEvent - key is released
'motion_notify_event' MouseEvent - mouse motion
'pick_event' PickEvent - an object in the canvas is selected
'resize_event' ResizeEvent - figure canvas is resized
'scroll_event' MouseEvent - mouse scroll wheel is rolled
'figure_enter_event' LocationEvent - mouse enters a new figure
'figure_leave_event' LocationEvent - mouse leaves a figure
'axes_enter_event' LocationEvent - mouse enters a new axes
'axes_leave_event' LocationEvent - mouse leaves an axes

話を戻して、以下のコードで現在のマウスのx, y座標に丸を描画します。0.01 は 丸の大きさです。何度か試して 0.01 が良さそうだったのでそうしています。

ax.add_patch(Circle((event.xdata, event.ydata), 0.01, color='r'))

以下のコードで、前のマウスのx, y座標と現在の座標に矢印を描画します。0.02 は線の太さです。これも試した結果これくらいがよさそうだったのでそうしています。

ax.add_patch(
    Arrow(pre_x, pre_y, event.xdata-pre_x, event.ydata-pre_y,
    linestyle='solid', color='blue', width=0.02))

これら 0.01 や 0.02 を計算でだせるようにするのは私の宿題としておきます。

5. 実行結果

f:id:keimina:20190705002816g:plain
matplotlibによるお絵かきアプリっぽいもの

6. まとめ・考察

matploblib でお絵かきアプリっぽいものを作りました。
マウスがクリックされたら〜するとか、マウスが動いたら〜するなどいろいろ応用できるのではないかと思います。筆跡を時系列データとして機械学習して分類するのもおもしろそうです、夏休みの自由研究にもできそうではないでしょうか。

7. 今後の予定

個人的にはこれを使って軌跡に囲まれた領域のデータの点だけを取得して、それらのデータの可視化したいです(実はもうしたけどそれは次回以降、書きます)。

8. 宿題

matplotlib の Arrow の太さとか Circle の大きさを計算できるように仕様を把握する

それでは、おやすみなさい。

勾配降下法により PCA(主成分分析)の第2主成分を数値解析的に求め、PCAについて理解を深める

1. はじめに

前回記事で、第2主成分は「第1主成分と直交する&分散が最大となる結合係数を求めるということだと理解しています」と書きました。本記事は、その理解で正しいか検証したときの記録です。検証した結果、その理解で正しいということがわかりました。

2. 前回記事

keimina.hatenablog.jp

3. サンプルデータについて

前回同様、使用するサンプルデータは 3次元の特徴量をもつ3つのデータとします。
現実にこのサイズのデータで PCA を行うことは稀かと思われますが、PCA の理解には十分だと思います。

4. 実装の方針

第2主成分が第1主成分と直交するというのは第1主成分をベクトル Wa 、 第2主成分をベクトル Wb とすると Wa と Wb の内積が 0 になるということです。 Wa と Wb が以下の時

Wa = [Wa1, Wa2, Wa3]
Wb = [Wb1, Wb2, Wb3]

内積は以下のように定義されます。

inner_product = Wa1*Wa2 + Wa2*Wb2 + Wa3*Wb3


Wa と Wb が直交しているとき、内積が 0 なので

Wa1*Wb1 + Wa2*Wb2 + Wa3*Wb3 = 0

を Wa, Wb は満たします。ちなみに前回記事で第1主成分 Wa は求めることができましたので、 Wa1, Wa2, Wa3 には数値が入ることがわかります。上記の式を変形して Wb1 について解くと

Wb1 = -(Wa2*Wb2 + Wa3*Wb3) / Wa1

となります。 これを Wb の Wb1 に代入すると

Wb = [-(Wa2*Wb2 + Wa3*Wb3) / Wa1, Wb2, Wb3]

と書くことができます(細かいですが Wa1が 0 でないときに限ります)。Wa と Wb が直交するという条件があるので、それを反映させたときの Wb がこの式になります。この Wb をサンプルデータに適用し適用後のデータ間の分散が最大になる Wb が第2主成分であると理解していますので、その Wb と sklearn.decomposition.PCA を使用して求めた第2主成分を比較し、同じ方向を向いているか検証します。

※当たり前ですが Wb を求めるときに Wa2, Wa3 すなわち Wa の情報が必要です。このように依存関係がありますのでコーディング時に注意します。

5. ソースコード

ソースコードについてですが、前回は行列のシンボルとして sympy.MatrixSymbol を使用していたのですが、今回は sympy.Symbol と np.array を使うことにしました。理由は、前回記事のソースコードでは、Matrixによるキャストが大量発生してしまい、期待していたほど可読性が上がるわけではなかったからです。(MatrixSymbol を使用すると行列の各要素の変数を自動的に作成してくれて便利だったのですが残念です)そのため本記事では MatrixSymbol の代わりに、 Symbol と np.array を使うことにしました。

# -*- coding:utf-8 -*-
# MatrixSymbol を使わないで np.array で統一した & Second component 追加したバージョン

import numpy as np
from sympy import symbols, diff, Eq, solve
from sklearn.decomposition import PCA
from itertools import chain

# 1行目が観測データ1、2行目が観測データ2、3行目が観測データ3
# 1列目が特徴量1、2列目が特徴量2、3列目が特徴量3 とする
data = np.array([[3,5,7],[11,13,17],[19,23,29]], dtype=np.float64)

################################################################
# 独自手法で
# first component (第一主成分)を求めるための処理
################################################################

# 観測データ入力用の 3x3 の行列のシンボルを作成
# 1列目が特徴量1、2列目が特徴量2、3列目が特徴量3
Data = np.array([[symbols('Data%d%d'%(i, j)) for j in range(1, 4)] for i in range(1, 4)])

# first component 用のシンボルを作成(サイズ=3x1)
Wa = np.array([[symbols('Wa%d1'%i)] for i in range(1, 4)])

# second component 用のシンボルを作成(サイズ=3x1)
Wb = np.array([[symbols('Wb%d1'%i)] for i in range(1, 4)])

# second component は first component と直行する(内積=0となる)条件を使用し、1変数を置き換える
inner_prod = Wa.T.dot(Wb)[0][0]
# solve 関数を使用して Wa と 直行する Wb の Wb[0][0] を求める
# solve を使用するとリストが戻るため、要素を取り出すため インデックスを指定して取り出す
new_Wb11 = solve(Eq(inner_prod), Wb[0, 0])[0]
# 全要素の数式に含まれる Wb11 を new_Wb11 と置き換える関数を定義
replace_Wb11 = np.vectorize(lambda x: x.subs({Wb[0, 0]: new_Wb11}))

# Data を以下のように変換し、次元を1次元に削減する
Y = Data.dot(Wa) # shape=(3, 1)
Y2 = Data.dot(Wb)

# Y の特徴量の分散を求める
Ymean = np.mean(Y, axis=0, keepdims=True)
Yhat = (Y - Ymean)   # 重心を (0, 0, 0) にする、shape=(3, 1)
V = Yhat.T.dot(Yhat) # 分散を求める、 shape=(1, 1)

# Y2 特徴量の分散を求める
Y2mean = np.mean(Y2, axis=0, keepdims=True)
Y2hat = (Y2 - Y2mean)   # 重心を (0, 0, 0) にする、shape=(3, 1)
V2 = Y2hat.T.dot(Y2hat) # 分散を求める、 shape=(1, 1)

# 分散 V が最大になるようにする
# ここでは、数値解析的に勾配降下法を使用しVの最大値を求めるため微分する
diff_V = np.vectorize(lambda x: diff(V[0, 0], x))
diff_V2 = np.vectorize(lambda x: diff(V2[0, 0], x))

dVdWa = diff_V(Wa) #shape=(3, 1)
dV2dWb = diff_V2(Wb) #shape=(3, 1)

# first components の値を保持する変数を用意しランダムな数値で初期化
np.random.seed(51) # 他の人がやっても同じ結果になるようにする
Wa_val = np.random.random((3,1))
Wb_val = np.random.random((3,1))

# ハイパーパラメータ(学習係数)を定義
b = 0.05

# 勾配を求める関数を定義
dVdWa_func = lambda Wa_val: np.array([w[0].subs(dict(zip(chain(Data.flatten(), Wa.flatten()),
                                       chain(data.flatten(), Wa_val.flatten()))))
                                      for w in dVdWa], dtype=np.float64).reshape((3, 1))
dV2dWb_func = lambda Wb_val, Wa_val: np.array([w[0].subs(dict(zip(chain(Data.flatten(), Wb.flatten(), Wa.flatten()),
                                       chain(data.flatten(), Wb_val.flatten(), Wa_val.flatten()))))
                                               for w in replace_Wb11(dV2dWb)], dtype=np.float64).reshape((3, 1))

print('===========================================')
print('================ 第1主成分 ================')
print('===========================================')
# 第1主成分を学習する
for i in range(101):
    # 勾配を求める
    dVdWa_val = dVdWa_func(Wa_val)
    # ハイパーパラメータ b をかけて Wa_val を更新する
    Wa_val += b*dVdWa_val
    # 正規化する
    Wa_val /= np.linalg.norm(Wa_val)
    if i%10==0:
        print('================ 第1主成分 学習 i = ', i, "================")
        v = V[0, 0].subs(dict(zip(chain(Data.flatten(), Wa.flatten()),
                                       chain(data.flatten(), Wa_val.flatten()))))
        print("分散V = ", v)
        print("Wa = ", Wa_val.T)



print('===========================================')
print('================ 第2主成分 ================')
print('===========================================')

# 第2主成分を学習する
calc_Wb = np.vectorize(lambda x: x.subs(dict(zip(chain(Wb.flatten(), Wa.flatten()),
                                                 chain(Wb_val.flatten(), Wa_val.flatten())))))
for i in range(101):
    # Wb を Wa と直交するようにする
    Wb_val = calc_Wb(replace_Wb11(Wb)).astype(np.float64)
    # 勾配を求める
    dV2dW_val = dV2dWb_func(Wb_val, Wa_val)
    # ハイパーパラメータ b をかけて Wb_val を更新する
    Wb_val += b*dV2dW_val
    # 正規化する
    Wb_val /= np.linalg.norm(Wb_val)
    if i%10==0:
        print('================ 第2主成分 学習 i = ', i, "================")
        v2 = V2[0, 0].subs(dict(zip(chain(Data.flatten(), Wb.flatten(), Wa.flatten()),
                                       chain(data.flatten(), Wb_val.flatten(), Wa_val.flatten()))))
        print("分散V2 = ", v2)
        print("Wb = ", Wb_val.T)

################################################################
# sklearn.decomposition.PCA で
# first component (第一主成分)を求めるための処理
################################################################
# solving by sklearn.decomposition.PCA
pca = PCA(n_components=3) # 1-3 components を求める
pca.fit(data) # 学習する

################################################################
# 独自手法 と sklearn.decomposition.PCA で
# 求めた first component (第一主成分)を表示するための処理
################################################################
print('\n================ 検証 ================')
# 独自手法 で求めた first component を表示(比較しやすいようにflattendを使用)
print('--- 独自手法 ---')
print('第1主成分 = ', Wa_val.flatten())
print('第2主成分 = ', Wb_val.flatten())

# sklearn.decomposition.PCA で求めた frist component を表示
w1_pca, w2_pca, w3_pca = pca.components_
print('--- sklearn.decomposition.PCA ---')
print('第1主成分 = ', w1_pca)
print('第2主成分 = ', w2_pca)

# sklearn.decomposition.PCA で求めた frist component と
# コサイン類似度を求める、w1_pca と Wa_val はノルムが1なので内積を
# とるだけで cosθ が求まる
print('--- コサイン類似度 ---')
cos_theta1 = w1_pca.dot(Wa_val)
cos_theta2 = w2_pca.dot(Wb_val)
print(' 第1主成分 cosθ = {} ({}°)'.format(*cos_theta1, *np.rad2deg(np.arccos(cos_theta1))))
print(' 第2主成分 cosθ = {} ({}°)'.format(*cos_theta2, *np.rad2deg(np.arccos(cos_theta2))))

# 第1主成分と第2主成分における分散の比率
print('--- 第1主成分と第2主成分における分散の比率(sklearn.decomposition.PCA) ---')
print('pca.explained_variance_ratio_[0] / pca.explained_variance_ratio_[1] \n =', pca.explained_variance_ratio_[0] / pca.explained_variance_ratio_[1])
print('--- 第1主成分と第2主成分における分散の比率(独自手法) ---')
print('V / V2 \n = ', v / v2)

6. ソースコードの説明

ソースコードの説明を以下に記載します。

diff_V = np.vectorize(lambda x: diff(V[0, 0], x))
diff_V2 = np.vectorize(lambda x: diff(V2[0, 0], x))

dVdWa = diff_V(Wa) #shape=(3, 1)
dV2dWb = diff_V2(Wb) #shape=(3, 1)

上記は、V の Wa に関する微分と V2 の Wb に関する微分を求めています。微分を行う関数である diff 関数をブロードキャストできるようにするため np.vectorize 関数を使用しています。

new_Wb11 = solve(Eq(inner_prod), Wb[0, 0])[0]

上記は、 Eq(inner_prod) で inner_prod = 0 となる数式を構築し、 solve(数式, Wb[0, 0]) で 数式をWb[0, 0] について解いています。solve 関数は解が複数ある場合に備えてリストで解を戻します。今回は、数式をみれば解が一つであることは自明なので、インデックスに 0 を指定してリストから要素をを取り出しています。ここで求めた new_Wb11 を使用して Wb の要素 Wb_11 を置き換えます。

replace_Wb11 = np.vectorize(lambda x: x.subs({Wb[0, 0]: new_Wb11}))

上記は、全要素の数式に含まれる Wb11 を new_Wb11 と置き換える関数を定義しています。

dVdWa_func = lambda Wa_val: np.array([w[0].subs(dict(zip(chain(Data.flatten(), Wa.flatten()),
                                       chain(data.flatten(), Wa_val.flatten()))))
                                      for w in dVdWa], dtype=np.float64).reshape((3, 1))

上記は Wa の勾配を求める関数を定義しています。シンボルに値を代入するために、key にシンボルを value に シンボルに対応する値を入れた辞書を渡しています。(TensorFlow がわかる人には placeholder に値を代入しているというと伝わるかもしれません。)また、前回記事では、微分をするのにjacobian を使用していましたが今回は for文で回して、各要素を微分しています。

dV2dWb_func = lambda Wb_val, Wa_val: np.array([w[0].subs(dict(zip(chain(Data.flatten(), Wb.flatten(), Wa.flatten()),
                                       chain(data.flatten(), Wb_val.flatten(), Wa_val.flatten()))))
                                               for w in replace_Wb11(dV2dWb)], dtype=np.float64).reshape((3, 1))

上記は Wb の勾配を求める関数を定義しています。 Wb を求めるときに Wa の情報が必要でコーディング時に注意と冒頭で説明しましたが、 上記がそのコードになります。 dict の中身をみると Wb だけでなく Wa の情報も使用していることがわかります。 また、 dV2dWb に直交条件を反映させるため replace_Wb11 関数を適用しています。これにより、 Wa に直交する Wb における分散の勾配を求めることができます。

7. 出力結果

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

===========================================
================ 第1主成分 ================
===========================================
================ 第1主成分 学習 i =  0 ================
分散V =  532.873871148815
Wa =  [[0.50212335 0.54326168 0.67285874]]
================ 第1主成分 学習 i =  10 ================
分散V =  533.003128892760
Wa =  [[0.48958803 0.55232128 0.67471829]]
================ 第1主成分 学習 i =  20 ================
分散V =  533.003128892760
Wa =  [[0.48958803 0.55232128 0.67471829]]
================ 第1主成分 学習 i =  30 ================
分散V =  533.003128892760
Wa =  [[0.48958803 0.55232128 0.67471829]]
================ 第1主成分 学習 i =  40 ================
分散V =  533.003128892760
Wa =  [[0.48958803 0.55232128 0.67471829]]
================ 第1主成分 学習 i =  50 ================
分散V =  533.003128892760
Wa =  [[0.48958803 0.55232128 0.67471829]]
================ 第1主成分 学習 i =  60 ================
分散V =  533.003128892760
Wa =  [[0.48958803 0.55232128 0.67471829]]
================ 第1主成分 学習 i =  70 ================
分散V =  533.003128892760
Wa =  [[0.48958803 0.55232128 0.67471829]]
================ 第1主成分 学習 i =  80 ================
分散V =  533.003128892760
Wa =  [[0.48958803 0.55232128 0.67471829]]
================ 第1主成分 学習 i =  90 ================
分散V =  533.003128892760
Wa =  [[0.48958803 0.55232128 0.67471829]]
================ 第1主成分 学習 i =  100 ================
分散V =  533.003128892760
Wa =  [[0.48958803 0.55232128 0.67471829]]
===========================================
================ 第2主成分 ================
===========================================
================ 第2主成分 学習 i =  0 ================
分散V2 =  0.282286132046230
Wb =  [[-0.85629139  0.15854369  0.49155768]]
================ 第2主成分 学習 i =  10 ================
分散V2 =  0.303318086521561
Wb =  [[-0.86855083  0.24061133  0.43327317]]
================ 第2主成分 学習 i =  20 ================
分散V2 =  0.315596549516175
Wb =  [[-0.87190187  0.30165001  0.38573876]]
================ 第2主成分 学習 i =  30 ================
分散V2 =  0.322411743461676
Wb =  [[-0.87100083  0.34594949  0.3488216 ]]
================ 第2主成分 学習 i =  40 ================
分散V2 =  0.326088827550761
Wb =  [[-0.86849373  0.37775612  0.32096564]]
================ 第2主成分 学習 i =  50 ================
分散V2 =  0.328042467143946
Wb =  [[-0.8656847   0.40051315  0.30029855]]
================ 第2主成分 学習 i =  60 ================
分散V2 =  0.329071954623274
Wb =  [[-0.86312883  0.41679385  0.28511667]]
================ 第2主成分 学習 i =  70 ================
分散V2 =  0.329612106418076
Wb =  [[-0.86100686  0.42845606  0.27403029]]
================ 第2主成分 学習 i =  80 ================
分散V2 =  0.329894869159100
Wb =  [[-0.85933015  0.43682358  0.26596403]]
================ 第2主成分 学習 i =  90 ================
分散V2 =  0.330042715591538
Wb =  [[-0.85804391  0.44283655  0.26010852]]
================ 第2主成分 学習 i =  100 ================
分散V2 =  0.330119970975593
Wb =  [[-0.85707561  0.44716326  0.25586408]]

================ 検証 ================
--- 独自手法 ---
第1主成分 =  [0.48958803 0.55232128 0.67471829]
第2主成分 =  [-0.85707561  0.44716326  0.25586408]
--- sklearn.decomposition.PCA ---
第1主成分 =  [0.48958803 0.55232128 0.67471829]
第2主成分 =  [ 0.85440069 -0.45835872 -0.24475855]
--- コサイン類似度 ---
 第1主成分 cosθ = 1.0 (0.0°)
 第2主成分 cosθ = -0.9998720868190456 (179.0835687095191°)
--- 第1主成分と第2主成分における分散の比率(sklearn.decomposition.PCA) ---
pca.explained_variance_ratio_[0] / pca.explained_variance_ratio_[1] 
 = 1614.1609966447247
--- 第1主成分と第2主成分における分散の比率(独自手法) ---
V / V2 
 =  1614.57402082519

8. 考察

出力結果を見ると、第1主成分は前回と同じで、問題なく動作していることが確認できます。第2主成分については i = 0 のとき分散が 0.28 で i = 100 のとき 0.33 と増加していることがわかりますので、学習はうまく進んだといえます。検証の箇所で sklearn.decomposition.PCA で求めた第2主成分とのコサイン類似度を出していますがおおよそ cosθ = -1 となりました。これは、 ベクトルが反対方向を向いているということで、 分散の計算でこのベクトルの符号は打ち消されますので、この方向でも反対方向でも分散が最大なのは変わりません。ちなみにどちらの方向になるかは Wb_val の初期状態できまります。最後に、 第2主成分の分散の値が 0.33 と少なく思えるため、分散の計算があっているか確かめるのに、 第1主成分の分散と第2主成分の分散の比率を sklearn.decomposition.PCA で導出したものと比較しています。これにより、分散の比率は 1614 とほとんど同じですので、分散の計算はあっていることがわかります。このサンプルデータにおいては第2主成分は分散が少ないですので、ほとんど意味がない情報といえそうです。以上にて、第2主成分について私の理解について問題ないことが確認できました。最後に、 3x3 のサンプルデータにおける最後の主成分の第3主成分ですが、これは第1主成分と第2主成分に直行するノルム1のベクトルですので連立方程式をとくと解が求まります。連立方程式を解いて第3主成分( Wc とします)を求めるコードを以下に記載し終わりとします。

9. 第3主成分を求めるコード

連立方程式を解いて第3主成分( Wc とします)を求めるコードを以下に記載し終わりとします。

# 連立方程式を解いて 最後の主成分である第3主成分を求める
Wc = np.array([[symbols('Wc%d1'%i)] for i in range(1, 4)])
wc_wa_inner = Wc.T.dot(Wa)[0, 0]
wc_wb_inner = Wc.T.dot(Wb)[0, 0]
wc_norm = Wc.T.dot(Wc)[0, 0]

eq1 = Eq(wc_wa_inner)
eq2 = Eq(wc_wb_inner)
eq3 = Eq(wc_norm - 1)

# eq1 かつ eq2 かつ eq3 の条件を満たす Wc の各要素の変数についてシンボリックに解く
kais = np.array(solve([eq1, eq2, eq3], *Wc.flatten()))

# 求めたシンボリックな解に Wa と Wb の値を代入して実際の値を求める関数を定義
subs = np.vectorize(lambda x: x.subs(dict(zip(chain(Wa.flatten(), Wb.flatten()),chain(Wa_val.flatten(), Wb_val.flatten())))))

# 解を表示
print('第3主成分 = ', subs(kais)[0])

# 連立方程式で求めた第3主成分と sklearn.decomposition.PCA で求めたものとの内積を求めて値が 1 に近いこと(方向が同じであること)を確認する
cos_theta3 = float(np.inner(subs(kais)[0], w3_pca))
print('--- コサイン類似度 ---')
print(' 第3主成分 cosθ = {} ({}°)'.format(cos_theta3, np.rad2deg(np.arccos(cos_theta3))))

10. 出力結果(第3主成分を求めるコード)

3主成分 =  [0.160390052085319 0.703552574067647 -0.692306873225203]
--- コサイン類似度 ---
 第3主成分 cosθ = 0.9998720868190454 (0.9164312904816843°)

第3主成分についても、コサイン類似度をみると、少し誤差はあるものの、sklearn.decomposition.PCA のものとほとんど同じであることが確認できました。誤差を少なくするには、勾配降下法で第1主成分と第2主成分を求めるときに学習率を大きくするなどの工夫が必要です。

11. まとめ

3x3 のサンプルデータについて第1主成分と第2主成分を数値解析的に求め、第3主成分を代数的に求めました。それぞれ sklearn.decomposition.PCA で求めた値ととほぼ同じ値になることを確認しました。PCAの理解を深めました。

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