2019-05-18

fashion-mnistにAutoencoderを適用してみる

目的

fashion-mnistにAutoencoderを適用し, 画像データがどのようにエンコードされているかを調べる.

今回は, fashion-mnistデータセットに対してAutoencoderを適用してみた. fashion-mnistはデータセット - Keras Documentationにも書かれているように, 28x28のサイズの服や靴といったファッション関連の画像を60000枚集めたデータセットである. 各データは10種類のラベルのうちのいずれかに分類される. それぞれのラベルを持つデータの代表例を挙げる.
f:id:gakushukun1:20190518144104p:plain

左上から

(0) Tシャツ/トップス
(1) ズボン
(2) プルオーバー
(3) ドレス
(4) コート
(5) サンダル
(6) シャツ
(7) スニーカー
(8) バッグ
(9) アンクルブーツ

に属するデータセットの一例である.

これらを含む600000枚のデータセットのうち, 半分の30000枚を使ってAutoencoderの学習を行う. 実際に学習に使ったネットワークの構造は次の通り.

エンコーダー

f:id:gakushukun1:20190518144615p:plain

デコーダー

f:id:gakushukun1:20190518144636p:plain

エンコーダー側は入力が28x28=784次元, 中間層が50次元, 出力が5次元であり, デコーダーはそれを反転させた構造を持つ. 学習の際には, エンコーダーとデコーダーを繋げて, もとの入力になるべく近い情報を出力するよう結合の重みを変えていく.

学習後に, 残りの30000枚のデータのうち, 10種類のラベルそれぞれに属しているデータをAutoencoderに通してみると, 次のように元の出力をある程度再現できていることが分かる. (各行1, 3番目の画像が元データ, 各行2, 4番目の画像が Autoencoderを通した画像)
f:id:gakushukun1:20190518171720p:plain

また, 30000枚のテストデータ全てについてエンコーダーに通し,その出力(5次元の情報)を各ラベルごとに平均をとった結果をデコーダーに通すと次のようになった.
f:id:gakushukun1:20190518175334p:plain

おおむね元のラベルが示す物の画像が出力されていることが確認できる. さらに, 前述したエンコーダーの出力を各ラベルごとに平均をとったものに対して, 5次元ある情報のうち1次元だけを0に変えてからデコーダーに渡す実験を行った. 得られた出力は次のようになる.
f:id:gakushukun1:20190518213916p:plain

縦の表示'input[i] == 0'は, デコーダーに渡したラベル0~9それぞれの平均ベクトル全てに対し, i番目の値だけ0に変えたことを意味する. この方法を用いることで, デコーダーの入力の各次元がどのような役割を担っているかを知ることができる.
i=0番目を0にする(図の1行目)と, 靴(サンダル, スニーカー, アンクルブーツ)かバッグを表していることが分かる. i=1番目を0にする(図の2行目)と, ズボン, プルオーバー, スニーカーを表していることが分かる. i=2番目を0にする(図の3行目)と, シャツ, プルオーバー, サンダルを表していることが分かる. i=3番目を0にする(図の4行目)と, シャツ, ズボン, プルオーバー, バッグを表していることが分かる. i=4番目を0にする(図の4行目)と, すべての図がぼやけ, どれともいえない画像となる.

Autoencoderは, このように情報を圧縮し, その組み合わせで画像を復元しているようだ.

まとめ

fashion-mnistを用いて, Autoencoderによる画像の符号化, 復号化の性質を調べた.

2019-05-16

ニューラルネットを用いた毒キノコの判定における入力次元削減

Keras

目的

キノコの外見的特徴のベクトルからニューラルネットを用いて毒キノコかの判定を行う. その際, 入力ベクトルのうちL1正則化を用いて不要な次元を削減する.

今回使用したデータセットはMachine Learning RepositoryのMushroom Data Setである. このデータセットはキノコの外見的特徴を表す22次元のベクトルと, そのキノコが毒キノコか食べられるキノコかを示す真のラベルのセットで構成されている.

この識別を行うにあたって構築したニューラルネットは次の図の通り.
f:id:gakushukun1:20190516235550p:plain

