引入
先验分布取beta分布是个好主意（转化率在0到1之间，刚好和beta分布值域一致）。而当样本是二项分布时，参数的共轭先验分布族是beta分布，所以不需要进行MCMC。
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as stats
import pymc3 as pm
from scipy.stats import beta
from IPython.core.pylabtools import figsize
visitors_to_A = 1300
visitors_to_B = 1300
conversions_from_A = 120
conversions_from_B = 125
alpha_prior = 1
beta_prior = 1
posterior_A = beta(alpha_prior + conversions_from_A,
beta_prior + visitors_to_A - conversions_from_A)
posterior_B = beta(alpha_prior + conversions_from_B,
beta_prior + visitors_to_B - conversions_from_B)
samples = 2000
samples_posterior_A = posterior_A.rvs(samples)
samples_posterior_B = posterior_B.rvs(samples)
print((samples_posterior_A > samples_posterior_B).mean())
figsize(12.5, 4)
plt.rcParams['savefig.dpi'] = 300
plt.rcParams['figure.dpi'] = 300
x = np.linspace(0, 1, 500)
plt.plot(x, posterior_A.pdf(x), label='posterior of A')
plt.plot(x, posterior_B.pdf(x), label='posterior of B')
plt.xlim(0.05, 0.15)
plt.xlabel('value')
plt.ylabel('density')
plt.show()
第一部分
1.导入所需要的库
import pymc3 as pm
import scipy.stats as stats
import numpy as np
from IPython.core.pylabtools import figsize
import matplotlib.pyplot as plt
2.定义Model
with pm.Model() as model:
p = pm.Uniform('p', lower=0, upper=1)
p_true = 0.05
N = 1500
occurrences = stats.bernoulli.rvs(p_true, size=N)
print("What is the observed frequency in Group A? %.4f" % np.mean(occurrences))
print("Does this equal the true frequency? %s" % (np.mean(occurrences) == p_true))
PyMC3中的模型以Model类为中心。它引用了所有随机变量(RVs)并计算了模型logp及其梯度。通常，您会将其实例化为具有上下文的一部分，在这个上下文中定义的变量都被添加到这个模型。
在这里我们定义了p的一个先验分布：[0, 1]的均匀分布
然后设定p_true的值还有样本的大小
使用scipy.stats.bernoulli.rvs获得了基于伯努利分布的人工数据 用来进行模拟
两个print语句中第一个是为了获得A里的观测频率，也就是产生的人工数据的一个均值
然后进行和自设定的真实的频率对比，是否相同
而答案是False，因为是随机产生的，会在p_true的数值进行上下波动
3.定义变量
with model:
obs = pm.Bernoulli('obs', p=p, observed=occurrences)
step = pm.Metropolis()
trace = pm.sample(12000, tune=1200step=step, cores=1)
burned_trace = trace[1000:]
定义一个伯努利分布的观察变量
需要将数据传递到观察到的关键字参数中
所以讲p这个随机性变量传入obs里面，并给定观测值为occurrences
通过 observed 参数来告诉这个变量其值是已经被观测到了的（就是 obs变量），不会被拟合算法改变。
使用MCMC采样方法获得后验分布来计算进行模型拟合
使用Metropolis方法
虽然极大后验估计是一个简单快速的估计未知参数的办法，但是没有对不确定性进行估计是其缺陷。相反，比如MCMC这类基于采样的方法可以获得后验分布接近真实的采样值。为了使用MCMC采样以获得后验分布的样本，PyMC3需要指定和特定MCMC算法相关的步进方法（采样方法）,如Metropolis, Slice sampling, or the No-U-Turn Sampler (NUTS)。PyMC3中的 step_methods 子模块提供了如下采样器： NUTS, Metropolis, Slice, HamiltonianMC, and BinaryMetropolis。可以由PyMC3自动指定也可以手动指定。自动指定是根据模型中的变量类型决定的：
- 二值变量：指定 BinaryMetropolis
- 离散变量：指定 Metropolis
- 连续变量：指定 NUTS手动指定优先，可覆盖自动指定。
MCMC采样算法的主要入口点是通过pm.sample()函数。默认情况下，这个函数尝试自动分配正确的采样器，如果不传递任何内容，则自动初始化
1.在这里，我们从后部抽取样本，采样器在迭代中调整其参数。然后获取1000以后的trace数据
2.还可以使用内核kwarg并行运行多个链，可以用修改参数cores，但是以默认值或大于1的值的话会报错，因为自己的的PC机无法支持并行运行，所以使用cores=1即可
sample 函数用指定的迭代器（采样算法）抽取了1800个，收集到的采样值存储在 Trace 对象中，按照采样值获得的先后顺序存储。Trace对象是一个字典中包含了多个map。
pm.traceplot(trace)
4.画后验的图
figsize(12.5, 4)
plt.title("Posterior distribution of $p_A$, the true effectiveness of site A")
plt.vlines(p_true, 0, 90, linestyles="--", label='true $p_A$ (unknown')
plt.hist(burned_trace['p'], bins=25, histtype='stepfilled', density=True)
plt.legend()
plt.show()
每次运行的图都会一定的差别，但是大致图像就是这个样子的
pm.summary(trace)
mean sd hpd_3% hpd_97% ... ess_sd ess_bulk ess_tail r_hat
p 0.042 0.005 0.032 0.052 … 620.0 625.0 574.0 1.0
第二部分
1.导入库
import pymc3 as pm
import scipy.stats as stats
import numpy as np
from IPython.core.pylabtools import figsize
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings("ignore")
这里使用warings的库函数是为了忽略掉会显示的警告，而这些警告是对结果没有任何影响的
2.定义人工数据
true_p_A = 0.05
true_p_B = 0.04
N_A = 1500
N_B = 750
observations_A = stats.bernoulli.rvs(true_p_A, size=N_A)
observations_B = stats.bernoulli.rvs(true_p_B, size=N_B)
3.定义模型与变量
with pm.Model() as model:
p_A = pm.Uniform('p_A', 0, 1)
p_B = pm.Uniform('p_B', 0, 1)
delta = pm.Deterministic('delta', p_A-p_B)
obs_A = pm.Bernoulli('obs_A', p_A, observed=observations_A)
obs_B = pm.Bernoulli('obs_B', p_B, observed=observations_B)
step = pm.Metropolis()
trace = pm.sample(20000, step=step, cores=1)
burned_trace = trace[1000:]
这里出现了新的一个函数Deterministic，定义了一个确定性变量，用来度量A与B转化率的大小比较
也可以手工定义一个函数,道理相同
def Deterministic(p_A, p_B):
return p_A - p_B
同样使用Metropolis方法调用sample函数
4.得到A与B和差值Delata的后验分布
p_A_samples = burned_trace["p_A"]
p_B_samples = burned_trace["p_B"]
delta_samples = burned_trace["delta"]
figsize(12.5, 10)
plt.figure(1)
ax = plt.subplot(311)
plt.xlim(0, .1)
plt.hist(p_A_samples, histtype='stepfilled', bins=25, alpha=0.85,
label="posterior of $p_A$", color="#A60628", normed=True)
plt.vlines(true_p_A, 0, 80, linestyle="--", label="true $p_A$ (unknown)")
plt.legend(loc="upper right")
plt.title("Posterior distributions of $p_A$, $p_B$, and delta unknowns")
ax = plt.subplot(312)
plt.xlim(0, .1)
plt.hist(p_B_samples, histtype='stepfilled', bins=25, alpha=0.85,
label="posterior of $p_B$", color="#467821", normed=True)
plt.vlines(true_p_B, 0, 80, linestyle="--", label="true $p_B$ (unknown)")
plt.legend(loc="upper right")
ax = plt.subplot(313)
plt.hist(delta_samples, histtype='stepfilled', bins=30, alpha=0.85,
label="posterior of delta", color="#7A68A6", normed=True)
plt.vlines(true_p_A - true_p_B, 0, 60, linestyle="--",
label="true delta (unknown)")
plt.vlines(0, 0, 60, color="black", alpha=0.2)
plt.legend(loc="upper right")
plt.figure(2)
plt.xlim(0, .1)
plt.hist(p_A_samples, histtype='stepfilled', bins=30, alpha=0.80,
label="posterior of $p_A$", color='#A60628', normed=True)
plt.hist(p_B_samples, histtype='stepfilled', bins=30, alpha=0.80,
label="posterior of $p_B$", color='#467821', normed=True)
plt.legend(loc='upper right')
plt.xlabel("Value")
plt.ylabel("Density")
plt.title('Posterior distribution of $p_A$ and $p_B$')
plt.ylim(0, 80)
plt.show()
5.结果
print("Probability site A is WORSE than site B: %.3f" % np.mean(delta_samples < 0))
print("Probability site A is BETTER than site B: %.3f" % np.mean(delta_samples > 0))
判断一下A与B转化率更高的概率
Probability site A is WORSE than site B: 0.086
Probability site A is BETTER than site B: 0.914
280
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
