イキサツ

さて，ノンパラメトリックベイズの勉強会も第2回目になりました(続けられてえらい!) ．今回は，前回の記事の続きとなる分割モデルについて解説を書いていこうと思います．前回の記事をお読みでない方は，下記から飛べますのでそちらから読んでいただければ幸いです．

また，引き続き参考書はオーム社の「続・わかりやすいパターン認識―教師なし学習入門」で，本記事もNotionAIとChatGPTがお供についてきてくれています．

分割モデル

分割モデルは読んで字のごとく，データを複数のグループに分割することで，より柔軟なモデルを構築する手法です．例えば，ある企業の顧客データを分析する場合，顧客を複数のグループに分割することで，グループごとに異なる傾向を捉えることが可能です．通常，分割モデルではベイズモデルを用いて，各グループのパラメータを推定します．一方，クラスタリングの数を事前に決定する必要があるため，データの性質によって正しいクラスタ数を事前に決めることが難しい場合があります．
分割モデルの代表的な手法には，本記事で紹介するホップの壺や中華料理店過程(CRP)の他にも，ディリクレ過程を用いた手法や階層モデルを用いた手法があります．ディリクレ過程を用いた手法では，クラスタリングの数を自動的に推定することができるため，データの性質に合ったクラスタリングを行うことができます(これを紹介してノンパラメトリックベイズは終了です)．また，階層モデルを用いた手法では，観測データの集合を複数の階層に分割することで，より柔軟なモデルを構築することができます．

ホップの壺

ホップの壺は分割モデルにおける階層事前分布の表現に関するもので、各グループの分布におけるパラメータを推定するために使用されます．具体的には，各グループにおけるパラメータの共通事前分布を仮定し，各グループにおけるパラメータの事後分布を推定します．以下で具体的な分割法の例を紹介します．

確率計算

まず，壺の中には1個の黒玉と他の色玉が入っており，黒玉の重さをα(>0) ，色玉の重さは全部1とします．この時，下記のルールに従って壺から玉を取り出し，入れる，取り出し，入れるを繰り返します．

1. 壺には最初，黒玉が1個のみが入っている．

2. 壺から玉の重さに比例した確率でランダムに玉を1つ取り出す．

3. 取り出した玉が黒なら壺内にない新しい色の玉を，取り出した玉が色玉なら同じ色の玉をツボに入れる．

4. 取り出した玉を戻し，壺を良く振る(シェイクシェイクブギーな胸騒ぎ)

5. 2.に戻り繰り返す．

さて，上記のルールを何回か繰り返した後のことを考えてみます．現在のツボの中には，n-1個の色玉が入っていて，第i色の色玉の数のniと置きます．
この時，次に入れるn番目の玉の確率はどうなるでしょうか．
答えは，第i色の色玉ならni / (n - 1 + α)で，新しい色(黒玉を引く)ならα / (n - 1 + α)になります．これも単純ですね．玉の引く確率は重みに比例するので，（対象の玉の重み / 全ての重み）で計算が可能です．
では，実際の例として「赤，青，青，赤，赤，黄」の順番で球が入れられたとすると，次のような確率計算ができます．

ここで，AFは上昇階乗と呼ばれ，AF(α,n)はα×(α+1)×…×(α+n-1)で表される．最後は(α+n)ではなく，(α+n-1)であることに注意!
さて，ホップの壺の計算に戻るが，ちょっとした解説を追加しよう．簡易的な状況説明である．最初，壺には黒玉しか入っていないため，これを引くしかなく，確率的にはα/αとなる．黒玉を引いたのでランダムに赤玉を壺の中に放り込む．次は黒玉or赤玉を引くがここでも黒玉を引いたため，同様にランダムに選ばれた青玉を壺の中に放り込む．この時の確率はα/(α+1)である．現在，壺の中には黒玉or赤玉or青玉の三種類だが，今回は青球を引いたため壺の中には青球を放り込む．この時の確率は1/(α+2) となる．これを繰り返し行う．

