Python3とPyMC3でベイズ推定でPyMC3の環境を作成したので、もう少しベイズ推定をやってみる。参考にしたのはBayesian Methods for HackersというGitHubにあるPyMCによる確率的プログラミングのオンラインドキュメント。PyMC3用のChapter 1: Introduction to Bayesian Methodsでは、日ごとの受信メッセージ数の推移から転換点をベイズ推移する方法が解説されている。これを参考に、去年の流行語大賞である「インスタ映え」の転換点を推定してみる。
Bash on Ubuntu on Windows
Jupyter Notebook
Google TrendsでGoogleで検索されたワードの週ごとの検索数（最大100に正規化された数）を調べることができる。これを利用して「インスタ映え」の検索数を取得する。期間は2017年の1年間。
結果をcsvで保存する。
csvの中身は以下の通り。Google Trendsの出力は日付と1週間ごとの検索数（最大100に正規化されている）。
Google Trendsで取得したデータをグラフ化してデータを確認する。
グラフは以下のようになる。ベイズ推定しなくても転換点がだいたいどの辺りかわかりそうだけど、練習なのでベイズ推定してみる。
Chapter 1: Introduction to Bayesian MethodsのExample: Inferring behaviour from text-message dataでは、1日ごとのメッセージ受信数はポアソン分布に従い、転換点tauの前後でメッセージ受信数の平均λが変化するとしてモデル化している。同様に1週間ごとの検索数がポアソン分布に従うとして同じモデルを使ってベイズ推定する。
結果は以下のようになった。tau、lambda_1、lambda_2のRハットはすべて1.1より小さいので収束している。tauの事後分布の平均は28なので29週目が転換点と推定された。転換点前後の平均検索数は7.8と66.7回で、29週目（7月）の辺りで検索数が増えていると言える。この辺りから「インスタ映え」がメディアに取り上げられるようになったのだろうか。
環境
Bash on Ubuntu on Windows
Jupyter Notebook
Google Trendsで「インスタ映え」の検索数を取得
Google TrendsでGoogleで検索されたワードの週ごとの検索数（最大100に正規化された数）を調べることができる。これを利用して「インスタ映え」の検索数を取得する。期間は2017年の1年間。
csvの中身は以下の通り。Google Trendsの出力は日付と1週間ごとの検索数（最大100に正規化されている）。
データの確認
Google Trendsで取得したデータをグラフ化してデータを確認する。
%matplotlib inline
from IPython.core.pylabtools import figsize
import matplotlib
import matplotlib.pyplot as plt
import codecs
import pandas as pd
import numpy as np
with codecs.open('インスタ映え.csv', 'r', encoding='utf-8') as f:
# skiprows=3: データは4行目からなので3行目までをスキップ
# 2列目（dt）をインデックスに指定し、datetime型で読み込み
df = pd.read_csv(f, skiprows=3, parse_dates=True, names=['dt', 'count'], index_col='dt')
# グラフでデータを確認
count_data = df['count'].as_matrix()
n_count_data = len(count_data)
plt.subplots(figsize = (11, 9))
plt.bar(np.arange(n_count_data), count_data, color='#348ABD')
plt.xlabel('Weeks')
plt.ylabel('Count')
plt.xlim(0, n_count_data);
plt.show()
グラフは以下のようになる。ベイズ推定しなくても転換点がだいたいどの辺りかわかりそうだけど、練習なのでベイズ推定してみる。
転換点を推定する
Chapter 1: Introduction to Bayesian MethodsのExample: Inferring behaviour from text-message dataでは、1日ごとのメッセージ受信数はポアソン分布に従い、転換点tauの前後でメッセージ受信数の平均λが変化するとしてモデル化している。同様に1週間ごとの検索数がポアソン分布に従うとして同じモデルを使ってベイズ推定する。
%matplotlib inline
import matplotlib
import matplotlib.pyplot as plt
import codecs
import pandas as pd
import numpy as np
import pymc3 as pm
with codecs.open('インスタ映え.csv', 'r', encoding='utf-8') as f:
# skiprows=3: データは4行目からなので3行目までをスキップ
# 2列目（dt）をインデックスに指定し、datetime型で読み込み
df = pd.read_csv(f, skiprows=3, parse_dates=True, names=['dt', 'count'], index_col='dt')
# 検索数のデータ
count_data = df['count'].as_matrix()
# データ数
n_count_data = len(count_data)
# モデリングとサンプリング
with pm.Model() as model:
# 転換点前後の平均検索数の分布（指数分布）
# lam: Inferring behaviour from text-message dataに習って平均検索数の逆数にする
# 経験的にこう近似すると良いらしい
alpha = 1.0 / count_data.mean()
lambda_1 = pm.Exponential('lambda_1', lam=alpha)
lambda_2 = pm.Exponential('lambda_2', lam=alpha)
# 転換点（週）の分布（一様分布）
tau = pm.DiscreteUniform('tau', lower=0, upper=n_count_data - 1)
idx = np.arange(n_count_data)
# 平均検索数
# tauより小さいときはlambda_1、そうでないときはlambda_2
lambda_ = pm.math.switch(tau > idx, lambda_1, lambda_2)
# 検索数の分布（ポアソン分布）
# mu: 平均検索回数
observation = pm.Poisson('obs', mu=lambda_, observed=count_data)
# サンプルの生成
# 10000: 求めるサンプル数（バーンインは含まない）
# tune: バーンイン
# step: MCMC法のアルゴリズム
# njobs: チェーン数
trace = pm.sample(10000, tune=2000, step=pm.Metropolis(), njobs=3)
# サンプルのプロット
pm.traceplot(trace)
plt.show()
# サンプルの要約統計量
print(pm.summary(trace))
結果は以下のようになった。tau、lambda_1、lambda_2のRハットはすべて1.1より小さいので収束している。tauの事後分布の平均は28なので29週目が転換点と推定された。転換点前後の平均検索数は7.8と66.7回で、29週目（7月）の辺りで検索数が増えていると言える。この辺りから「インスタ映え」がメディアに取り上げられるようになったのだろうか。
0 件のコメント:
コメントを投稿