22次元の入力に対し, 中間層のユニット数を5とし, 最終的な出力は食べられるなら0, 毒キノコならば1となるようにした. なお各層の活性化関数はシグモイド関数である.
学習の際には, 全8124個あるデータのうち, 半分の4062個を訓練データとし, 残りの6500個のデータはテストデータとした. このデータセットは, いくつかのルールベースの方法によって簡単に分類できることが知られており(同梱の説明書参照), ニューラルネットで学習する場合はL1正則化のハイパーパラメータをかなり大きくしても, テストデータに対して100%近い識別率を出すことができる.
そこで、今回はハイパーパラメータ $\lambda$ を $0.01$ とかなり大きめに設定し, なるべく結合がスパースになるよう学習させた.
f:id:gakushukun1:20190516235133p:plain

100エポックだけ学習を繰り返した結果, 学習後のネットワークは上の図のようになった. ただし, 赤が正の結合, 青が負の結合を表しており, 結合荷重の絶対値が0.1より大きい結合のみ表示している. 最終的な訓練データに対する識別率は98.65%, テストデータに対する識別率は99.24%となった.
テストデータについての説明書によると, 上の図で入力19にあたるspore-print-color()がgreen=緑かどうかで, 99.41%の識別率を出すことができるとある. 今回のニューラルネットの学習結果にも入力19から出力の結合が残っているため, この部分の情報が推論に強い影響を与えている可能性がある.

まとめ

L1正則化を用いて, 入力次元の削減をする方法により, 毒キノコの判定を行った.

2019-05-13

車の評価の推測に用いた学習済みネットワークの構造の分析

Keras

目的

階層型ニューラルネットの学習結果を表示して, その構造の考察を行う。

gakushukun1.hatenablog.com

前回の記事で車の評価の推測を行った学習済みネットワークについて, 具体的な構造を図で確認し, 分析してみた.

前回の記事で示したネットワークの大まかな構造は次のようになっている.
f:id:gakushukun1:20190507232430p:plain

前回の実験では, このネットワークに対して, L1正則化のハイパーパラメータ $\lambda = 1.0 \times 10^{-5}$ を設定し, ユニット間の結合がスパースになるよう学習した.
今回は, その学習済みネットワークの各ユニット間の結合の重みを図示することで, 入力・出力間にどのような関係があるかを探ってみる.

具体的に学習済みネットワークを図示するための関数は次の通り. なお, 動作にはpydotおよびGraphvizが必要である.

def draw_graph(model, filename, weight=True, threshold=0, input_label=[], output_label=[]):

    input_shape = model.get_input_shape_at(0)[1]
    output_shape = 0

    graph = pydot.Dot(graph_type='digraph', rankdir='LR', splines='false', dpi='200', ranksep='1.5')
    graph.set_node_defaults(label='', shape='circle', width='0.3', height='0.3', fontsize='7', penwidth='0.3', margin='0')
    graph.set_edge_defaults(arrowsize='0.3', penwidth='0.3')
    idx = 0

    if len(input_label) == input_shape:
        input_label = []

        for i in range(input_shape):
            input_label.append(str(i))

    for i in range(input_shape):
        graph.add_node(pydot.Node('i.' + str(i), label=str(i)))
        graph.add_edge(pydot.Edge('i.' + str(i), str(idx) + '.' + str(i)))

    for layer in model.layers:
        if type(layer) != Dense:
            continue

        weights = layer.get_weights()[0]

        wx, wy = weights.shape

        max_w = np.max(weights.flatten())
        min_w = np.min(weights.flatten())

        for i in range(wx):
            for j in range(wy):
                edge = pydot.Edge(str(idx) + '.' + str(i), str(idx + 1) + '.' + str(j))
                if weight:
                    if weights[i][j] > threshold:
                        edge.set_penwidth(weights[i][j] / max_w)
                        edge.set_color('red')
                        graph.add_edge(edge)
                    elif weights[i][j] < -threshold:
                        edge.set_penwidth(weights[i][j] / min_w)
                        edge.set_color('blue')
                        graph.add_edge(edge)
                    else:
                        edge.set_style('invis')
                        graph.add_edge(edge)
                else:
                    graph.add_edge(edge)

        idx += 1
        output_shape = wy

    if len(output_label) == output_shape:
        output_label = []

        for i in range(output_shape):
            output_label.append(str(i))

    for i in range(output_shape):
        graph.add_node(pydot.Node('o.' + str(i), label=str(i)))
        graph.add_edge(pydot.Edge(str(idx) + '.' + str(i), 'o.' + str(i)))

    graph.write_png(filename)

この関数に対して, 学習済みのKerasのモデルを次のように渡す.

draw_nn.draw_graph(model, 'car_model.png', threshold=0.5)

