3

我正在尝试使用MCMCpack两个转换率之间的差异来获得后验分布，类似于本 PyMC 教程的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))
4

1 回答 1

3

你只是在挣扎，因为你还没有完全采用贝叶斯思维方式。没关系，刚开始时我有很多相同的概念问题。（这个问题已经很老了，所以你可能已经想通了）。

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

在这种情况下，由于参数之间的差异是一个常数，并且您有大量数据，因此 delta 几乎没有不确定性。估计平均值为 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()
于 2017-03-30T18:54:56.427 回答