今回は、日毎の販売回数がポアソン分布に従うとして事後分布を推定していきます。
(1) 共役自然分布を使ったポアソン分布の事後分布推定の概要
#日毎の販売回数の集計 | |
# Convert the 'DateTime' column to datetime objects | |
bakery_data['DateTime'] = pd.to_datetime(bakery_data['DateTime']) | |
# Extract the date from the 'DateTime' column | |
bakery_data['Date'] = bakery_data['DateTime'].dt.date | |
# Group by the date and count the number of transactions per day | |
daily_sales = bakery_data.groupby('Date')['TransactionNo'].nunique() | |
# Display the first few rows of the daily sales | |
daily_sales |
from scipy.stats import gamma, poisson | |
import numpy as np | |
# Define the prior parameters for the gamma distribution (alpha and beta) | |
alpha_prior = 1 | |
beta_prior = 1 | |
# Calculate the posterior parameters | |
# For gamma-poisson conjugate, the posterior alpha is prior alpha + sum of observed counts | |
# And posterior beta is prior beta + number of observations | |
alpha_posterior = alpha_prior + daily_sales.sum() | |
beta_posterior = beta_prior + len(daily_sales) | |
# Generate a range of lambda values for plotting the posterior | |
lambda_values = np.linspace(0, daily_sales.max(), 100) | |
# Calculate the posterior gamma distribution for these lambda values | |
posterior_distribution = gamma.pdf(lambda_values, a=alpha_posterior, scale=1/beta_posterior) | |
# Plot the posterior distribution | |
import matplotlib.pyplot as plt | |
plt.figure(figsize=(10, 6)) | |
plt.plot(lambda_values, posterior_distribution, label="Posterior Distribution") | |
plt.title("Posterior Distribution of the Poisson Parameter (λ)") | |
plt.xlabel("λ (Average Number of Sales per Day)") | |
plt.ylabel("Probability Density") | |
plt.legend() | |
plt.grid(True) | |
plt.show() |
AP推定は事後分布を最大化するλの値を見つけることです。ガンマ分布の場合、a>1 のとき、λのMAP推定値は次のようになります:
# Calculate the mode (MAP estimate) of the posterior gamma distribution | |
if alpha_posterior > 1: | |
lambda_map_estimate = (alpha_posterior - 1) / beta_posterior | |
else: | |
lambda_map_estimate = 0 | |
lambda_map_estimate |
(4)事後予測分布
事後予測分布(Posterior Predictive Distribution)は、ベイズ統計学の枠組みで使用される重要な概念の1つで、ベイズ推論を行った結果から未知のデータポイントの分布を予測するための確率分布です。事後予測分布は以下のように表されます:
P(Xnew|X)=∫P(Xnew|λ)P(λ|X)dλ
ここで、ポアソン分布のパラメータ λ がガンマ分布に従う場合、事後予測分布は負の二項分布(またはポリヤ分布)になります。
負の二項分布の確率質量関数(PMF)は以下のように表されます:
P(xnew|apost,bpost,X)=Γ(xnew+apost)xnew!Γ(apost)(bpostbpost+1)apost(1bpost+1)xnew
これをやろうとすると、先ほど求めたパラメタαpost=9466,β=160(α=2,β=1)となってしまい、累乗の計算が困難です。以下のように実施します。
from scipy.stats import nbinom | |
# 新しいデータ点 x_new が取り得る値の範囲 | |
x_new_values = range(0, daily_sales.max()) | |
# 事後予測分布を計算(負の二項分布を使用) | |
p_pred = [nbinom.pmf(k=x_new, n=alpha_posterior, p=beta_posterior/(1+beta_posterior)) for x_new in x_new_values] | |
# 結果をプロット | |
import matplotlib.pyplot as plt | |
plt.bar(x_new_values, p_pred, color='blue', alpha=0.7) | |
plt.title('Posterior Predictive Distribution') | |
plt.xlabel('x_new') | |
plt.ylabel('Probability') | |
plt.show() |
最後にMCMCで事後予測分布を可視化します。
import pymc as pm | |
import arviz as az | |
import numpy as np | |
import matplotlib.pyplot as plt | |
from scipy.stats import gamma | |
# 観測データ | |
x = daily_sales.values # 観測データの例 | |
with pm.Model() as model: | |
mu = pm.Gamma('mu', alpha=2, beta=1) #事前分布 | |
obs = pm.Poisson('obs', mu=mu, observed=x) #尤度 | |
with model: | |
trace = pm.sample(3000, return_inferencedata=False) #定義されたモデルを使って、観測データxからサンプリングされるmuの値を3000個算出している | |
idata = pm.to_inference_data(trace) #複数のchainを並行してサンプリング | |
az.plot_trace(idata) #複数のchainでのmuの分布 | |
print(az.summary(idata)) #統計値の確認 | |
az.plot_posterior(idata) #muの変動の様子 | |
with model: #事後予測分布(事後分布を使って実際のデータの予測) | |
ppc = pm.sample_posterior_predictive(idata, return_inferencedata=False) #このidataからサンプルされたパラメータ(この場合はmu)を使って、新しい観測データのサンプルを生成 | |
idata_ppc = pm.to_inference_data(posterior_predictive=ppc) | |
# 記述統計量(平均・分散)によるチェック | |
ppc_samples = ppc['obs'] | |
ppc_samples = ppc_samples.reshape(-1, 50) | |
ppc_mean = ppc_samples.mean(axis=1) | |
ppc_var = ppc_samples.var(axis=1) | |
print('PPC_Mean = {0}'.format(ppc_mean)) | |
print('PPC_Var = {0}'.format(ppc_var)) | |
# 分布によるチェック | |
az.plot_ppc(idata_ppc, kind='kde', num_pp_samples=50, figsize=(12, 4)); | |
az.plot_ppc(idata_ppc, kind='kde', num_pp_samples=3000, figsize=(12, 4)); |