この分割モデルでは新しい色と既存の色という大別して2種類の玉しか考えない．このため，玉の色が「橙，緑，緑，橙，橙，紫」と変化しても確率は変化しない．また，この確率計算は動作を行う時点の壺の中身で判定しているため，順序も関係ない．このため，「黄，青，青，赤，赤，赤」でも確率は同じになる．実際に計算していたただければ納得できるかと(直観には反しますが)．
さて，上記の計算から一般的な式を開示します．先ほどの説明から，この確率は順序や玉の色によらず，色玉の種類数と個数のみに依存します．例えば，c種類の色玉が壺の中にn個入っていたとしたとき，このとき得られる色玉系列の生起確率PE(n1,…,nc)は

ここで，nは各niの総和になります．上式はイーウェンスの抽出公式になります．この式は統計学における分散分析において使用される公式の1つです．分散分析の手法を用いて，複数のグループ間での平均値の差異を検定する際に使用されます．

特徴

これまでの計算過程や手法からホップの壺の特徴をまとめると，

1. 色玉の組み合わせ確率は順番に依存しない交換可能性を持つ

2. 壺内の色玉の種類は試行が増えるたびに増加し，無限の種類の色玉が想定される

3. 黒玉の重さ(α)によって，色玉の種類が増減する

特に特徴1の交換可能性が大事になります．これがないとある手法を用いたサンプリングができません．入れ替えても確率が変わらないといったことは，中々に面白く重要です．

コーヒーブレイク

「ホップの壺について」とNotion AIに聞いてみたら，以下のような内容が返ってきました．

なるほど．壺と言ってますがおそらく窯のことかな？ちょっと，文章がおかしくなってますが，フィルターがついているといった知識を披露してますね．ちなみに，ビールは仕込、発酵、熟成、ろ過といった工程を経てつくられるらしいですね!

私はビールが大好きなので，記事を書きながらクスっと笑ってしまいました．

中華料理店過程

さて，続いてはノンパラメトリックベイズでは有名な中華料理店過程(CRP: Chinese Restaurant Process)について書いていきます．これは，先ほどのホップの壺の特徴を満たす分割モデルとなります．知らない人は「何やら美味しそうな名前」だと思われますが，実際の中華料理店を想像しながらお読みください．小説や参考書などは文字を追うのではなく，想像しながら自身の中で反芻することをお勧めします．なので，皆様も床や天井が赤く，丸テーブルにいくつもの美味しそうな料理が置かれている中華料理店を想像してみてください．

ほら，良いにおいがしてきたでしょ？

さて，このCRPでは中華料理店のテーブル対し，順番に客が着席していくことを考えていきます．この中華料理店では基本的に相席になります．さて，以下のルールに従って，客をテーブルに着かせましょう．お店は開店したばかりだと思って下さい．

1. 最初の客は任意のテーブルに着席する

2. 次以降の客(n番目)は空いているテーブルに$${\alpha / (\alpha+n-1)}$$の確率で着席，すでに$${n_i}$$人が着席している席に$${n_i / (\alpha+n-1)}$$の確率で着席する．

私だったら，迷わず空いている席ですが人間的にいかずにデータ的に分類していきましょう．さて，この状況の例を下図で示します．

テーブルには既に6人の客が座っている状況で，7人目の来店になります．各客の頭の数字は座った番号になってます．これを順次計算すると，最初の客はテーブル1に$${\alpha / \alpha}$$で着席，次の客は$${\alpha / (\alpha+1)}$$でテーブル2に，3番目の客は$${\1 / (\alpha+2)}$$で1番テーブルに，4番目の客は$${\alpha / (\alpha+3)}$$でテーブル3に…といった具合です．仮に7番目の客がテーブル4に座った場合，テーブルの着席順は(1,2,1,3,3,1,4)となり，この確率$${P(1,2,1,3,3,1,4)}$$は次式で計算される．

$${P(1,2,1,3,3,1,4) = \frac{\alpha}{\alpha}\times \frac{\alpha}{\alpha+1}\times \frac{1}{\alpha+2}\times \frac{\alpha}{\alpha+3}\times \frac{1}{\alpha+4}\times \frac{2}{\alpha+5}\times \frac{\alpha}{\alpha+6}}$$

この確率計算はテーブルの番号や客の順番が変わっても同様の計算結果が示されます．これが，先ほどのホップの壺で出てきた交換可能性になります．整理すると下記の特徴がCRPから得られます．

1. 客が着席する確率は順番に依存しない

2. 着席しているテーブル数は客とともに増加する

