Bakery Sales Dataset : (2) 日毎の販売回数をポアソン分布に従うとして、パラメタの事後分布を推定

 今回は、日毎の販売回数がポアソン分布に従うとして事後分布を推定していきます。

(1) 共役自然分布を使ったポアソン分布の事後分布推定の概要

ポアソン分布の尤度関数は次のように与えられる：

$$p(X|\lambda) = \prod_{i=1}^{n} \frac{e^{-\lambda} \lambda^{x_i}}{x_i!}$$

ここで、$X=(x_1,\dots,x_n)$ は観測データである。

ガンマ事前分布は次のように与えられる：

$$p(\lambda|a,b) = \frac{b^a}{\Gamma(a)} \lambda^{a-1} e^{-b\lambda}$$

ベイズの定理により、

$$p(\lambda|X) = \frac{p(\lambda,X)}{p(X)} = \frac{p(X|\lambda)p(\lambda|a,b)}{p(X)} \propto p(X|\lambda)p(\lambda|a,b)$$

$$= \prod_{i=1}^{n} \frac{e^{-\lambda} \lambda^{x_i}}{x_i!} \cdot \frac{b^a}{\Gamma(a)} \lambda^{a-1} e^{-b\lambda}$$

$$\propto \lambda^{\left(\sum_{i=1}^{n} x_i + a - 1\right)} e^{-(n+b)\lambda}$$

以上より、事後分布はガンマ分布となり、そのパラメータは次のようになる：

$$a_{\text{post}} = \sum_{i=1}^{n} x_i + a$$

$$b_{\text{post}} = n + b$$

したがって、$\lambda$ の事後分布は次のようになる：

$$p(\lambda|X) = \frac{b_{\text{post}}^{a_{\text{post}}}}{\Gamma(a_{\text{post}})} \lambda^{a_{\text{post}}-1} e^{-b_{\text{post}}\lambda}$$

つまり、ポアソン分布のパラメタ$\lambda$の事前分布をガンマ分布とおくと、$\lambda$の事後分布は同じくガンマ分布になる。

(2) 販売回数がポアソン分布に従うときのパラメタ推定

ここでは、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

view raw gistfile1.txt hosted with ❤ by GitHub

ここでは、事前分布としてガンマ分布を使用し、事後分布を計算してみましょう。ガンマ分布は、形状パラメータ（$\alpha$）と尺度パラメータ（$\beta$）を持ちます。これらのパラメータを適切に選ぶことが重要です。仮に、$\alpha=\beta=1$（一様事前分布に近い）として計算してみます。

	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()

view raw gistfile1.txt hosted with ❤ by GitHub

ポアソン分布のパラメータ（$\lambda$）の事後分布をプロットしました。このグラフは、観測データとガンマ分布を事前分布とした場合の$\lambda$の事後確率密度を示しています。

この事後分布から、平均販売回数の最も確からしい値や、その不確実性を評価することができます。例えば、事後分布のピークは$\lambda$の最も可能性の高い値を示し、分布の幅は不確実性の程度を示します。

(3) 得られた$\lambda$の事後分布からのMAP（Maximum A Posteriori）推定
AP推定は事後分布を最大化する$\lambda$の値を見つけることです。ガンマ分布の場合、$a > 1$ のとき、$\lambda$のMAP推定値は次のようになります： 

$$\hat{\lambda}_{\text{MAP}} = \arg\max_\lambda (p(X|\lambda) \cdot p(\lambda|a,b))$$  

$$\hat{\lambda}_{\text{MAP}} = \frac{a_{\text{post}} - 1}{b_{\text{post}}}$$  

ガンマ分布の最大値はそのモード（最頻値）となる性質があります。

実際求めてみます。

	# 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

view raw gistfile1.txt hosted with ❤ by GitHub

ポアソン分布のパラメータ（λ）の最大事後確率（MAP）推定値は約59.16です。これは、観測データをもとに、平均販売回数が一日あたり約59回であると推定されることを意味します。この値は、データと事前分布から導き出された最も確からしい販売回数の推定値です。

今度は、$\alpha=2, \beta=1$でやって、先ほどの結果と比較してみます。新しい事前分布に基づいた$\lambda$の最大事後確率（MAP）推定値は約59.16で同じでした。これは、事前分布を変更しても、観測データの影響が大きいため、MAP推定値に大きな変化が見られないことを示しています。この結果は、データが豊富であればあるほど、事前分布の影響が小さくなるというベイズ推定の特性を反映しています。

(4)事後予測分布

事後予測分布（Posterior Predictive Distribution）は、ベイズ統計学の枠組みで使用される重要な概念の1つで、ベイズ推論を行った結果から未知のデータポイントの分布を予測するための確率分布です。事後予測分布は以下のように表されます: 

$$P(X_{\text{new}}|X) = \int P(X_{\text{new}}|\lambda)P(\lambda|X)d\lambda$$  

ここで、ポアソン分布のパラメータ $\lambda$ がガンマ分布に従う場合、事後予測分布は負の二項分布（またはポリヤ分布）になります。 

負の二項分布の確率質量関数（PMF）は以下のように表されます: 

$$P(x_{\text{new}}|a_{\text{post}}, b_{\text{post}}, X) = \frac{\Gamma(x_{\text{new}} + a_{\text{post}})}{x_{\text{new}}!\Gamma(a_{\text{post}})} \left(\frac{b_{\text{post}}}{b_{\text{post}} + 1}\right)^{a_{\text{post}}} \left(\frac{1}{b_{\text{post}} + 1}\right)^{x_{\text{new}}}$$

これをやろうとすると、先ほど求めたパラメタ$\alpha_post=9466, \beta=160 (\alpha=2, \beta=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()

view raw gistfile1.txt hosted with ❤ by GitHub

最後に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));

view raw gistfile1.txt hosted with ❤ by GitHub

ここまでのソースコード