Menyiapkan algoritma simulasi untuk memeriksa kalibrasi probabilitas posterior Bayesian

Mencari tahu bagaimana mensimulasikan sesuatu sering kali merupakan cara terbaik untuk memahami prinsip-prinsip yang mendasarinya. Saya sedikit bingung bagaimana tepatnya mensimulasikan yang berikut ini.

Misalkan dan bahwa memiliki distribusi sebelumnya yaitu . Berdasarkan sampel pengamatan disingkat dengan hanya , saya tertarik untuk menunjukkan kepada non-Bayesian bahwa probabilitas posterior bahwa dikalibrasi dengan baik, misalnya, Prob mana adalah probabilitas posterior. Diskusi terkait ada di sini$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$

Yang saya benar-benar ingin tunjukkan adalah bahwa jika seseorang melakukan pengujian sekuensial dan menghentikan pengambilan sampel ketika probabilitas posterior melebihi beberapa level seperti 0,95 probabilitas bahwa tidak .$\mu > 0$$< 0.95$

Saya mencoba untuk meyakinkan sering bahwa probabilitas Bayesian bermakna tanpa masuk ke diskusi tentang kesalahan tipe I. Saya kira ada masalah filosofis ketika berbicara dengan seorang frequentist yang menghibur hipotesis nol dalam bahwa jika sebelumnya adalah berkesinambungan (seperti di atas) probabilitas bahwa adalah nol dan simulasi tidak diperlukan. Saya akan menghargai beberapa saran tentang bagaimana memikirkan seluruh masalah dan bagaimana merancang simulasi demonstrasi. Saya terbiasa melakukan simulasi frequentist di mana hanya disetel ke konstanta tunggal; Bayesians tidak mengkondisikan pada .$\mu = 0$$\mu$$\mu$

Untuk situasi sekuensial kami menetapkan ukuran sampel maksimum yang mungkin, misalnya, .$n=1000$

Ada sedikit masalah dengan masalah yang selalu saya pikirkan. Seorang skeptis sejati kadang-kadang khawatir tentang klaim salah efektif ( ) ketika proses tersebut benar-benar tidak memiliki efek ( ). Subltynya adalah bahwa orang yang skeptis "memilih" nol sebagai nilai khusus, dan mungkin memberikan probabilitas bukan nol pada peristiwa (?). Metode kami untuk menunjukkan bahwa posisi dikalibrasi mungkin tidak membuat senang skeptis karena skeptis benar-benar tampaknya ingin mengkondisikan pada dan sebagai Bayesians kami hanya mengkondisikan pada apa yang bisa diketahui. Mungkin ini adalah kasus di mana distribusi sebelumnya yang menggunakan statistik konflik dengan distribusi sebelumnya terputus yang skeptis gunakan?$\mu > 0$$\mu=0$$\mu = 0$$\mu = 0$

Hasil simulasi akan tergantung pada bagaimana parameter sampel dalam simulasi. Saya tidak berpikir ada perselisihan tentang apakah probabilitas posterior akan dikalibrasi (dalam arti frekuensi) jika probabilitas sebelumnya, jadi saya curiga simulasi tidak akan meyakinkan siapa pun tentang sesuatu yang baru.

Lagi pula, dalam kasus sampling berurutan yang disebutkan dalam pertanyaan (paragraf ketiga) dapat disimulasikan "sebagaimana adanya" dengan menggambar dari sebelumnya, menggambar sampel yang diberikan ini hingga atau kriteria terminasi lainnya terjadi (kriteria terminasi lain diperlukan karena ada probabilitas positif bahwa probabilitas posterior yang sedang berjalan tidak akan pernah melebihi ). Kemudian, untuk setiap klaim , periksa apakah sampel dasar -parameter positif dan hitung jumlah positif sejati vs positif palsu. Jadi, untuk :$\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$

• Contoh$\mu_i \sim N(\gamma, \tau^2)$
• Untuk : $j=1,\ldots^\ast$
  • Contoh$y_{i,j} \sim N(\mu_i, \sigma^2)$
  • Hitung$p_{i,j} := P(\mu_i>0 \mid y_{i,1:j})$
  • Jika$p_{i,j}>0.95$
    • Jika , penghitung benar positif$\mu_i>0$
    • Jika , increment false positive counter$\mu_i\leq0$
    • Istirahat dari dalam untuk loop
  • $\ast$ beberapa kondisi melanggar lainnya, seperti$j\geq j_{\max}$

Rasio positif sebenarnya dengan semua positif akan setidaknya , yang menunjukkan kalibrasi klaim.$0.95$$P(\mu>0 \mid D)>0.95$

Implementasi Python lambat-dan-kotor (bug sangat mungkin + ada bias menghentikan potensial dalam bahwa saya debugged sampai saya melihat memegang properti kalibrasi yang diharapkan).

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

Saya mendapat untuk proporsi positif sejati untuk semua klaim. Ini lebih dari karena probabilitas posterior tidak akan mencapai tepat . Untuk alasan ini kode melacak juga probabilitas posterior "diklaim" rata-rata, yang saya dapatkan .$0.9807$$0.95$$0.95$$0.9804$

Kita juga dapat mengubah parameter sebelumnya untuk setiap untuk menunjukkan kalibrasi "atas semua kesimpulan" (jika prior dikalibrasi). Di sisi lain, seseorang dapat melakukan pembaruan posterior mulai dari hiperparameter sebelumnya "salah" (berbeda dari apa yang digunakan dalam menggambar parameter ground-truth), dalam hal ini kalibrasi mungkin tidak berlaku.$\gamma,\tau$$i$

Memperluas jawaban terbaik oleh @ juho-kokkala dan menggunakan R di sini adalah hasilnya. Untuk distribusi sebelum rata-rata populasi, saya menggunakan campuran yang sama dari dua normals dengan rata-rata nol, salah satunya sangat skeptis tentang rata-rata besar.

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