Configuración del algoritmo de simulación para verificar la calibración de probabilidades posteriores bayesianas

A menudo, descubrir cómo simular algo es la mejor manera de comprender los principios subyacentes. No sé exactamente cómo simular lo siguiente.

Suponer que $Y \sim N(\mu, \sigma^{2})$ y eso $\mu$ tiene una distribución previa que es $N(\gamma, \tau^{2})$. Basado en una muestra de$n$ observaciones $Y_{1}, \dots, Y_{n}$ abreviado por solo $Y$, Estoy interesado en mostrarle a un no bayesiano que la probabilidad posterior de que $\mu > 0 | Y$ está bien calibrado, por ejemplo, Prob$(\mu > 0 | P) = P$ dónde $P$es la probabilidad posterior Una discusión relacionada está aquí

Lo que realmente quiero mostrar es que si uno hiciera pruebas secuenciales y dejara de tomar muestras cuando la probabilidad posterior excede algún nivel, como 0,95, la probabilidad de que $\mu > 0$ no es $< 0.95$.

Estoy tratando de convencer a los frecuentistas de que las probabilidades bayesianas son significativas sin entrar en ninguna discusión sobre el error tipo I. Supongo que hay un problema filosófico al hablar con un frecuentista que considera hipótesis nulas en el sentido de que si lo anterior es continuo (como arriba), la probabilidad de que$\mu = 0$es cero y no se necesitan simulaciones. Agradecería algunas sugerencias sobre cómo pensar en todo el problema y cómo diseñar simulaciones de demostración. Estoy acostumbrado a hacer simulaciones frecuentas donde$\mu$solo se establece en una sola constante; Los bayesianos no condicionan$\mu$.

Para la situación secuencial establecemos un tamaño de muestra máximo posible, por ejemplo, $n=1000$.

Hay un problema sutil en el que siempre tengo problemas para pensar. Un escéptico real a veces está preocupado por un falso reclamo de efectividad ($\mu > 0$) cuando el proceso realmente no tiene exactamente efecto ($\mu=0$) La sutileza es que el escéptico está "señalando" cero como un valor especial, y tal vez está dando una probabilidad distinta de cero al evento$\mu = 0$(?) Nuestro método de demostrar que las partes posteriores están calibradas puede no hacer feliz a un escéptico porque el escéptico realmente parece querer condicionar$\mu = 0$y como bayesianos solo condicionamos lo que se puede conocer. ¿Quizás este es un caso donde la distribución previa que el estadístico está usando está en conflicto con una distribución previa discontinua que está usando el escéptico?

Los resultados de la simulación dependerán de cómo se muestrea el parámetro en la simulación. No creo que haya ninguna disputa sobre si las probabilidades posteriores se calibrarán (en el sentido de la frecuencia) si las probabilidades anteriores lo son, por lo que sospecho que una simulación no convencerá a nadie de nada nuevo.

De todos modos, en el caso de muestreo secuencial mencionado en la pregunta (tercer párrafo) se puede simular "tal cual" dibujando $\mu$ del anterior, tomando muestras dadas esto $\mu$ hasta $p(\mu>0\mid \textrm{samples})>0.95$ o se produce algún otro criterio de terminación (se necesita otro criterio de terminación ya que existe una probabilidad positiva de que la probabilidad posterior de ejecución nunca exceda $0.95$) Entonces, por cada$p(\mu>0\mid \textrm{samples})>0.95$ reclamo, compruebe si el subyacente muestreado $\mu$-parámetro es positivo y cuenta el número de positivos verdaderos frente a falsos positivos. Entonces, para$i=1,2,\ldots$:

• Muestra $\mu_i \sim N(\gamma, \tau^2)$
• por $j=1,\ldots^\ast$:
  • Muestra $y_{i,j} \sim N(\mu_i, \sigma^2)$
  • Calcular $p_{i,j} := P(\mu_i>0 \mid y_{i,1:j})$
  • Si $p_{i,j}>0.95$
    • Si $\mu_i>0$, incrementa el contador positivo verdadero
    • Si $\mu_i\leq0$, incremente el contador falso positivo
    • Romper desde el interior para bucle
  • $\ast$ alguna otra condición de ruptura, como $j\geq j_{\max}$

La relación de verdaderos positivos a todos los positivos será al menos $0.95$, que demuestra la calibración de $P(\mu>0 \mid D)>0.95$ reclamación (es.

Una implementación de Python lenta y sucia (errores muy posibles + hay un posible sesgo de detención en el que depuré hasta que vi que se mantenía la propiedad de calibración esperada).

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

tengo $0.9807$por la proporción de verdaderos positivos a todas las afirmaciones. Se acabó$0.95$ ya que la probabilidad posterior no golpeará exactamente $0.95$. Por esta razón, el código también rastrea la probabilidad posterior media "reclamada", para lo cual obtuve$0.9804$.

También se podrían cambiar los parámetros anteriores. $\gamma,\tau$ para cada $i$para demostrar una calibración "sobre todas las inferencias" (si las anteriores están calibradas). Por otro lado, uno podría realizar las actualizaciones posteriores a partir de hiperparámetros anteriores "incorrectos" (diferentes de los que se utilizan para dibujar el parámetro de verdad fundamental), en cuyo caso la calibración podría no ser válida.

Ampliando la excelente respuesta de @ juho-kokkala y usando R aquí están los resultados. Para una distribución previa para la población media mu, utilicé una mezcla igual de dos normales con media cero, una de ellas muy escéptica sobre las medias grandes.

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