ポアソン混合モデル　その6：ポアソン混合モデルによる学習、EMアルゴリズムの適用<実装編-2：一般のデータで試してみる>

クラスの紹介

下記のMixed_Poissonクラスを用いて、実データのクラスタリングを実行してみようと思う。

今回は、ヒストグラムの色を連続的に変化させるよう、描画用メソッドを追加した。

クラスのコードを以下に示す。

# クラス内でimportするのもあれだよね？
import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from matplotlib.pyplot import cm, rc
from matplotlib import colors
from sklearn.decomposition import PCA
from sklearn.preprocessing import MinMaxScaler


# Google Colabではこの分が無いとアニメーションが再生されない。
rc('animation', html='jshtml')

class Mixed_Poisson:
  def __init__(self, pi_vec, lam_vec, X=None, criteria=1e-5):
    if len(pi_vec) != len(lam_vec):
        raise ValueError(
            "shapes of pi_vec and lam_vec didn't match."
        )
    self.K = len(pi_vec)
    self.pi_history = np.array(pi_vec)[np.newaxis]
    self.lam_history = np.array(lam_vec)[np.newaxis]
    self.Q_history = np.array([0])[np.newaxis]
    self.gamma_matrix = None
    self.X = None
    self.color_map = None

  def execute_EM(self, X, n_iter=1):
    self.X = X
    for i in range(n_iter):
      self.execute_E_Step()
      self.execute_M_Step()
      self.check_convergence()

  def execute_E_Step(self):
    self.gamma_matrix = self.gamma_all()

  def execute_M_Step(self):
    self.lam_history = np.vstack((self.lam_history, self.update_lambda()))
    self.pi_history = np.vstack((self.pi_history, self.update_pi()))

  def check_convergence(self):
    self.Q_history = np.vstack((self.Q_history, self.calc_Q()))

  # 単一データ x_n を入力にとる
  def update_gamma(self, x):
    # k番目のクラスタでの確率質量とクラスタ重みの積。lambda関数として型だけ作る。
    partial_poisson = lambda k, x : self.pi_history[-1][k] * stats.poisson.pmf(x, mu=self.lam_history[-1][k])
    # k=1,2,3 について足し合わせる
    all_poisson = np.array([partial_poisson(k=k, x=x) for k in range(self.K)]).sum()
    # 各クラスタでのγの値を要素にもつnp.arrayをつくる
    gamma_n = [partial_poisson(k=k, x=x) / all_poisson for k in range(self.K)]
    return gamma_n

  # 複数データ X = (x_1, .., x_n) を入力にとる
  def gamma_all(self):
    return np.array([self.update_gamma(x_n) for x_n in self.X])

  def update_lambda(self):
    weighted_x = lambda k : np.sum(self.gamma_matrix[:, k] * self.X)
    weight_sum = lambda k : np.sum(self.gamma_matrix[:, k])
    lam_vec = [weighted_x(k=k) / weight_sum(k=k) for k in range(self.K)]
    return lam_vec

  def update_pi(self):
    pi_k = lambda k : np.sum(self.gamma_matrix[:, k]) / len(self.gamma_matrix)
    pi_vec = [pi_k(k=k) for k in range(self.K)]
    return pi_vec

  def calc_Q(self):
    Q = np.array([self.calc_q_n(n, x_n) for n, x_n in enumerate(self.X)]).sum()
    return Q

  # 単一データ x_n に対するQ関数の値
  def calc_q_n(self, n, x):
    # avoid ZeroDivisionError
    partial_log_poisson = lambda k : np.log1p(self.pi_history[-1][k] * stats.poisson.pmf(x, mu=self.lam_history[-1][k]))
    gamma = self.gamma_matrix[n]
    q_n = np.array([partial_log_poisson(k=k) * gamma[k] for k in range(self.K)]).sum()
    return q_n

  # 描画用メソッド
  # nth番目の推論結果を描画する。デフォルトでは最新の結果(ヒストリーのインデックス'-1'の値)を表示させる
  def show_one_inference(self, nth=-1):
    x_min = min(self.X)
    x_max = max(self.X)
    x = np.arange(start=x_min, stop=x_max)
    kth_poisson = lambda k : self.pi_history[nth][k] * stats.poisson.pmf(x, mu = self.lam_history[nth][k])

    mixed_poisson = np.zeros(shape=len(x))

    for k in range(self.K):
      poisson = kth_poisson(k=k)
      plt.plot(x, poisson)
      mixed_poisson += poisson

    plt.plot(x, mixed_poisson, color='orange')
    plt.hist(x=self.X, density=True, bins=100, color='silver')
  
  def show_inference_process(self):
    # colorはK次元ベクトル
    color = cm.brg_r(np.linspace(0, 1, self.K))

    x_min = min(self.X)
    x_max = max(self.X)
    x = np.arange(start=x_min, stop=x_max)

    fig = plt.figure()
    ims = []  

    for n_iter in range(len(self.pi_history)):
      kth_poisson = lambda k : self.pi_history[n_iter][k] * stats.poisson.pmf(x, mu = self.lam_history[n_iter][k])
      mixed_poisson = np.zeros(shape=len(x))

      nth_imgs = []

      for k, c in zip(range(self.K), color):
        poisson = kth_poisson(k=k)
        im, = plt.plot(x, poisson, color=c)
        nth_imgs.append(im)
        mixed_poisson += poisson

      im1, = plt.plot(x, mixed_poisson, color='black')
      nth_imgs.append(im1)

      ims.append(nth_imgs)
    
    N, bins, patches = plt.hist(x=self.X, density=True, bins=100, )
    self.gradate_histo(N, bins, patches)

    # 10枚のプロットを 100ms ごとに表示
    ani = animation.ArtistAnimation(fig, ims, interval=100)
    plt.close()
    plt.plot(np.arange(len(self.Q_history)), self.Q_history)
    # Google Colaboratoryを使用していて、gifとして保存する場合、以下のコマンド１．を別セルで実行後、２．をコメントアウトする。
    # １．!apt-get update && apt-get install imagemagick
    # ２．ani.save("output.gif", writer="imagemagick")
    # 参照：https://qiita.com/shinmura0/items/ed96863281637e4fa10c
    return ani

  # 取得したN,bins,patchesを用いてヒストグラムの色を変更する(戻り値なし)
  def gradate_histo(self, N, bins, patches):
    # k番目のbin以上でk+1番目のbinより小さいか」で真偽値を得る。唯一のTrueがあるインデックス番号がbinの番号に対応している O(NK)
    class_num = len(bins) - 1

    # [bins[k] <= x < bins[k+1]で作った場合、最大値がbins[k+1]と一致するので、戻り値が空のnp.arrayになってしまう。
    # np.where(条件)で、条件を満たすインデックス番号のnp.arrayを格納したサイズ1のタプルが帰る
    # np.whereの戻り値に[0][-1]とすることで、タプルを外し、最後のインデックスを取得する。
    class_indices = np.array([np.where([bins[k] <= x for k in range(class_num)])[0][-1] for x in self.X])

    # binのインデックスをkeyにとり、そのbinに含まれるデータのデータ番号(x_nのn)をvalueとする辞書を作成
    class_indices_dict = {i:np.where(class_indices==i)[0] for i in range(class_num)}

    # 一応ゼロ割回避
    bin_sums = [np.sum(self.gamma_matrix[class_indices_dict[i]], axis=0) / len(class_indices_dict[i]) if len(class_indices_dict[i]) != 0 else np.zeros(self.K) for i in range(class_num)]

    # 今回はPCAでπを1次元に次元削減する
    pca = PCA(n_components=1)
    scaler = MinMaxScaler(feature_range=(0,1))
    bin_color_1d = pca.fit_transform(X=bin_sums).flatten()
    norm = colors.Normalize(bin_color_1d.min(), bin_color_1d.max())

    for bin, patch in zip(bin_color_1d, patches):
      color = plt.cm.brg_r(norm(bin))
      patch.set_facecolor(color)

