Einrichten eines Simulationsalgorithmus zur Überprüfung der Kalibrierung der Bayes'schen posterioren Wahrscheinlichkeiten

Herauszufinden, wie man etwas simuliert, ist oft der beste Weg, um die zugrunde liegenden Prinzipien zu verstehen. Ich weiß nicht genau, wie ich Folgendes simulieren soll.

Angenommen, und hat eine vorherige Verteilung, die . Basierend auf einer Stichprobe von Beobachtungen abgekürzt mit nur , möchte ich einem Nicht-Bayesianer zeigen, dass die hintere Wahrscheinlichkeit ist gut kalibriert, z. B. Prob wobei die hintere Wahrscheinlichkeit ist. Eine verwandte Diskussion ist hierμ 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$

Was ich wirklich zeigen möchte, ist, dass wenn man sequentielle Tests durchführt und die Probenahme stoppt, wenn die hintere Wahrscheinlichkeit ein Niveau wie 0,95 überschreitet, die Wahrscheinlichkeit, dass nicht .< $\mu > 0$$< 0.95$

Ich versuche, Frequentisten davon zu überzeugen, dass Bayes'sche Wahrscheinlichkeiten sinnvoll sind, ohne über Typ-I-Fehler zu diskutieren. Ich nehme an, es gibt ein philosophisches Problem, wenn man mit einem Frequentisten spricht, der Nullhypothesen dahingehend unterhält, dass bei kontinuierlichem Prior (wie oben) die Wahrscheinlichkeit, dass Null ist und keine Simulationen erforderlich sind. Ich würde mich über einige Vorschläge freuen, wie man über das gesamte Problem nachdenkt und wie man Demonstrationssimulationen entwirft. Ich bin es gewohnt, häufig auftretende Simulationen durchzuführen, bei denen nur auf eine einzelne Konstante gesetzt ist. Bayesianer konditionieren nicht auf .μ $\mu = 0$$\mu$$\mu$

Für die sequentielle Situation legen wir eine maximal mögliche Stichprobengröße fest, z. B. .$n=1000$

Das Problem ist subtil, und ich habe immer Probleme, darüber nachzudenken. Ein echter Skeptiker ist manchmal besorgt über eine falsche Behauptung der Wirksamkeit ( ), wenn der Prozess wirklich genau keine Wirkung hat ( ). Die Subtilität ist, dass der Skeptiker Null als speziellen Wert "herausgreift" und dem Ereignis (?) Möglicherweise eine Wahrscheinlichkeit ungleich Null gibt . Unsere Methode zu zeigen, dass die Posterioren kalibriert sind, macht einen solchen Skeptiker möglicherweise nicht glücklich, weil der Skeptiker wirklich auf konditionieren zu wollen scheint und wir als Bayesianer nur auf das bedingen, was erkennbar ist. Vielleicht ist dies ein Fall, in dem die vorherige Verteilung, die der Statistiker verwendet, mit einer diskontinuierlichen vorherigen Verteilung, die der Skeptiker verwendet, in Konflikt steht?μ = 0 μ = 0 μ = $\mu > 0$$\mu=0$$\mu = 0$$\mu = 0$

Die Simulationsergebnisse hängen davon ab, wie der Parameter in der Simulation abgetastet wird. Ich glaube nicht, dass es Streit darüber gibt, ob die posterioren Wahrscheinlichkeiten (im Frequenzsinn) kalibriert werden, wenn die vorherigen Wahrscheinlichkeiten vorliegen, daher vermute ich, dass eine Simulation niemanden von etwas Neuem überzeugen wird.

$\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$

• $\mu_i \sim N(\gamma, \tau^2)$
• $j=1,\ldots^\ast$
  • $y_{i,j} \sim N(\mu_i, \sigma^2)$
  • $p_{i,j} := P(\mu_i>0 \mid y_{i,1:j})$
  • $p_{i,j}>0.95$
    • $\mu_i>0$
    • $\mu_i\leq0$
    • Brechen Sie von der inneren for-Schleife ab
  • $\ast$$j\geq j_{\max}$

Das Verhältnis von echten Positiven zu allen Positiven beträgt mindestens , was die Kalibrierung der Ansprüche .P ( μ > 0 ≤ D ) > $0.95$$P(\mu>0 \mid D)>0.95$

Eine langsame und schmutzige Python-Implementierung (Fehler sehr wahrscheinlich + es besteht eine potenzielle Stopp-Verzerrung darin, dass ich debuggt habe, bis ich die erwartete Kalibrierungseigenschaft gesehen habe).

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

Ich habe für den Anteil der echten Positiven an allen Ansprüchen erhalten. Dies ist über da die hintere Wahrscheinlichkeit nicht genau . Aus diesem Grund verfolgt der Code auch die mittlere "beanspruchte" hintere Wahrscheinlichkeit, für die ich .0,95 0,95 $0.9807$$0.95$$0.95$$0.9804$

Man könnte auch die vorherigen Parameter für jedes ändern , um eine Kalibrierung "über alle Schlussfolgerungen" zu demonstrieren (wenn die Prioritäten kalibriert sind). Andererseits könnte man die posterioren Aktualisierungen ausgehend von "falschen" vorherigen Hyperparametern durchführen (anders als beim Zeichnen des Grundwahrheitsparameters verwendet). In diesem Fall könnte die Kalibrierung nicht zutreffen.$\gamma,\tau$$i$

Die ausgezeichnete Antwort von @ juho-kokkala zu erweitern und R hier zu verwenden, sind die Ergebnisse. Für eine vorherige Verteilung für den Populationsmittelwert mu habe ich eine gleiche Mischung aus zwei Normalen mit dem Mittelwert Null verwendet, von denen einer sehr skeptisch gegenüber großen Mittelwerten ist.

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