混合ベルヌーイ分布でTravel Reviews Data Setを分類

**目的
混合ベルヌーイ分布によるクラスタリングを実データに適用して、どんな情報が発見できるかを試してみる。

gakushukun1.hatenablog.com

前回のエントリで実装した混合ベルヌーイ分布を用いて, Machine Learning RepositoryのTravel Reviews Data Setの分類を行ってみた.

Travel Reviews Data Setは, TripAdvisor.comから収集された東アジアの観光名所のレビューデータから, 各名所を10個のカテゴリに分類し, 980人のユーザーがそれぞれのカテゴリの名所に平均何点つけたかを表すデータである. 例として,データの一行目を見てみると,

	Category 1	Category 2	Category 3	Category 4	Category 5	Category 6	Category 7	Category 8	Category 9	Category 10
User 1	0.93	1.8	2.29	0.62	0.8	2.42	3.19	2.79	1.82	2.42

User1はカテゴリ1の名所には平均0.93点, カテゴリ2の名所には平均1.8点, カテゴリ3の名所には平均2.29点…をつけているといった情報が読み取れる.
なお, 各名所には0〜4点の範囲の点数がつけられており, 0が最低, 4が最高である. また, 名所のカテゴリはそれぞれ,

category1 … 美術館(art galleries)
category2 … ダンスクラブ(dance clubs)
category3 … ジュースバー(juice bars)
category4 … レストラン(restaurants)
category5 … 博物館(museums)
category6 … リゾート地(resorts)
category7 … 公園(parks/picnic spots)
category8 … ビーチ(beaches)
category9 … 劇場・映画館(theaters)
category10 … 宗教施設(religious institutions)

を表している.
さて, このデータを混合ベルヌーイ分布を用いて分類するにあたって, データの数値を0と1で表されるように前処理する必要がある. 今回は各カテゴリについて980個あるデータの中央値を調べ, 中央値より大きいデータは1, そうでなければ0とした.
実装したPythonのコードは次の通り. 分類するクラス数は6とした.

import pystan
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

# 所属クラスごとのmuの値をプロット
def plot_result(mu, category, k):
    plt.figure(figsize=(10,15))
    for i in range(k):
        class_label = 'class' + str(i)
        plt.plot(category, mu[i], label=class_label)
    plt.legend
    plt.show()


if __name__ == '__main__':
    category = ['art galleries', 'dance clubs', 'juice bars',
                'restaurants', ' museums', 'resorts', 'parks/picnic spots',
                'beaches', ' theaters', 'religious institutions']

    df = pd.read_csv("tripadvisor_review.csv", sep=",")
    data = df.values

    n = data.shape[0]  # サンプル数
    N = data.shape[1]  # サンプル一つあたりの次元
    k = 6  # クラス数
    print(data.shape)
    X = np.zeros((n, N))

    
    for i in range(n):
        for j in range(N):
            if data[i, j] >= np.median(data[:, j]):
                X[i, j] = 1

    X = X.astype(int)

    # Stanのコードに渡すデータ
    stan_data = {'n': n, 'N': N, 'k': k, 'X': X}

    # Stanのコードで定義されたモデルを生成
    sm = pystan.StanModel(file='mixture_bernoulli.stan')

    # サンプリング
    fit = sm.sampling(data=stan_data)

    fit.plot()
    samples = fit.extract()

    d = pd.DataFrame(data=samples['mu'][-1], index=np.arange(1, k+1), columns=category)

    print(samples['theta'][-1])
    plot_result(samples['mu'][-1], category, k)

    d.to_csv('trip_result_mu.csv', float_format='%.2f')

Stanのコードは次の通り. おおむね前回と同様だが, 混合比の事前分布のハイパーパラメータを $\alpha = \{5.0\}$ とした. この値が大きいと, 混合比が0のクラスができにくくなる.

data {
    int<lower=0> n;
    int<lower=0> N;
    int<lower=0> k;
    int<lower=0, upper=1> X[n, N];
}

parameters {
    simplex[k] theta;
    real<lower=0, upper=1> mu[k, N];
}

transformed parameters {
    vector[k] ps[n];

    for(i in 1:n){
        for(r in 1:k){
            ps[i, r] = log(theta[r]);
            for(j in 1:N){
                ps[i, r] += bernoulli_lpmf(X[i, j] | mu[r, j]);
            }
        }
    }
}

model {
    theta ~ dirichlet(rep_vector(5.0, k));
    for(r in 1:k){
        for(j in 1:N){
            mu[r, j] ~ beta(1.0, 3.0);
        }
    }
    for(i in 1:n){
        target += log_sum_exp(ps[i]);
    }

}

generated quantities {
    simplex[k] y[n];

    for(i in 1:n){
        y[i] = softmax(ps[i]);
    }
}

各クラスに所属するユーザーが, それぞれのカテゴリの観光名所にどれだけ中央値より高い点数をつける傾向があるかを示す $\mu$ をプロットすると次のようになった.
各折れ線が一つのクラスに対応している.
f:id:gakushukun1:20190506000047p:plain

また, 各ユーザーがそれぞれのクラスに所属する割合 $\theta$ は次のようになった.

$\theta = \{0.11, 0.14, 0.14, 0.40, 0.09, 0.12 \}$

$\theta$ の値をみると, 混合比の最も大きいのはクラス3であり, 図では赤の折れ線に対応している. このクラスに所属しているユーザーは, 飲食店と公園は高く評価する傾向にあるが, 宗教施設は低く評価している. 逆に２番目に大きいクラス1のユーザーは, 飲食店と公園を低く評価する傾向にあり, 宗教施設は高く評価している.