ここで, thresholdは重みの絶対値がそれ以下の結合を図に表示させないための閾値である. 結果として得られる出力は次にようになる.
f:id:gakushukun1:20190513230023p:plain

与えられた閾値より, 図では絶対値0.5以下の重みの結合は省略されている. また, 重みが正ならば赤, 負ならば青で表示し, 重みの絶対値が大きい結合ほど太く表示するようにしている.
図中で, 0, 1, 2, 3, 4, 5で表されている入力は, それぞれ

0 ... buying(価格) 数値が大きいほど高い
1 ... maint(管理費) 数値が大きいほど高い
2 ... doors(ドアの数) 数値が大きいほど多い
3 ... persons(乗員数) 数値が大きいほど多い
4 ... lug_boot(載せられる荷物の量) 数値が大きいほど多い
5 ... safety(安全性) 数値が大きいほど高い

を表している.

一方, 0, 1, 2, 3で表されている出力は, 4つある評価ラベルそれぞれに所属する確率を表す. 評価ラベルの詳細は次の通り.

0 ... 良くない
1 ... まあまあ
2 ... 良い
3 ... とても良い

図を見ると, 入力3の「荷物の載せられる量」が, (正の結合-正の結合-正の結合)をたどるの2つのルートで出力3につながっていることが見て取れる. このことから, 「荷物の載せられる量」が車の高い評価に正の影響を与えていることが分かる. また, 入力0の「価格」が, (負の結合-正の結合-正の結合)をたどるルートで出力3に, (負の結合-負の結合-負の結合)をたどるルートで出力1,2につながっていることも見て取れる. このことからは, 「価格」の高さが上位３個の評価(まあまあ, 良い, とても良い)に負の影響を与えていることが分かる. さらに, 入力2の「ドアの数」はあまりどの出力にも影響を与えていないことも分かる.

まとめ

階層型ニューラルネットの学習結果を図示し, 結合の正負・重みを確認することで, 入力から出力へどのような推論が行われているかの分析を行った。

2019-05-08

ニューラルネットを用いた車の評価の推測

Keras

目的

階層型ニューラルネットをKerasで実装し, 実データに適用した学習の例を紹介する。

今回は, Kerasを用いて構築したニューラルネットのモデルを用いて, Car Evaluation Data Setから車の評価を推測してみる.

Car Evaluation Data Setには, それぞれの車が持つ6つの特徴と, その車に対する評価を示すラベルがセットとなったデータが1728個含まれている. このデータのうち半分を学習に用い, もう半分をテストデータとして検証に用いることとする.

各車が持つ6つの特徴は次の通り.

buying(価格) とても高い, 高い, 普通, 安いの四段階
maint(管理費) とても高い, 高い, 普通, 安いの四段階
doors(ドアの数) 2枚, 3枚, 4枚, 5枚以上の四段階
persons(乗員数) 2人, 4人, 5人以上の四段階
lug_boot(載せられる荷物の量) 少ない, 普通, 多いの三段階
safety(安全性) 高い, 普通, 低いの三段階

各車の評価を示すラベルは, unacceptable(良くない), acceptable(まあまあ), good(良い), very-good(とても良い)の四段階に分けられている.

さて, 今回はKerasを用いて次のようなネットワークを構築した.
f:id:gakushukun1:20190507232430p:plain
二つの中間層では, 活性化関数としてsigmoid関数を用いている. また, 過学習を抑えるために, 活性化の前にBatchNormalizaionを挟んだほか, L1正則化によるパラメータの罰則項を損失関数に加えている. 出力層では, softmax関数を用いて, 4個あるクラスそれぞれに所属する確率(合計で1になる)を出力するようにしている. モデルの学習には, Adamを用い, 交差エントロピーを損失関数とする.

6つの特徴を入力するにあたって, 各入力値が-1.0から1.0の範囲に含まれるように調整する. 例えば, 車の価格は[とても高い, 高い, 普通, 安い]の四段階で示されているが, [-1.0, -0.33, 0.33, 1.0]の四段階の値に対応づけることとする.

これらの条件を用いて, L1正則化のハイパーパラメータ $\lambda$ を変化させながら, 864個の学習データを2000epochミニバッチで学習させた. 2000epochループさせた後の, 残りのテストデータ864個についてのテスト誤差と正答率は次のように変化した.
f:id:gakushukun1:20190508155251p:plain

上のグラフが $\lambda$ ごとのテスト誤差の変化, 下のグラフが $\lambda$ ごとの正答率の変化である. この結果から, $\lambda = 1.0 \times 10^{-5}$ の時最も汎化性能が高くなったと判断できる.

