Setting up simulation algorithm to check calibration of Bayesian posterior probabilities Figuring out how to simulate something is often the best way to understand the underlying principles.  I am a bit at a loss on exactly how to simulate the following.
Suppose that $Y \sim N(\mu, \sigma^{2})$ and that $\mu$ has a prior distribution that is $N(\gamma, \tau^{2})$.  Based on a sample of $n$ observations $Y_{1}, \dots, Y_{n}$ abbreviated by just $Y$, I am interested in showing to a non-Bayesian that the posterior probability that $\mu > 0 | Y$ is well calibrated, e.g., Prob$(\mu > 0 | P) = P$ where $P$ is the posterior probability.  A related discussion is here
What I really want to show is that if one were to do sequential testing and stop sampling when the posterior probability exceeds some level such as 0.95 the probability that $\mu > 0$ is not $< 0.95$.
I am trying to convince frequentists that Bayesian probabilities are meaningful without getting into any discussion about type I error.  I suppose there is a philosophical problem when talking to a frequentist who entertains null hypotheses in that if the prior is continuous (as above) the probability that $\mu = 0$ is zero and simulations are not needed.  I would appreciate some suggestions about how to think of the whole problem and how to design demonstration simulations.  I am used to doing frequentist simulations where $\mu$ is just set to a single constant; Bayesians don't condition on $\mu$.
For the sequential situation we set a maximum possible sample size, e.g., $n=1000$.
There is a subtlty to the problem that I always have trouble thinking about.  A real skeptic is sometimes worried about a false claim of effectiveness ($\mu > 0$) when the process really has exactly no effect ($\mu=0$).  The subtlty is that the skeptic is "singling out" zero as a special value, and perhaps is giving non-zero probability to the event $\mu = 0$ (?).  Our method of showing that the posteriors is calibrated may not make such a skeptic happy because the skeptic really seems to want to condition on $\mu = 0$ and as Bayesians we only condition on what is knowable.  Perhaps this is a case where the prior distribution that the statistician is using conflicts with a discontinuous prior distribution the skeptic is using?
 A: Simulation results will depend on how the parameter is sampled in the simulation. I don't think there is any dispute over whether the posterior probabilities will be calibrated (in the frequency sense) if the prior probabilities are, so I suspect a simulation will not convince anyone of anything new. 
Anyway, in the sequential sampling case mentioned in the question (third paragraph) can be simulated "as is" by drawing $\mu$ from the prior, drawing samples given this $\mu$ until $p(\mu>0\mid \textrm{samples})>0.95$ or some other termination criterion occurs (another termination criterion is needed since there is positive probability that the running posterior probability will never exceed $0.95$). Then, for every $p(\mu>0\mid \textrm{samples})>0.95$ claim, check whether the underlying sampled $\mu$-parameter is positive and count the number of true positives vs. false positives. So, for $i=1,2,\ldots$:


*

*Sample $\mu_i \sim N(\gamma, \tau^2)$

*For $j=1,\ldots^\ast$:


*

*Sample $y_{i,j} \sim N(\mu_i, \sigma^2)$

*Compute $p_{i,j} := P(\mu_i>0 \mid y_{i,1:j})$

*If $p_{i,j}>0.95$


*

*If $\mu_i>0$, increment true positive counter

*If $\mu_i\leq0$, increment false positive counter

*Break from the inner for loop


*$\ast$ some other breaking condition, such as $j\geq j_{\max}$



The ratio of true positives to all positives will be at least $0.95$, which demonstrates calibration of the $P(\mu>0 \mid D)>0.95$ claims. 
A slow-and-dirty Python implementation (bugs very possible + there is a potential stopping bias in that I debugged until I saw the expected calibration property holding).
# (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))

I got $0.9807$ for the proportion of true positives to all claims. This is over $0.95$ as the posterior probability will not hit exactly $0.95$. For this reason the code tracks also the mean "claimed" posterior probability, for which I got $0.9804$.
One could also change the prior parameters $\gamma,\tau$ for every $i$ to demonstrate a calibration "over all inferences" (if the priors are calibrated). On the other hand, one could perform the posterior updates starting from "wrong" prior hyperparameters (different than what are used in drawing the ground-truth parameter), in which case the calibration might not hold.
