本文是學(xué)習(xí)《貝葉斯方法 概率編程與貝葉斯推斷》重點(diǎn)參考了這個(gè)教程
import numpy as np
import pymc3 as pm
import matplotlib.pyplot as plt
import os
# 加載數(shù)據(jù)并且展示
count_data=np.loadtxt("data/txtdata.csv")
plt.figure(figsize=(12.5,3.5))
n_count_data=len(count_data)
plt.bar(np.arange(n_count_data),count_data,color="#348ABD")
plt.xlabel("days")
plt.ylabel("number of messages")
plt.title("did the model change?")
plt.xlim(0,n_count_data)
plt.show()
bayesian_opitimization_2_0.png
初步假設(shè)财松，在時(shí)間之前短信數(shù)量
服從泊松分布
,在時(shí)間
之后短信數(shù)量
服從泊松分布
.這樣就有了三個(gè)參數(shù)
次氨，
，
母市，我們的目標(biāo)呢就是根據(jù)證據(jù)給出這三個(gè)參數(shù)的概率分布岸霹。（注意栖忠！是參數(shù)的概率分布ｒ蚪洹召烂！這是貝葉斯推斷的核心＜罟ぁ！）
為了進(jìn)行推斷呢奏夫，我們首先給這三個(gè)參數(shù)一個(gè)初始的先驗(yàn)分布怕篷，因?yàn)橄闰?yàn)分布是要根據(jù)證據(jù)不斷更新的，相當(dāng)于一個(gè)初始值酗昼，會(huì)在后續(xù)的調(diào)整中不斷地接近真實(shí)值廊谓，所以我們先驗(yàn)可以給的隨意一些。是個(gè)離散值麻削，最常見是用一個(gè)離散的均勻分布就完事了蒸痹，
春弥，
是連續(xù)值，我們用一個(gè)指數(shù)分布
叠荠。有人就說啦匿沛，憑啥呀？為啥用這些做先驗(yàn)?zāi)亻欢Γ繘]有憑啥俺祠，我們只是想給一個(gè)先驗(yàn)分布而已，你要是想用別的分布完全也可以借帘，甚至蜘渣，你可以不用已知的分布，直接一個(gè)個(gè)值的指定肺然，也沒毛病蔫缸，自創(chuàng)一個(gè)分布。都行际起。
接下來的東西呢拾碌，為了方便編程，我們用到了Pymc3這個(gè)庫街望，這是一個(gè)專門針對(duì)概率相關(guān)問題進(jìn)行編程的庫校翔。還是那句話，用了方便灾前，不用也行
# 先驗(yàn)分布
with pm.Model() as model:
alpha=1.0/count_data.mean()
lambda_1=pm.Exponential("lambda_1",alpha)
lambda_2=pm.Exponential("lambda_2",alpha)
tau=pm.DiscreteUniform("tau",lower=0,upper=n_count_data)
# get likelihood
lbds=np.zeros(n_count_data)
lbds[:tau]=lambda_1
lbds[tau:]=lambda_2
target=pm.Poisson("target",observed=lbds)
# 學(xué)習(xí)方法1： MAP防症。返回的是一個(gè)參數(shù)的數(shù)值，而不是分布哎甲！ 這玩意就是極大似然求導(dǎo)吧（不是蔫敲，要區(qū)分MLE和MAP），但是這個(gè)tau是怎么求出來的炭玫，不可導(dǎo)吧奈嘿。 參考這篇文章https://zhuanlan.zhihu.com/p/40024110
map_estimate=pm.find_MAP(model=model)
print(map_estimate)
{'lambda_1_log__': array(-9.1772865), 'lambda_2_log__': array(-9.1772865), 'tau': array(37), 'lambda_1': array(0.00010336), 'lambda_2': array(0.00010336)}
# 學(xué)習(xí)方法2：MCMC
import arviz as az
with model:
trace=pm.sample(500,return_inferencedata=False)
az.plot_trace(trace)
display(az.summary(trace,round_to=2))
Multiprocess sampling (4 chains in 4 jobs)
CompoundStep
>NUTS: [lambda_2, lambda_1]
>Metropolis: [tau]
Sampling 4 chains for 1_000 tune and 500 draw iterations (4_000 + 2_000 draws total) took 10 seconds.
There was 1 divergence after tuning. Increase `target_accept` or reparameterize.
The number of effective samples is smaller than 25% for some parameters.
bayesian_opitimization_6_4.png