$\lambda = 1.0 \times 10^{-5}$ とした時のepochごとの訓練誤差, テスト誤差, 訓練データとテストデータに対する正答率を示す.
f:id:gakushukun1:20190508161423p:plain

$\lambda$ を上記のように選んだところ, 学習誤差もテスト誤差もepochごとに減っていき, 過学習するような挙動は見られなかった.

今回の実験で用いたコードを以下に示す.

import tensorflow as tf
import numpy as np
import csv
from keras.models import Sequential
from keras.layers import Dense, Input, BatchNormalization, Activation
from keras.activations import sigmoid, softmax
from keras.regularizers import l1
from keras.optimizers import Adam
from keras.utils import np_utils, plot_model
import matplotlib.pyplot as plt


information = [
    ['vhigh', 'high', 'med', 'low'],
    ['vhigh', 'high', 'med', 'low'],
    ['2', '3', '4', '5more'],
    ['2', '4', 'more'],
    ['small', 'med', 'big'],
    ['low', 'med', 'high'],
]


def convert_to_value(str_vec):
    size = len(str_vec)
    conv_dict = {}

    for i in range(size):
        conv_dict[str_vec[i]] = -1.0 + (2.0 / (size - 1)) * i

    return conv_dict


def convert_to_data(row):
    data = [0] * (len(row) - 1)

    for i in range(len(row) - 1):
        data[i] = convert_to_value(information[i])[row[i]]

    label = {'unacc': 0 , 'acc': 1, 'good': 2, 'vgood': 3}[row[- 1]]
    label = np_utils.to_categorical(label, 4)

    return np.array(data), label


def read_file():
    f = open('car.data', 'r')
    reader = csv.reader(f)
    inputs = []
    labels = []

    for row in reader:
        data, label = convert_to_data(row)
        inputs.append(data)
        labels.append(label)

    f.close()

    train_size = int(len(inputs))

    rand_idx = np.random.permutation(range(len(inputs)))

    inputs = np.array(inputs)[rand_idx]
    labels = np.array(labels)[rand_idx]

    return inputs[0:train_size], labels[0:train_size]


def plot_loss_accuracy(fit, l1_lambda):
    fig, (loss_f, acc_f) = plt.subplots(2, 1, figsize=(15, 20))
    loss_f.plot(fit.history['loss'], label="train")
    loss_f.plot(fit.history['val_loss'], label="test")
    loss_f.set_xlabel("epoch")
    loss_f.set_ylabel("loss")
    loss_f.legend()
    acc_f.plot(fit.history['acc'], label="train")
    acc_f.plot(fit.history['val_acc'], label="test")
    acc_f.set_xlabel("epoch")
    acc_f.set_ylabel("accuracy")
    acc_f.legend()

    filename = 'car_fit_plot' + str(l1_lambda) + '.png'

    fig.savefig(filename)
    plt.close()


if __name__ == '__main__':

    train_inputs, train_labels = read_file()

    l1_lambda = 0.000001

    model = Sequential()
    model.add(Dense(10, activity_regularizer=l1(l1_lambda), input_shape=(6,)))
    model.add(BatchNormalization())
    model.add(Activation(sigmoid))
    model.add(Dense(5, activity_regularizer=l1(l1_lambda)))
    model.add(BatchNormalization())
    model.add(Activation(sigmoid))
    model.add(Dense(4, activity_regularizer=l1(l1_lambda)))
    model.add(BatchNormalization())
    model.add(Activation(softmax))

    model.compile(optimizer=Adam(), loss='categorical_crossentropy', metrics=['accuracy'])

    model.summary()
    plot_model(model, to_file='model.png', show_shapes=True)

    fit = model.fit(train_inputs, train_labels, epochs=2000, validation_split=0.5)

    plot_loss_accuracy(fit, l1_lambda)
<||

**まとめ

2019-05-06

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

PyStan

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

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のユーザーは, 飲食店と公園を低く評価する傾向にあり, 宗教施設は高く評価している.

まとめ

個人による旅行のレビューの実データをクラスタリングすることで, 飲食店に対する評価と宗教施設に関する評価は逆転する傾向があることがわかった。

2019-05-04

PyStanで混合ベルヌーイ分布を実装

PyStan

目的