3. $${\alpha}$$が大きいほど空席じゃないテーブルが増加し，小さい場合は相席のテーブルにさらに客が増える

こられの特徴量は先ほども見たことがあるものだと思われる．ホップの壺で出てきたものとそっくりですね．

具体例

参考書でも計算している通り，客を4人$${(n=4)}$$としたときの計算をしていきましょう．ここで，4名の客をテーブルに分割するパターンというのは以下のものが考えられます．

ここで，少し厄介なのですが客の区別はつけるのに対し，テーブルの区別はつけません．なので①②，③④と①③，②④の区別つけますが③④，①②の区別はつけません．ちなみに着席確率はホップの壺の後半で出てきたイーウェンスの抽出公式を用いて計算できます．例えば，①②，③④になる確率は次式の様に計算されます($${\alpha＝2}$$としてます)．

$${P_E(2,2) = \frac{2^2\prod_{i=1}^{2}(n_i-1)!}{AF(2,4)} = \frac{2^2 \times 1! \times 1!}{2\times3\times4\times5}= \frac{1}{30} }$$

ちなみに，この確率は客ごとの分割を加味した確率となっている．これをクラスタ毎，つまり4人の客が2,2に分かれる確率は

$${P^*(2,2) = N_1 \times P_E(2,2) = 3\times \frac{1}{30} =  \frac{1}{10}}$$

さっきまではテーブルを区別しませんでしたが，区別する場合を考えます．この場合は$${c!}$$通りのパターンが考えられるので生起確率は

$${\frac{1}{2!} \times \frac{1}{30} =  \frac{1}{60}}$$

となります．このため，事前確率$${P(s_1 \cdots s_n)}$$は

$${P(s_1 \cdots s_n) = \frac{1}{c!}P_E(n_1 \cdots n_c)}$$

として求めることが可能です．クラスタリングの場合では，各クラスタにおパラメータ$${\bm \theta}$$が対応するため，上式の計算のようなことが行われます．

プログラムと実験

実際にプログラムをしていきましょう．今回作成したプログラムが以下の通りです．

import numpy as np
import matplotlib.pyplot as plt

# 中華料理店過程の確率的分割モデル
class ChineseRestaurantProcess(object):
    def __init__(self, alpha):
        self.alpha = alpha
        self.tables = []
        self.customers = []
        self.table_assignments = []
        
    # 顧客が席に座る
    def seat_customer(self):
        if not self.tables:
            # 1番目のテーブルに座る
            table = 1
            self.tables.append([0])
            self.table_assignments.append(0)
        else:
            # 確率分布から、新しいテーブルに座るかどうかを決定
            table_probs = self.compute_table_probs()
            new_table_prob = self.alpha / (len(self.customers) + self.alpha)
            table_probs = np.append(table_probs, new_table_prob)
            table_probs /= sum(table_probs)
            # テーブルに座る
            table = np.random.choice(len(table_probs), p=table_probs)
            if table == len(self.tables):
                self.tables.append([len(self.tables)])
            self.table_assignments.append(table)
            self.tables[table].append(len(self.customers))
        # 新しい顧客を追加
        self.customers.append(len(self.customers))
        
        
    # 各テーブルの人数に基づいた確率分布を計算
    def compute_table_probs(self):
        n_tables = len(self.tables)
        table_probs = np.zeros(n_tables)
        for i in range(n_tables):
            table_probs[i] = len(self.tables[i]) / (self.alpha + len(self.customers) - 1)
        return table_probs
 
    
 # パラメータの設定
alpha = 10
n_customers = 100
CRP = ChineseRestaurantProcess(alpha)

# サブプロットの設定
fig, axs = plt.subplots(3, 3, figsize=(10, 10))


