Configuration d'un algorithme de simulation pour vérifier l'étalonnage des probabilités postérieures bayésiennes

Déterminer comment simuler quelque chose est souvent le meilleur moyen de comprendre les principes sous-jacents. Je ne sais pas exactement comment simuler ce qui suit.

Supposer que $Y \sim N(\mu, \sigma^{2})$ et cela $\mu$ a une distribution antérieure qui est $N(\gamma, \tau^{2})$. Basé sur un échantillon de$n$ observations $Y_{1}, \dots, Y_{n}$ abrégé par juste $Y$, Je souhaite montrer à un non-bayésien que la probabilité postérieure que $\mu > 0 | Y$ est bien calibré, par exemple, Prob$(\mu > 0 | P) = P$ où $P$est la probabilité postérieure. Une discussion connexe est ici

Ce que je veux vraiment montrer, c'est que si l'on devait faire des tests séquentiels et arrêter l'échantillonnage lorsque la probabilité postérieure dépasse un certain niveau tel que 0,95, la probabilité que $\mu > 0$ n'est pas $< 0.95$.

J'essaie de convaincre les habitués que les probabilités bayésiennes sont significatives sans entrer dans une discussion sur l'erreur de type I. Je suppose qu'il y a un problème philosophique lorsque l'on parle à un fréquentateur qui nourrit des hypothèses nulles en ce que si le prieur est continu (comme ci-dessus) la probabilité que$\mu = 0$est nul et aucune simulation n'est nécessaire. J'apprécierais quelques suggestions sur la façon de penser à l'ensemble du problème et de concevoir des simulations de démonstration. J'ai l'habitude de faire des simulations fréquentistes où$\mu$est juste réglé sur une seule constante; Les Bayésiens ne conditionnent pas$\mu$.

Pour la situation séquentielle, nous définissons une taille d'échantillon maximale possible, par exemple, $n=1000$.

Il y a une subtilité au problème à laquelle j'ai toujours du mal à penser. Un vrai sceptique s'inquiète parfois d'une fausse déclaration d'efficacité ($\mu > 0$) lorsque le processus n'a vraiment aucun effet ($\mu=0$). La subtilité est que le sceptique "choisit" zéro comme valeur spéciale, et donne peut-être une probabilité non nulle à l'événement$\mu = 0$(?). Notre méthode de montrer que les postérieurs sont calibrés peut ne pas rendre un tel sceptique heureux parce que le sceptique semble vraiment vouloir se conditionner$\mu = 0$et en tant que Bayésiens, nous ne conditionnons que ce qui est connaissable. Peut-être s'agit-il d'un cas où la distribution antérieure que le statisticien utilise entre en conflit avec une distribution antérieure discontinue utilisée par le sceptique?

Les résultats de la simulation dépendront de la façon dont le paramètre est échantillonné dans la simulation. Je ne pense pas qu'il y ait de contestation quant à savoir si les probabilités postérieures seront calibrées (au sens de la fréquence) si les probabilités antérieures le sont, donc je soupçonne qu'une simulation ne convaincra personne de quelque chose de nouveau.

Quoi qu'il en soit, dans le cas d'échantillonnage séquentiel mentionné dans la question (troisième paragraphe) peut être simulé "tel quel" en dessinant $\mu$ de la précédente, en tirant des échantillons compte tenu de cette $\mu$ jusqu'à $p(\mu>0\mid \textrm{samples})>0.95$ ou un autre critère de terminaison se produit (un autre critère de terminaison est nécessaire car il existe une probabilité positive que la probabilité de fonctionnement postérieur ne dépasse jamais $0.95$). Ensuite, pour chaque$p(\mu>0\mid \textrm{samples})>0.95$ vérifier si le sous-jacent échantillonné $\mu$-le paramètre est positif et compte le nombre de vrais positifs par rapport aux faux positifs. Donc pour$i=1,2,\ldots$:

• Échantillon $\mu_i \sim N(\gamma, \tau^2)$
• Pour $j=1,\ldots^\ast$:
  • Échantillon $y_{i,j} \sim N(\mu_i, \sigma^2)$
  • Calculer $p_{i,j} := P(\mu_i>0 \mid y_{i,1:j})$
  • Si $p_{i,j}>0.95$
    • Si $\mu_i>0$, incrémenter le vrai compteur positif
    • Si $\mu_i\leq0$, incrémenter le compteur de faux positifs
    • Rompre de l'intérieur pour la boucle
  • $\ast$ une autre condition de rupture, telle que $j\geq j_{\max}$

Le rapport des vrais positifs à tous les positifs sera au moins $0.95$, qui démontre l'étalonnage du $P(\mu>0 \mid D)>0.95$ réclamations.