0と1を要素とする長さNのベクトルが作るn個のデータがあるとき、その集合をクラスタリングするモデルとして混合ベルヌーイ分布がある。ここではそのベイズ推測をPythonによるStanのインターフェースであるPyStanで行う方法をまとめる。

混合ベルヌーイ分布とは, $\theta$ と $\mu$ をパラメータとする $x_i \in \{0, 1\}^N$ の分布であり, 次式で定義されるものである.

$\begin{eqnarray} \Pr(x_i | \boldsymbol{\theta}, \boldsymbol{\mu}) = \sum_{r=0}^{k - 1} \theta_r \prod_{j=0}^{N - 1} \mu_{rj}^{x_{ij}} ( 1 - \mu_{rj})^{1 - x_{ij}} \end{eqnarray}$

今回は, PyStanで混合ベルヌーイ分布のパラメータのサンプリングを行ってみた.

サンプリングを行うにあたって, パラメータ $\theta$ , $\mu$ の事前分布を

$\begin{eqnarray} \theta &\sim& {\rm Dirichlet}(\alpha, N) \\ \mu &\sim& {\rm Beta}(\beta_1, \beta_2) \end{eqnarray}$

とし, 入力するデータを, 既知の混合ベルヌーイ分布から生成された120個のサンプル $\{x_i\}$
(各 $x_i$ は $N=80$ 次元のベクトル)とする.

(ただし, $\theta = \{0.25, 0.25, 0.25, 0.25 \}$ ,
各 $r = 0, 1, 2, 3$ , $j = 0, 1, ..., 79$ について $20 * r \leq j < 20 * (r + 1)$ のとき $\mu_{rj} = 0.9$ , それ以外は $\mu_{rj} = 0.1$ とする)

また, サンプリングされた $\theta$ と $\mu$ を用いて, 各サンプル $x_i$ が $r$ に所属する事後確率 $y_{ir}$ を計算し, どのクラスタに所属しているかを推定する.

$y_{ir}$ は次のように表される.(ここで, $g_i$ は $x_i$ の所属するクラスタを表す )

$\begin{eqnarray} y_{ir} = \Pr(g_i = r | \{ x_i \}, \boldsymbol{\theta}, \boldsymbol{\mu}) = \frac{\theta_r \prod_j \mu_{rj}^{x_{ij}} ( 1 - \mu_{rj})^{1 - x_{ij}}}{\sum_s \theta_s \prod_j \mu_{sj}^{x_{ij}} ( 1 - \mu_{sj})^{1 - x_{ij}}} \end{eqnarray}$

データの生成と入力, 結果の出力を行うPythonのコードは次のように実装する.

import pystan
import numpy as np
import matplotlib.pyplot as plt
import collections


def make_cluster(result, n):
    cluster = collections.defaultdict(list)

    for i in range(n):
        cluster[np.argmax(result[i])].append(i)

    return cluster


def draw_matrix(X, n, cluster, k, title):
    plt.figure(figsize=(7,7))
    divide_pos = -0.5
    min_lim = -0.5
    max_lim = n - 0.5
    X_new = [X[i] for r in range(k) for i in cluster[r]]

    for r in range(k):
        if r == k - 1:
            break
        divide_pos += len(cluster[r])
        plt.hlines([divide_pos], min_lim, max_lim, linestyles="dashed")

    plt.imshow(X_new, cmap="Greys", interpolation="none")
    plt.tick_params(labelbottom='off')
    plt.tick_params(labelleft='off')
    plt.savefig(title)
    plt.show()


if __name__ == '__main__':
    n = 120  # サンプル数
    N = 80  # サンプル一つあたりの次元
    k = 4  # クラスタ数

    # テストデータのthetaの値
    real_theta = [1. / k for i in range(k)]

    # テストデータのmuの値
    real_mu = [[0.9 if (N / k) * i <= j < (N / k) * (i + 1)
                else 0.1 for j in range(N)] for i in range(k)]

    # テストデータの生成
    X = np.zeros((n, N))

    for i in range(n):
        c = int(np.random.choice(k, 1, p=real_theta))
        for j in range(N):
            mu = real_mu[c][j]
            X[i][j] = np.random.choice(2, 1, p=[1.0 - mu, mu])

    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, chains=4)

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

    # 集めたサンプルからクラスタ構造を推定
    cluster = make_cluster(samples['y'].mean(axis=0), n)

    # 結果の表示
    draw_matrix(X, n, cluster, k, 'y_sample')

Pythonのコードに読み込ませるStanのコードは次の通り.
事前分布のハイパーパラメータは次のように設定した.

