はじめに

Pythonのベイズ推定ライブラリPyMC3を用いて、二峰性のある離散確率分布からポアソン分布の期待値$\lambda_1,\lambda_2$、並びに各峰に属する比率を推定します。

動作環境&使用ライブラリ

動作環境

#OS
import platform
print(platform.platform())
#Windows-10-10.0.15063-SP0

#pythonバージョン
print(platform.python_version())
#version 3.7.7

#ライブラリ
import numpy as np #version 1.19.1
import pandas as pd #version 1.1.0
from scipy.stats import poisson #version 1.5.0
import matplotlib.pyplot as plt #version 3.2.2
import seaborn as sns #version 0.10.1
import pymc3 as pm #version 3.8
import theano.tensor as T #1.0.4

やってみる

ベルヌーイ試行回数はいずれも同じ回数として、成功回数期待値$\lambda$が3と10の郡が混ざっている離散型確率分布を考えます。

#期待値3を900、期待値10を1000回試行したときの成功回数をランダムに生成します
obs = np.hstack((poisson.rvs(3,size=900), poisson.rvs(10,size=1000)))

sns.countplot(obs)

出力結果

(これの二峰の期待値$\lambda$と比率を推定したい)

PyMC3のモデルを考える

期待値3の郡をA、期待値10の郡をBとした時、
観測した値$n_x$がどちらの郡に属しているのか事前の情報ではわからない為、
観測した値$n_x$がAである確率$p(A)$は0~1の一様分布に従うとします。
そうすると、観測した値$n_x$がBである確率$p(B)$は$1 - p(A)$ということになります。

更に、A,B真の期待値も分からないとしてAとBの期待値を推定します。
事前分布は指数分布に従うとし、指数分布の期待値$\frac{1}{\lambda}$の$\lambda$に観測した値の平均値を当てはめます。この時、二峰あることも考えます。
事前分布にて得られた値を期待値mu,観測値を先程ランダムに生成した値として、事後分布を推定します。

with pm.Model() as model:
    p1 = pm.Uniform('p1',0,1) #p(A)である確率
    p2 = 1 - p1 #p(B)である確率
    p = T.stack([p1, p2])

    #p1=0, p2=1とするカテゴリ変数をサンプル数だけ用意します。
    assignment = pm.Categorical('assignment', p, shape=obs.shape[0], testval = np.random.randint(0, 2, obs.shape[0]))

    # 事前分布を 期待値 = 1 / 観測値平均とした 指数分布として設定
    lambda_ = pm.Exponential('lambda_', 1./obs.mean(), shape=2)
    lambda_i = pm.Deterministic('lambda_i', lambda_[assignment])
    # muをlambda_ 生成した値をobservedに設定
    observations = pm.Poisson('obs', mu=lambda_i, observed=obs)

    # 探索
    step1 = pm.Metropolis(vars=[p, lambda_]) #離散的な値はMetropolisに設定
    step2 = pm.ElemwiseCategorical(vars=[assignment])
    trace = pm.sample(5000, step=[step1, step2])

# 結果をプロット
plt.figure(figsize=(12, 8))

plt.subplot(211)
plt.plot(trace['lambda_'][:,1], label='A郡', c='r')
plt.plot(trace['lambda_'][:,0], label='B郡', c='c')
plt.ylim(2, 11)
plt.xlabel('探索回数')
plt.ylabel('期待値')
plt.title('推定結果')
leg = plt.legend(loc='upper right')
leg.get_frame().set_alpha(0.7)

plt.subplot(212)
plt.plot(1 - trace['p1'], label='A郡に選ばれたサンプルの割合', c='r')
plt.legend(loc='upper right')
plt.xlabel('探索回数')
plt.ylabel('割合(%)')

plt.show()plt.figure(figsize=(12, 8))

# 1000回目以降の探索平均値を出力
trace_lambda_1 = trace['lambda_'][:,1][1000:].mean()
trace_lambda_2 = trace['lambda_'][:,0][1000:].mean()
p1 = (1 - trace['p1'][1000:]).mean()

print('A郡 期待値推定結果 %.2f' % trace_lambda_1)
print('B郡 期待値推定結果 %.2f' % trace_lambda_2)
print('A郡である割合 %.3f' % p1)

A郡 期待値探索結果 3.01
B郡 期待値探索結果 9.84
A郡 である割合 0.461

結果

	真の値	推定結果
A郡期待値	3	3.01
B郡期待値	10	9.84
A郡割合	0.473	0.46

となりました。

ポアソン分布期待値が3.01と9.84
サンプルが46:54の比で構成されている離散型確率分布は以下となります。

# 推定結果の期待値にて0.001~0.999の間の確率を取る成功回数k
a_obs = np.arange(poisson.ppf(0.001, trace_lambda_1),
             poisson.ppf(0.999, trace_lambda_1))
b_obs = np.arange(poisson.ppf(0.001, trace_lambda_2),
             poisson.ppf(0.999, trace_lambda_2))
# 期待値3,10にて0.001~0.999の間の確率を取る成功回数k
a_true = np.arange(poisson.ppf(0.001, 3),
             poisson.ppf(0.999, 3))
b_true = np.arange(poisson.ppf(0.001, 10),
             poisson.ppf(0.999, 10))
# 推定結果での成功回数kの確率
A_obs = pd.DataFrame({'lambda': a_obs,
                 'percent': poisson.pmf(a_obs, trace_lambda_1) * (46 / 100)})
B_obs = pd.DataFrame({'lambda': b_obs,
                'percent': poisson.pmf(b_obs, trace_lambda_2) * (54 / 100)} )
# 期待値3,10での成功回数kの確率
A_true = pd.DataFrame({'lambda': a_true,
                 'percent': poisson.pmf(a_true, 3) * (9 / 19)})
B_true = pd.DataFrame({'lambda': b_true,
                'percent': poisson.pmf(b_true, 10) * (10 / 19)} )

#ランダムに生成したものと比べてみる
plt.figure(figsize=(20, 6) )
plt.subplot(121)
sns.barplot(data=pd.merge(A_obs,B_obs,how='outer', on=['lambda','percent']).groupby('lambda', as_index=False)['percent'].sum(),
            x='lambda', y='percent')
plt.title('推定結果')
plt.ylabel('確率')
plt.xlabel('成功回数k')
plt.xticks(np.int8(np.arange(0,21)),np.int8(np.arange(0,21)))
plt.subplot(122)
sns.barplot(data=pd.merge(A_true,B_true,how='outer', on=['lambda','percent']).groupby('lambda', as_index=False)['percent'].sum(),
            x='lambda', y='percent')
plt.title('期待値3,期待値10 比率9:10 の確率分布')
plt.ylabel('確率')
plt.xlabel('成功回数k')
plt.xticks(np.int8(np.arange(0,21)),np.int8(np.arange(0,21)))
plt.show()

二峰性のある分布の期待値と比率を推定することができました。