Comment obtenir un postérieur d'une différence en utilisant MCMCpack?

Question

Je suis en train d'obtenir une distribution postérieure à l'aide MCMCpack d'une différence entre les deux taux de conversion, semblable à la A et B Ensemble section this PyMC tutorial.

je peux obtenir les dents postérieures des deux taux échantillonnés très bien, mais je me bats comment mettre en œuvre le delta échantillonné .. Des idées?

Modifier Le vrai delta (qui aurait été inconnu si nous n'avions pas fabriqué les données et est ce que nous voulons estimer en utilisant MCCM) est la différence entre les deux taux true_p_a et true_p_b dire 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)) 

Answer 1

vous luttez seulement parce que vous ne l'avez pas pleinement adopté un état d'esprit bayésienne. C'est tout à fait OK, j'ai eu beaucoup de problèmes conceptuels quand j'ai commencé. (Cette question est ancienne, vous l'avez probablement déjà compris).

Une densité postérieure bayésienne incorpore toutes les informations disponibles sur la distribution des paramètres du modèle. Ainsi, pour calculer une fonction de l'un des paramètres du modèle, vous calculez simplement cette fonction pour chaque tirage à partir de la distribution a posteriori. Vous n'avez pas besoin de vous soucier des erreurs standard et de l'inférence asymptotique car vous disposez déjà de toutes les informations dont vous avez besoin. Dans ce cas, parce que la différence entre les paramètres est constante et que vous avez beaucoup de données, il y a peu d'incertitude sur le delta. Il est estimé à une moyenne de 0,014 avec un SD (non SE) de 0,009.

Votre code avec l'analyse terminée:

# 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()