$\begin{eqnarray} \alpha &=& \{ 1.0 \} \\ \beta_1 &=& 1.0 \\ \beta_2 &=& 3.0 \end{eqnarray}$

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(1.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]);
    }
}

推測されたパラメータは次のようになった.
f:id:gakushukun1:20190504151136p:plain
パラメータから推測されたクラスタ分割から, サンプル $\{x_i\}$ の各ベクトル(0=白, 1=黒)を縦に並べなおすと次のようになる. おおむねあらかじめ設定したサンプル生成元の分布と一致することがわかる.
f:id:gakushukun1:20190504155058p:plain

今回は成功した例だが, Stanでこのモデルのパラメータを推測すると失敗する場合もあった.

まとめ

混合ベルヌーイ分布を用いてクラスタリングを行うコードを紹介した. 実際に混合数kやハイパーパラメータを定めるには, 交差検証や情報量規準などが必要になる。

2019-04-29

PythonでBox-Cox変換を試してみる

前処理

今回の目的

正の実数に値をとるサンプルの分析のための前処理を行う方法として Box-Cox 変換を紹介する。

サンプル $\{y_i\}$ が正規分布でない分布から生成されている場合に, Box-Cox変換を用いることで, その分布を正規分布に近づけることができる.
Box-Cox変換は, 適当な $\lambda$ が与えられた上で, 次の式で表される:

$y_i^{(\lambda)} = \begin{cases} \frac{y_i^{\lambda} - 1}{\lambda} & \lambda \neq 0\\ \log y_i & \lambda = 0 \end{cases}$

scipy.stats.boxcoxを使うと, 適切な $\lambda$ を得た上で, 与えたサンプルをBox-Cox変換してくれる.
例として,

$\begin{eqnarray*} y \sim 2\exp(-2x) \\ y \sim {\rm Beta}(2, 5) \end{eqnarray*}$

の二つの分布から得られたサンプルをBox-Cox変換してみる.

import numpy as np
from scipy import stats
import matplotlib.pyplot as plt


if __name__ == '__main__':

    plt.close('all')
    np.random.seed(0)

    y1 = np.random.exponential(scale=0.5, size=1000000)
    y1_lambda, lmax1 = stats.boxcox(y1)
    y2 = np.random.beta(2.0, 5.0, size=1000000)
    y2_lambda, lmax2 = stats.boxcox(y2)

    bins = 100
    plt.rcParams["font.size"] = 15
    plt.figure(figsize=(10, 10))
    plt.subplot(2, 2, 1)
    plt.hist(y1, bins=bins)
    plt.title(r'$y_1$')
    plt.subplot(2, 2, 2)
    plt.hist(y1_lambda, bins=bins)
    plt.title(r'$y_1^{(\lambda)} (\lambda=%f)$' % lmax1)
    plt.subplot(2, 2, 3)
    plt.hist(y2, bins=bins)
    plt.title(r'$y_2$')
    plt.subplot(2, 2, 4)
    plt.hist(y2_lambda, bins=bins)
    plt.title(r'$y_2^{(\lambda)} (\lambda=%f)$' % lmax2)
    plt.tight_layout()
    plt.savefig('figure.png')
    plt.show()

結果:
f:id:gakushukun1:20190429125931p:plain

内部的には, 対数尤度関数

$llf(\lambda) = (\lambda - 1) \sum_i(\log(y_i)) - \frac{N}{2} \log(\frac{\sum_i (y_i ^{(\lambda)}- \bar{y^{(\lambda)}})^2 }{N})$

を最大化するような $\lambda$ を $-2.0 < \lambda < 2.0$ の範囲で見つけて返している.

まとめ

Box-Cox変換は最尤推定によってデータを正規分布からのサンプルに近づけるものであることが分かった. また, 実際にBox-Coxを使うことで, サンプルの分布を正規分布に近づけることができた.

Gakushukun1’s diary

20代エンジニア, 統計的機械学習勉強中　twitter: @a96665004

fashion-mnistにAutoencoderを適用してみる

目的

まとめ

ニューラルネットを用いた毒キノコの判定における入力次元削減

目的

まとめ

車の評価の推測に用いた学習済みネットワークの構造の分析

目的

まとめ

ニューラルネットを用いた車の評価の推測

目的

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

まとめ

PyStanで混合ベルヌーイ分布を実装

目的

まとめ

PythonでBox-Cox変換を試してみる

今回の目的

まとめ