Une implémentation Python lente et sale (bogues très possibles + il y a un biais d'arrêt potentiel en ce que j'ai débogué jusqu'à ce que je vois la propriété de calibrage attendue).

# (C) Juho Kokkala 2016
# MIT License 

import numpy as np

np.random.seed(1)

N = 10000
max_samples = 50

gamma = 0.1
tau = 2
sigma = 1

truehits = 0
falsehits = 0

p_positivemus = []

while truehits + falsehits < N:
    # Sample the parameter from prior
    mu = np.random.normal(gamma, tau)

    # For sequential updating of posterior
    gamma_post = gamma
    tau2_post = tau**2

    for j in range(max_samples):
        # Sample data
        y_j = np.random.normal(mu, sigma)

        gamma_post = ( (gamma_post/(tau2_post) + y_j/(sigma**2)) /
                       (1/tau2_post + 1/sigma**2) )
        tau2_post = 1 / (1/tau2_post + 1/sigma**2)

        p_positivemu = 1 - stats.norm.cdf(0, loc=gamma_post,
                                          scale=np.sqrt(tau2_post))

        if p_positivemu > 0.95:
            p_positivemus.append(p_positivemu)
            if mu>0:
                truehits += 1
            else:
                falsehits +=1
            if (truehits+falsehits)%1000 == 0:
                print(truehits / (truehits+falsehits))
                print(truehits+falsehits)
            break

print(truehits / (truehits+falsehits))
print(np.mean(p_positivemus))

j'ai eu $0.9807$pour la proportion de vrais positifs par rapport à toutes les allégations. C'est terminé$0.95$ car la probabilité postérieure ne frappera pas exactement $0.95$. Pour cette raison, le code suit également la probabilité postérieure "revendiquée" moyenne, pour laquelle j'ai obtenu$0.9804$.

On pourrait également changer les paramètres antérieurs $\gamma,\tau$ pour chaque $i$pour démontrer un étalonnage "sur toutes les inférences" (si les priors sont étalonnés). D'un autre côté, on pourrait effectuer les mises à jour postérieures à partir de «mauvais» hyperparamètres antérieurs (différents de ceux utilisés pour dessiner le paramètre de vérité terrain), auquel cas l'étalonnage pourrait ne pas tenir.

En développant l'excellente réponse de @ juho-kokkala et en utilisant R, voici les résultats. Pour une distribution antérieure de la moyenne de la population mu, j'ai utilisé un mélange égal de deux normales avec une moyenne nulle, l'une d'entre elles étant très sceptique quant aux moyennes importantes.

## Posterior density for a normal data distribution and for
## a mixture of two normal priors with mixing proportions wt and 1-wt
## and means mu1 mu2 and variances v1 an
## Adapted for LearnBayes package normal.normal.mix function

## Produces a list of 3 functions.  The posterior density and cum. prob.
## function can be called with a vector of posterior means and variances
## if the first argument x is a scalar

mixpost <- function(stat, vstat, mu1=0, mu2=0, v1, v2, wt) {
  if(length(stat) + length(vstat) != 2) stop('improper arguments')
  probs      <- c(wt, 1. - wt)
  prior.mean <- c(mu1, mu2)
  prior.var  <- c(v1,  v2)

  post.precision <- 1. / prior.var + 1. / vstat
  post.var       <- 1. / post.precision
  post.mean <- (stat / vstat + prior.mean / prior.var) / post.precision
  pwt       <- dnorm(stat, prior.mean, sqrt(vstat + prior.var))
  pwt       <- probs * pwt / sum(probs * pwt)

  dMix <- function(x, pwt, post.mean, post.var)
    pwt[1] * dnorm(x, mean=post.mean[1], sd=sqrt(post.var[1])) +
    pwt[2] * dnorm(x, mean=post.mean[2], sd=sqrt(post.var[2]))
  formals(dMix) <- z <-
    list(x=NULL, pwt=pwt, post.mean=post.mean, post.var=post.var)

  pMix <- function(x, pwt, post.mean, post.var)
    pwt[1] * pnorm(x, mean=post.mean[1], sd=sqrt(post.var[1])) +
    pwt[2] * pnorm(x, mean=post.mean[2], sd=sqrt(post.var[2]))
  formals(pMix) <- z

  priorMix <- function(x, mu1, mu2, v1, v2, wt)
    wt * dnorm(x, mean=mu1, sd=sqrt(v1)) +
    (1. - wt) * dnorm(x, mean=mu2, sd=sqrt(v2))
  formals(priorMix) <- list(x=NULL, mu1=mu1, mu2=mu2, v1=v1, v2=v2, wt=wt)
  list(priorMix=priorMix, dMix=dMix, pMix=pMix)
}

## mixposts handles the case where the posterior distribution function
## is to be evaluated at a scalar x for a vector of point estimates and
## variances of the statistic of interest
## If generates a single function

mixposts <- function(stat, vstat, mu1=0, mu2=0, v1, v2, wt) {
  post.precision1 <- 1. / v1 + 1. / vstat
  post.var1       <- 1. / post.precision1
  post.mean1      <- (stat / vstat + mu1 / v1) / post.precision1

  post.precision2 <- 1. / v2 + 1. / vstat
  post.var2       <- 1. / post.precision2
  post.mean2      <- (stat / vstat + mu2 / v2) / post.precision2

  pwt1 <- dnorm(stat, mean=mu1, sd=sqrt(vstat + v1))
  pwt2 <- dnorm(stat, mean=mu2, sd=sqrt(vstat + v2))
  pwt <- wt * pwt1 / (wt * pwt1 + (1. - wt) * pwt2)

  pMix <- function(x, post.mean1, post.mean2, post.var1, post.var2, pwt)
    pwt        * pnorm(x, mean=post.mean1, sd=sqrt(post.var1)) +
    (1. - pwt) * pnorm(x, mean=post.mean2, sd=sqrt(post.var2))
  formals(pMix) <-
    list(x=NULL, post.mean1=post.mean1, post.mean2=post.mean2,
         post.var1=post.var1, post.var2=post.var2, pwt=pwt)
 pMix
}

## Compute proportion mu > 0 in trials for
## which posterior prob(mu > 0) > 0.95, and also use a loess smoother
## to estimate prob(mu > 0) as a function of the final post prob
## In sequential analyses of observations 1, 2, ..., N, the final
## posterior prob is the post prob at the final sample size if the
## prob never exceeds 0.95, otherwise it is the post prob the first
## time it exceeds 0.95

sim <- function(N, prior.mu=0, prior.sd, wt, mucut=0, postcut=0.95,
                nsim=1000, plprior=TRUE) {
  prior.mu <- rep(prior.mu, length=2)
  prior.sd <- rep(prior.sd, length=2)
  sd1 <- prior.sd[1]; sd2 <- prior.sd[2]
  v1 <- sd1 ^ 2
  v2 <- sd2 ^ 2
  if(plprior) {
    pdensity <- mixpost(1, 1, mu1=prior.mu[1], mu2=prior.mu[2],
                        v1=v1, v2=v2, wt=wt)$priorMix
    x <- seq(-3, 3, length=200)
    plot(x, pdensity(x), type='l', xlab=expression(mu), ylab='Prior Density')
    title(paste(wt, 1 - wt, 'Mixture of Zero Mean Normals\nWith SD=',
                round(sd1, 3), 'and', round(sd2, 3)))
  }
  j <- 1 : N
  Mu <- Post <- numeric(nsim)
  stopped <- integer(nsim)

  for(i in 1 : nsim) {
    # See http://stats.stackexchange.com/questions/70855
    component <- sample(1 : 2, size=1, prob=c(wt, 1. - wt))
    mu <- prior.mu[component] + rnorm(1) * prior.sd[component]
    # mu <- rnorm(1, mean=prior.mu, sd=prior.sd) if only 1 component

    Mu[i] <- mu
    y  <- rnorm(N, mean=mu, sd=1)
    ybar <- cumsum(y) / j
    pcdf <- mixposts(ybar, 1. / j, mu1=prior.mu[1], mu2=prior.mu[2],
                     v1=v1, v2=v2, wt=wt)
    if(i==1) print(body(pcdf))
    post    <- 1. - pcdf(mucut)
    Post[i] <- if(max(post) < postcut) post[N]
               else post[min(which(post >= postcut))]
    stopped[i] <- if(max(post) < postcut) N else min(which(post >= postcut))
  }
  list(mu=Mu, post=Post, stopped=stopped)
}

# Take prior on mu to be a mixture of two normal densities both with mean zero
# One has SD so that Prob(mu > 1) = 0.1
# The second has SD so that Prob(mu > 0.25) = 0.05
prior.sd <- c(1 / qnorm(1 - 0.1), 0.25 / qnorm(1 - 0.05))
prior.sd
set.seed(2)
z <- sim(500, prior.mu=0, prior.sd=prior.sd, wt=0.5, postcut=0.95, nsim=10000)

mu   <- z$mu
post <- z$post
st   <- z$stopped
plot(mu, post)
abline(v=0, col=gray(.8)); abline(h=0.95, col=gray(.8))
hist(mu[post >= 0.95], nclass=25)
k <- post >= 0.95
mean(k)   # 0.44 of trials stopped with post >= 0.95
mean(st)  # 313 average sample size
mean(mu[k] > 0)  # 0.963 of trials with post >= 0.95 actually had mu > 0
mean(post[k])    # 0.961 mean posterior prob. when stopped early
w <- lowess(post, mu > 0, iter=0)
# perfect calibration of post probs 
plot(w, type='n',         # even if stopped early
     xlab=expression(paste('Posterior Probability ', mu > 0, ' Upon Stopping')),
     ylab=expression(paste('Proportion of Trials with ',  mu > 0)))
abline(a=0, b=1, lwd=6, col=gray(.85))
lines(w)