EMアルゴリズムのあたりのコードについては、別記事(ポアソン混合モデル　その5：ポアソン混合モデルによる学習、EMアルゴリズムの適用<実装編-1：やっと推論>)のコードにて多く注釈をつけているので、参照されたい。

gradate_histo メソッドでは、$ K $ 次元ベクトル $ \textbf{π} $ の情報からヒストグラムに色を付けている。

ヒストグラムのbinごとに色を付けることはあまり難しくないが、$ K $ 次元ベクトル $ \textbf{π} $ からスカラーに変換する方法がなかなか難しい。

今回は、単純にPCAによる次元削減を行った。

データの紹介

さて、今回使うデータの紹介だ。

data.worldというサイトで入手したBicycle Counts for East River Bridges- Historical (https://data.world/city-of-ny/gua4-p9wg) というデータをもとに、ポアソン混合モデルを用いたクラスタリングを実行してみようと思う。

このデータは、複数個所でカウントされた自転車の台数を日ごとに集計したデータである。

今回は簡単のため、Brooklyn Bridgeにおける自転車の台数についてクラスタリングを行ってみよう。

ということで、データの準備。Excel形式のファイル。

import pandas as pd
import glob

dfs = pd.DataFrame()

# Excelファイルを格納したディレクトリ
Path = '/content/drive/MyDrive/Colab Notebooks/tmp'

for excel in glob.glob(Path + '/*'):
  df = pd.read_excel(excel, header=5, usecols="C:K").dropna()
  dfs = pd.concat([dfs, df])

dfs.reset_index(inplace=True)

import datetime

dfs['Day of Week'] = dfs['Day'].apply(lambda x: x.strftime('%A'))

df1 = dfs.copy()[['Brooklyn Bridge', 'Day of Week']]

実は、このデータにMixed_Poissonクラスを用いてEMアルゴリズムを実行したら、エラーが生じた。

推測される原因は後ほど述べるが、これを回避するため、今回はデータの数値を1/60倍にスケーリングする。

自分でスケールを変えてどのようにグラフに変化があるか、見てみると面白い。

オリジナルスケールのデータのヒストグラムは以下である。

scale = 1/60

df1['Brooklyn Bridge Scaled'] = df1['Brooklyn Bridge'] * scale

plt.hist(df1['Brooklyn Bridge'], bins=100)

実験してみよう

クラスタ数 $ K = 2,3,4 $ について実験してみようと思う。上ヒストグラムを見て、$ K = 2,3,4 $ それぞれについて $ λ_k $ のそれっぽい初期値を作成した。

original_lam_vec_2 = np.array([2000, 3000])
original_lam_vec_3 = np.array([1500, 2500, 3000])
original_lam_vec_4 = np.array([1000,2000,3000,4000])

# #cluster = 2
pmm = Mixed_Poisson(pi_vec=[0.5, 0.5], 
                     lam_vec=original_lam_vec_2 * scale)
pmm.execute_EM(X=np.array(df1['Brooklyn Bridge Scaled'].astype(int)), n_iter=50)
pmm.show_inference_process()

# #cluster = 3
pmm = Mixed_Poisson(pi_vec=[0.3, 0.3, 0.4], 
                     lam_vec=original_lam_vec_3 * scale)
pmm.execute_EM(X=np.array(df1['Brooklyn Bridge Scaled'].astype(int)), n_iter=50)
pmm.show_inference_process()

$ K = 3 $ の場合、赤色のヒストグラム領域が不自然に少なくなっている。

$ \textbf{π} $ ベクトルに対するPCAによって不適切に情報が落ちたのかもしれない。

# #cluster = 4
pmm = Mixed_Poisson(pi_vec=[0.2, 0.2, 0.4,0.2], 
                     lam_vec=original_lam_vec_4 * scale)
pmm.execute_EM(X=np.array(df1['Brooklyn Bridge Scaled'].astype(int)), n_iter=50)
pmm1.show_inference_process()

$ K = 4 $ の場合、一つのポアソン分布がつぶれて消えていく様子が見て取れる。初期値の問題か、あるいは $ K ≧ 4 $ の仮定が不適切だったのか。

スケーリングの問題

もともと本データのスケールは、1000-3000であった。

このままEMアルゴリズムを実行すると、RuntimeWarningが出た。

実験として、データおよび $ λ_k (k = 1, 2, .., K) $ をスケーリングしてみた。

スケーリング1/10では、ポアソン分布は尖っているように見える。スケーリング1/60では、ポアソン分布はより裾を引いているように見える。だが、、、

例えばの話をしてみよう。

データおよび $ λ_k (k = 1, 2, .., K) $のスケールを1/10にしたら、データから計算される分散は $ V(1/10X)=1/100V(X) $ となる。

一方、ポアソン分布では、パラメータ $ λ $ は分布平均及び分散に等しいため、分散は $ 1/10 $ となる。

ということで、データはスケールを小さくしたらその分分散が小さくなるが、ポアソン分布はそれよりは分散が小さくならない。

1000-3000程度のスケールをもつ元データでは、ほとんど $ 0 $ の値となるような裾の先で値を計算していたため、計算が不安定になってRuntimeWarningが出たものと考えている。

正直あまり自信はないので、ご指摘を待っていたい。

ガウス分布を使っていれば分散も制御できたわけだし、なるほどポアソン分布はパラメータ数も少ないし学習向け教材としては優れているがガウス分布を使いたくなってきた。