Configurando algoritmo de simulação para verificar a calibração de probabilidades posteriores bayesianas

Descobrir como simular algo geralmente é a melhor maneira de entender os princípios subjacentes. Estou um pouco sem saber exatamente como simular o seguinte.

Suponha que e tenham uma distribuição anterior que seja . Com base em uma amostra de observações abreviadas por apenas , estou interessado em mostrar a um não Bayesiano que a probabilidade posterior de que está bem calibrado, por exemplo, Prob onde é a probabilidade posterior. Uma discussão relacionada está aquiμ N ( γ , τ 2 ) n Y 1 , … , Y n Y μ > 0 | Y ( μ > 0 | P ) = P $Y \sim N(\mu, \sigma^{2})$$\mu$$N(\gamma, \tau^{2})$$n$$Y_{1}, \dots, Y_{n}$$Y$$\mu > 0 | Y$$(\mu > 0 | P) = P$$P$

O que realmente quero mostrar é que, se alguém fizer testes sequenciais e parar a amostragem quando a probabilidade posterior exceder algum nível, como 0,95, a probabilidade de não será .< $\mu > 0$$< 0.95$

Estou tentando convencer os freqüentadores de que as probabilidades bayesianas são significativas sem entrar em nenhuma discussão sobre o erro do tipo I. Suponho que exista um problema filosófico ao conversar com um freqüentador que tenha hipóteses nulas, pois se o anterior for contínuo (como acima), a probabilidade de que é zero e simulações não são necessárias. Gostaria de receber algumas sugestões sobre como pensar em todo o problema e como projetar simulações de demonstração. Estou acostumado a fazer simulações freqüentes, onde é apenas definido como uma única constante; Os bayesianos não se condicionam a .μ $\mu = 0$$\mu$$\mu$

Para a situação seqüencial, definimos o tamanho máximo possível da amostra, por exemplo, .$n=1000$

Há uma sutileza no problema em que sempre tenho problemas para pensar. Um cético real às vezes se preocupa com uma alegação falsa de eficácia ( ) quando o processo realmente não tem exatamente nenhum efeito ( ). A sutileza é que o cético está "destacando" zero como um valor especial e talvez esteja dando probabilidade diferente de zero ao evento (?). Nosso método de mostrar que as partes posteriores são calibradas pode não deixar esse cético feliz porque o cético realmente parece querer condicionar em e, como bayesianos, apenas condicionamos o que é conhecido. Talvez este seja um caso em que a distribuição anterior que o estatístico está usando conflite com uma distribuição anterior descontínua que o cético está usando?μ = 0 μ = 0 μ = $\mu > 0$$\mu=0$$\mu = 0$$\mu = 0$

Os resultados da simulação dependerão de como o parâmetro é amostrado na simulação. Eu não acho que exista uma disputa sobre se as probabilidades posteriores serão calibradas (no sentido de frequência) se as probabilidades anteriores forem, então eu suspeito que uma simulação não convencerá ninguém de algo novo.

De qualquer forma, no caso de amostragem seqüencial mencionado na pergunta (terceiro parágrafo) pode ser simulado "como está", desenhando do anterior, desenhando amostras fornecidas até que ou ocorre algum outro critério de terminação (é necessário outro critério de terminação, pois existe uma probabilidade positiva de que a probabilidade posterior em execução nunca exceda ). Então, para cada afirmação , verifique se o parâmetro amostrado subjacente é positivo e conte o número de positivos verdadeiros versus falsos positivos. Então, para :$\mu$$\mu$$p(\mu>0\mid \textrm{samples})>0.95$$0.95$$p(\mu>0\mid \textrm{samples})>0.95$$\mu$$i=1,2,\ldots$

• Exemplo$\mu_i \sim N(\gamma, \tau^2)$
• Para : $j=1,\ldots^\ast$
  • Amostra$y_{i,j} \sim N(\mu_i, \sigma^2)$
  • Calcular$p_{i,j} := P(\mu_i>0 \mid y_{i,1:j})$
  • Se$p_{i,j}>0.95$
    • Se , aumente o contador positivo verdadeiro$\mu_i>0$
    • Se , aumente o contador de falso positivo$\mu_i\leq0$
    • Quebra do loop for interno
  • j ≥ j $\ast$ alguma outra condição de interrupção, como$j\geq j_{\max}$

A proporção de verdadeiros positivos para todos os positivos será de pelo menos , o que demonstra a calibração das reivindicações .$0.95$$P(\mu>0 \mid D)>0.95$

Uma implementação Python lenta e suja (bugs muito possíveis + existe um viés de parada potencial em que depurei até ver a propriedade de calibração 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))

Recebi para a proporção de verdadeiros positivos para todas as reivindicações. Isso é superior a pois a probabilidade posterior não atinge exatamente . Por esse motivo, o código também rastreia a probabilidade posterior "reivindicada" média, pela qual obtive .$0.9807$$0.95$$0.95$$0.9804$

Pode-se também alterar os parâmetros anteriores para cada demonstrar uma calibração "em todas as inferências" (se os anteriores forem calibrados). Por outro lado, pode-se realizar as atualizações posteriores a partir de hiperparâmetros anteriores "errados" (diferentes do que é usado no desenho do parâmetro ground-truth); nesse caso, a calibração pode não se manter.$\gamma,\tau$$i$

Expandindo a excelente resposta de @ juho-kokkala e usando R aqui estão os resultados. Para uma distribuição prévia da população, a média mu utilizou uma mistura igual de duas normais com média zero, uma delas muito cética em relação às médias 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)