r - 如何使用MCMCpack获得差异的后验？

我正在尝试使用MCMCpack获得两个转换率之间差异的后验分布，类似于this PyMC tutorial.的 A和B Together 部分

我可以很好地了解两个采样率的后验，但是我在努力实现采样增量。

编辑真正的增量(如果我们没有捏造数据，这是我们要使用MCMC估算的值，则是未知的)是true_p_a和true_p_b这两个比率之差，即 0.01 。

# define true success rates
true_p_a = 0.05
true_p_b = 0.04

# set sample sizes
n_samples_a = 1000
n_samples_b = 1000

# fabricate some data
set.seed(10);
obs_a = rbinom(n=n_samples_a, size=1, prob=true_p_a)
set.seed(1);
obs_b = rbinom(n=n_samples_b, size=1, prob=true_p_b)

# what are the observed conversion rates?
mean(obs_a) #0.056
mean(obs_b) #0.042

# convert to number of successes
successes_a = sum(obs_a) #56
successes_b = sum(obs_b) #42

# calculate the posterior
require(MCMCpack)

simulations = 20000

posterior_a = MCbinomialbeta(successes_a ,n_samples_a, alpha=1, beta=1,mc=simulations)
posterior_b = MCbinomialbeta(successes_b ,n_samples_b, alpha=1, beta=1,mc=simulations)

posterior_delta = ????

posterior_density_a = density(posterior_a)
posterior_density_b = density(posterior_b)


# plot the posteriors
require(ggplot2)
ggplot() +
  geom_area(aes(posterior_density_a$x, posterior_density_a$y), fill="#7ad2f6", alpha=.5) +
  geom_vline(aes(xintercept=.05), color="#7ad2f6", linetype="dotted", size=2) +
  geom_area(aes(posterior_density_b$x, posterior_density_b$y), fill="#014d64", alpha=.5) +
  geom_vline(aes(xintercept=.04), color="#014d64", linetype="dotted", size=2) +
  scale_x_continuous(labels=percent_format(), breaks=seq(0,0.1, 0.01))

最佳答案

您之所以努力，是因为您尚未完全采用贝叶斯的思维方式。完全没问题，刚开始时我遇到了很多相同的概念问题。 (这个问题已经很老了，所以您可能已经知道了)。

贝叶斯后验密度合并了有关模型参数分布的所有可用信息。因此，要计算模型中任何参数的函数，只需从后验分布中为每次绘制计算该函数。您无需担心标准错误和渐近推断，因为您已经拥有了所需的所有信息。

在这种情况下，因为参数之间的差异是一个常数，并且您有大量数据，所以增量的不确定性很小。估计的平均值为0.014，SD(非SE)为0.009。

完成分析的代码:

  # define true success rates
  true_p_a = 0.05
  true_p_b = 0.04

  # set sample sizes
  n_samples_a = 1000
  n_samples_b = 1000

  # fabricate some data
  set.seed(10);
  obs_a = rbinom(n=n_samples_a, size=1, prob=true_p_a)
  set.seed(1);
  obs_b = rbinom(n=n_samples_b, size=1, prob=true_p_b)

  # what are the observed conversion rates?
  mean(obs_a) #0.056
  mean(obs_b) #0.042

  # convert to number of successes
  successes_a = sum(obs_a) #56
  successes_b = sum(obs_b) #42

  # calculate the posterior
  require(MCMCpack)

  simulations = 20000

  posterior_a = MCbinomialbeta(successes_a ,n_samples_a, alpha=1, beta=1,mc=simulations)
  posterior_b = MCbinomialbeta(successes_b ,n_samples_b, alpha=1, beta=1,mc=simulations)

  # Subtract the posterior deltas, look at empirical summaries and plot the empirical density function

  posterior_delta = posterior_a - posterior_b

  summary(posterior_delta)

  require(ggplot2)

  ggplot(data.frame(delta=as.numeric(posterior_delta)),aes(x=delta)) + geom_density() + theme_minimal()