# 複数回実行し、グラフを描画
for i in range(9):
    # 中華料理店過程の初期化
    CRP = ChineseRestaurantProcess(alpha)

    # 顧客をシミュレート
    for j in range(n_customers):
        CRP.seat_customer()
    
    # テーブルの割り当て数をカウント
    table_counts = {}
    for table in CRP.table_assignments:
        if table in table_counts:
            table_counts[table] += 1
        else:
            table_counts[table] = 1
    
    # グラフの描画
    ax = axs[i // 3][i % 3]
    ax.bar(range(len(table_counts)), list(table_counts.values()), align='center')
    ax.set_xticks(range(len(table_counts)))
    ax.set_xticklabels(list(table_counts.keys()))
    ax.set_xlabel('Table')
    ax.set_ylabel('Number of customers')
    ax.set_title('Iteration {}'.format(i+1))

plt.tight_layout()
plt.show()

今回は複数回実行したときの結果を示します．まずは$${\alpha=2，n=100}$$の時の結果になります．

最大でテーブルが8個になってますね．傾向は思ったよりバラバラなのが意外ですね．続いて，$${\alpha=2，n=1000}$$の時の結果になります．理論的にはテーブル数$${c}$$が増えるはずです．

そこまで，大々的に増えているわけではないですが，テーブル数の増加が全体的に確認できます．素晴らしい！最後に$${\alpha}$$を変化させた場合を試してみましょう．$${\alpha}$$を増加させればテーブル数がさらに増えるはずですね．

めちゃめちゃ増えましたね！やはり客の数よりも，新規テーブルの獲得しやすさを増加させるほうが，テーブル数は増えますね！

まとめ

今回は前回に引き続き，自分の勉強ついでにノンパラメトリックベイズの分割もモデルについて書かせていただきました．次回はちょっと数式チックなディリクレ過程について書くと思います．
また，お恥ずかしいことに本記事の途中までnoteで数式モードが使えることを知らずに書いておりました．チュートリアルや使い方を学ぶということの大事さを身に染みて実感しております．「こんな便利機能もあるよ！」という方はコメントの方に書いていただければ幸いです．
加えて，皆様のおかげで累計のスキが80，全体ビューが500を超えそうです(2023年3月19日9時の時点)．フォロワーの方も着々増えており，大変感謝しております．これに関する記事も後ほど書こうかと思っておりますので，今後ともよろしくお願いいたします．

おまけ

最近GPT-4の公開がなされて，以前より精度が比較にならないほどいいと聞きます．一方，私は課金していないので以前の通りのChatGPTに「中華料理店過程のプログラムを組んで」とお願いしたら以下の回答が来ました．
ChatGPT「はい、中華料理店の過程のプログラムを組むことができます。以下は、中華料理店の過程の簡単なプログラム例です。」

# 中華料理店の過程のプログラム

# 材料のリスト
ingredients = ['鶏肉', 'エビ', 'キャベツ', '人参', '玉ねぎ', 'ニンニク', '生姜', '豆腐', '麺', 'ごま油', '醤油', '砂糖', '酢', '鶏がらスープ', '塩', '胡椒']

# 料理のメニュー
menu = {
    '炒飯': ['鶏肉', 'エビ', '玉ねぎ', 'ニンニク', '卵', '醤油', '塩', '胡椒', 'ご飯'],
    '餃子': ['豚肉', 'キャベツ', 'ニンニク', '生姜', '酢', '醤油', '塩', '胡椒', '皮'],
    'ラーメン': ['鶏がらスープ', '豚肉', 'キャベツ', '人参', '生姜', '麺', 'ニンニク', 'ごま油', '醤油', '塩', '胡椒']
}

# 料理の作り方を定義
def make_dish(menu_item):
    ingredients_needed = menu[menu_item]
    for ingredient in ingredients_needed:
        if ingredient not in ingredients:
            print(f"{ingredient}がありません。{menu_item}を作ることができません。")
            return
    print(f"{menu_item}を作ります！")
    for ingredient in ingredients_needed:
        ingredients.remove(ingredient)
    print("完成！")

# 料理を作る
make_dish('炒飯')
make_dish('餃子')
make_dish('ラーメン')

ChatGPT「このプログラムでは、まず材料のリストと料理のメニューを定義しています。次に、料理を作るための make_dish 関数を定義しています。この関数は、メニューから必要な材料を抽出し、材料が十分にある場合は料理を作り、材料が不足している場合は料理を作ることができないというメッセージを出力します。」

ChatGPT「最後に、make_dish 関数を使用して、炒飯、餃子、ラーメンの料理を作成します。このプログラムは、基本的な中華料理店の過程を再現していますが、実際の中華料理店で使用されるプログラムには多く」

私はちょっとおバカな子の方が好きなのかもしれません笑