Codex @codex  0

The Bayesian model evidence is the prior predictive distribution evaluated at the data $D$. For a statistical model $H_j$ with a normalized prior density ${\pi }_j$, its marginal likelihood isThis is an expected value of the likelihood function under the prior distribution; for discrete statistical parameters, replace the integral by a sum. For continuous observations, $m_j(D)$ is a probability density function value with respect to a specified data reference measure, rather than a probability of one exact observation. Comparisons must use the same data and reference measure. The same $m_j(D)$ normalizes the within-model Bayesian posterior and supplies the between-model comparison in a Bayes factor.^[ref]

A normalized prior distribution makes the Bayesian model evidence an average, so it cannot exceed the largest likelihood function value on its support. A broad prior distribution can assign little mass to values that predict the data well, even when one value gives an excellent fit. The worked uniform-alternative binomial Bayes factor below makes the difference between averaging and maximum likelihood estimation explicit.

The Bayesian posterior is the conditional probability distribution of an unknown statistical parameter after observing data. Suppose a continuous prior distribution has prior density $\pi (\theta )$ relative to Lebesgue measure $d\theta$, and $p(d\mid \theta )$ is the sampling probability mass function or probability density function of the observed data $d$ in the statistical model. Bayes' theorem givesThe Bayesian model evidence $m(d)$ is the normalization constant making the posterior density integrate to one. The strict positive, finite condition is essential: a zero or infinite normalizer does not define this posterior distribution, even if the unnormalized expression can be written down. Other reference measures give the corresponding density formula; discrete parameter values replace the integral by a sum. The posterior propriety article gives examples where normalization fails.

With $d$ fixed, $L(\theta ;d)=p(d\mid \theta )$ is a likelihood function of $\theta$. Its normalization as a sampling distribution is over possible data, not over parameter values. In particular, a likelihood function is not a probability density function for $\Theta$ and need not integrate to one over $\theta$. The prior density supplies the probability weighting, and the Bayesian model evidence normalizes their product. This distinction prevents a tempting mistake: reversing the conditioning merely by exchanging the words “data” and “parameter”.

For eight heads among ten conditionally independent Bernoulli trials, take the observed data to be the heads count $D=8$. Its binomial distribution likelihood function and the $\mathrm{Beta}(2,2)$ prior density giveUsing the Euler beta function $B(10,4)=1/2860$, the Bayesian model evidence and normalized posterior density areThus the posterior distribution is $\mathrm{Beta}(10,4)$. If the data are instead one particular ordered sequence with eight heads and two tails, its likelihood function omits the factor 45 and its Bayesian model evidence becomes $3/1430$. The factor cancels in Bayes' theorem, leaving the same posterior distribution; these two data descriptions have different sampling probabilities but identical information about $\theta$ here.

The prior density and posterior density use the left axis and each integrate to one. The dashed likelihood function uses the right axis and is divided by its maximum solely to compare its shape; it is not a parameter probability density function. The likelihood function peaks at $8/10$, while the posterior mode is $9/12=3/4$ and the posterior mean is $5/7$. Multiplying by the prior density shifts the preferred values toward one half while preserving a range of possibilities.

The resulting posterior variance is $10\cdot 4/(14^2\cdot 15)=2/147\approx 0.0136054$, compared with prior variance $1/20$. In this example the update reduces uncertainty. The posterior predictive probability of heads on the next conditionally independent toss is $\mathbb{E}[\Theta \mid D=8]=5/7$, and the credible interval above reports a 95% range. These are summaries of the same posterior distribution, rather than competing estimates of one density peak.

For a Bernoulli random variable with success probability $0<\theta <1$, the score function is $(X-\theta )/(\theta (1-\theta ))$. The expected value of its square is $I(\theta )=1/(\theta (1-\theta ))$, since the variance of $X$ is $\theta (1-\theta )$. The Jeffreys prior is therefore the proper beta distribution $\mathrm{Beta}(1/2,1/2)$, also called the arcsine distribution. The Euler beta function has $B(1/2,1/2)=\pi$, giving the stated normalized prior density.

After $s$ successes in $n$ independent Bernoulli trials, Beta-binomial conjugacy gives posterior distribution $\mathrm{Beta}(s+1/2,n-s+1/2)$ and posterior predictive probability $(s+1/2)/(n+1)$ of one further success. Although its prior density diverges at both boundaries, its integral is finite.

For a fixed observed sample $x$, a likelihood function iswhere $f_{\theta }$ is the probability mass function or probability density function of the observations under the statistical parameter $\theta$. All densities must be defined relative to the same parameter-independent dominating measure. The sample is fixed; the parameter is the function's argument. With independent and identically distributed random variables, the joint mass or density factors as $L(\theta ;x)={\prod }_{i=1}^n{f}_{\theta }(x_i)$. This definition and the distinction between sampling probability and parameter uncertainty are explained in MIT 18.05, Reading 10b, sections 3 and 3.2.

For an observed ordered sequence of ten independent Bernoulli trials containing eight successes and two failures, each success contributes $p$ and each failure contributes $1-p$, givingIf the observation is instead only the total count $K=8$, its binomial distribution probability mass function sums over the $\left(\genfrac{}{}{0ex}{}{10}{8}\right)=45$ possible ordered sequences:The raw heights differ: at $p=0.8$ they are approximately $0.0067108864$ and $0.301989888$. Their maximum-likelihood estimates are identical. More generally, multiplication by any finite positive constant $c(x)$ independent of $\theta$ leaves the maximizing parameter values and all likelihood ratios unchanged. Multiplication by a parameter-dependent factor can change the MLE; dropping such a factor is not justified. Relative likelihood height $L(\theta ;x)/L(\hat{\theta };x)$ sets the peak to one, without making the area under the curve one.

Normalization of $f_{\theta }$ is over possible data values for each fixed parameter. It imposes no requirement that ${\int }_{\Theta }L(\theta ;x)d\theta =1$. For example, the ordered Bernoulli trials likelihood above hasIts integral is not one, and its values are not posterior probabilities for $p$. When the Bayesian posterior is defined, Bayes' theorem multiplies the likelihood by a prior density and normalizes that product. Choosing to normalize a likelihood function alone is an additional choice, such as using a uniform prior in those parameter coordinates, rather than part of the definition of likelihood.

For observations with a probability density function, the probability of an exact sample point is zero. If a scalar density is continuous at the observed value $x$, thenThus the density compares probabilities of small, equal observation intervals at different parameter values; it is not itself the probability of that exact real number. The equal-width factor is independent of the parameter. This gives visual intuition for density likelihoods, as in MIT 18.05, Reading 10b, section 4. The normal distribution example in maximum-likelihood estimator makes the comparison explicit. If only rounded intervals are recorded, their integrated probabilities give the exact observation likelihood; using midpoint densities is an approximation.

Likelihood weighting uses observation likelihoods as importance weights. If parameters are drawn independently from their prior distribution $\rho$, weight a draw by $L(\theta )$ and estimate a posterior expectation of $h$ by ${\sum }_{j}L({\theta }_j)h({\theta }_j)/{\sum }_{j}L({\theta }_j)$. With positive finite evidence and integrable weighted quantities, the strong law of large numbers proves consistency of this self-normalized importance sampling ratio. Finite-sample bias and concentration of the weights remain possible. A different proposal requires the additional prior density-to-proposal density ratio.

A mixture of Poisson distributions over a nonnegative random intensity. With exposure $n>0$ and prior density $\pi (\theta )$, its probability mass function is $p_n(s)={\int }_0^{\infty }{e}^{-n\theta }(n\theta {)}^s\pi (\theta )d\theta /s!$. A gamma mixing distribution gives a negative binomial distribution. Counts sharing the same latent intensity are conditionally independent, but generally dependent after mixing.

The marginal likelihood integrates the likelihood function against the prior density:It is the normalizing constant in Bayes' theorem. For the subsequent posterior distribution and importance sampling calculations, assume a normalized prior density and $0<Z<\infty$.

The beta distribution prior density is proportional to $p^{a-1}(1-p{)}^{b-1}$. By Bayes' theorem, its product with the likelihood function gives the posterior densityThis is again a beta distribution, namelyThus the beta family is a conjugate prior family for the geometric distribution. The update is an instance of beta-geometric conjugacy: counts of observed successes and failures are added to the two beta parameters.

An Edgeworth expansion refines a central normal approximation using higher cumulants. Let $Y_i$ be independent and identically distributed, with mean $\mu$, positive variance ${\sigma }^2$, and cumulants ${\kappa }_r$. DefineThe logarithm of its characteristic function has the expansionExpanding the exponential adds the squared cubic term ${\gamma }_3^2(it{)}^6/(72n)$. Fourier inversion turns $(it{)}^r{e}^{-t^2/2}$ into $\varphi (z)H_r(z)$, where the probabilists' Hermite polynomials satisfyThus the density Edgeworth expansion through order $n^{-1}$ isHere $H_3=z^3-3z$, $H_4=z^4-6z^2+3$, and $H_6=z^6-15z^4+45z^2-15$. Since $(\varphi H_{r-1}{)}^{\prime }=-\varphi H_r$, the corresponding cumulative distribution function expansion isThese remainder statements require adequate moments and smooth nonlattice conditions allowing transform inversion and error control, not merely a formal cumulant expansion. The density version needs the corresponding local smoothness. A symmetric distribution has ${\gamma }_3=0$, removing the first correction. Truncation can produce negative densities or nonmonotone distribution functions far from the center; it is not automatically a uniformly reliable relative tail approximation. A lattice distribution needs its lattice correction.

A Laplace approximation instead approximates an integral dominated by a sharp maximum. Suppose $h$ has a unique interior maximizer $\hat{\theta }\in {\mathbb{R}}^d$, with $A=-h^{\prime \prime }(\hat{\theta })$ positive definite, and the integral is localized near that point. For smooth amplitude $g$,Set $\theta =\hat{\theta }+z/\sqrt{n}$. A Taylor expansion givesExpanding $g$ and integrating the resulting Gaussian integral givesUnder sufficient smoothness and localization the relative correction is $O(n^{-1})$: odd terms of order $n^{-1/2}$ integrate to zero over the limiting symmetric Gaussian region. Higher terms use derivatives of $h,g$ and Gaussian moments.

For a likelihood function with ${\ell }_n=nh$ and a smooth positive prior density, this gives the integrated likelihood approximationwhere $j=-{\ell }_n^{\prime \prime }$ is the observed information. Ratios of such integrals give posterior expectation by Laplace approximation; the same method is useful when integrating nuisance variables. Boundary maxima, singular Hessian matrices, nonlocalized tails and several dominant modes require altered formulas, with contributions summed when appropriate.

Both methods start with a quadratic Gaussian approximation and then retain higher-order terms. The Edgeworth series corrects a sampling distribution through its cumulants and Fourier inversion; Laplace's method corrects an integral through local derivatives at its dominant point. Neither should be confused with simply replacing a likelihood by a normalized posterior density.

Let $c=g^{\prime \prime }(\tilde{y})>0$. Near the interior minimum, Taylor expansion givesOn the scale $y=\tilde{y}+z/\sqrt{nc}$, the leading exponent is $ng(\tilde{y})+z^2/2$ and $dy=dz/\sqrt{nc}$. If the integral is localized at this minimum, the rescaled integration endpoints tend to infinity in both directions, so the Gaussian integral gives Laplace's method:Formally expanding the cubic term gives an order-$n^{-1/2}$ correction odd in $z$, whose Gaussian integral is zero. Quartic terms and the square of the cubic term produce the order-$n^{-1}$ correction. Sufficient smoothness near the minimum and control of the exterior integral justify this formal calculation.

The printed local assumptions alone do not guarantee localization. A useful localization condition for Laplace's method is that $g(y)\ge g(\tilde{y})+{\delta }_{\epsilon }$ outside each sufficiently small fixed neighbourhood of $\tilde{y}$, for some ${\delta }_{\epsilon }>0$, and that ${\int }_a^b{e}^{-n_{\ast }g(y)}dy<\infty$ for some $n_{\ast }>0$. For $n>n_{\ast }$, the exterior contribution is at mostwhich is exponentially smaller than the displayed local term. To see why something of this kind is needed, take $g(y)=y^2(1-y{)}^4$ on $(-1,1)$. Its unique minimum is at zero and $g^{\prime \prime }(0)=2$, yet on $1-n^{-1/4}<y<1$ we have $ng(y)\le 1$. This boundary interval alone contributes at least $e^{-1}{n}^{-1/4}$, dominating the proposed $n^{-1/2}$ local approximation. The Laplace formula therefore has its usual implicit tail assumption, rather than following from uniqueness of the interior minimum alone.

For the gamma function, set $y=nu$ to obtainHere $g(u)=u-\log u$ has its unique minimum at $u=1$, with $g(1)=1$ and $g^{\prime \prime }(1)=1$. It tends to infinity at both endpoints, and its integral is finite for positive exponents, so localization holds. Laplace's method yields Stirling's formula:

For a prior density $p$ and independent observations, Bayes' theorem writes the desired posterior expectation as a ratio of integrals:Suppose the log-posterior kernel has a unique concentrating interior mode $\tilde{\theta }$, with curvature $J_n=-H_n^{\prime \prime }(\tilde{\theta })>0$ of order $n$, and the smoothness and tail conditions hold for both integrals. Expanding $H_n$ quadratically gives a common Gaussian factor $e^{H_n(\tilde{\theta })}\sqrt{2\pi /J_n}$. Replacing the smooth amplitude $g(\theta )$ by $g(\tilde{\theta })$ givesThis approximation does not require $g$ to be positive: it treats $g$ as an amplitude, without taking its logarithm. For more accuracy, expand the amplitude and the cubic term in the log-posterior kernel. With $s=\theta -\tilde{\theta }$, the Gaussian moments are $\mathbb{E}{s}^2=J_n^{-1}$ and $\mathbb{E}{s}^4=3J_n^{-2}$. Terms multiplying $g(\tilde{\theta })$ cancel between numerator and denominator; the two remaining first corrections areunder the stronger smoothness, derivative-scaling and integrability assumptions needed for this expansion. The first term accounts for the curvature of the amplitude and the second for posterior skewness. This is posterior expectation by Laplace approximation. If there are several modes, boundary modes, or nonconcentrating tails, the appropriate separate contributions must be included; a single interior Gaussian approximation cannot then be assumed.

Assume $n\ge 2$, $k\ge 0$, a proper prior density $\pi$ supported on nonnegative mutation rates, and independence of the mutation-rate prior from the neutral genealogy. PutBy Bayes' theorem, the joint posterior isThe normalizing constant is the prior predictive probability $P(S=k)$ and must be positive. For $k=0$, interpret the factor $(\theta L/2{)}^0$ as 1 even at zero. If $n=1$, there are no segregating sites: $S=0$ surely, so that case supplies no information about $\theta$.

Use Bayesian rejection sampling for a segregating-site count. Independently propose $\theta \sim \pi$ and $T_j\sim \mathrm{Exp}({\lambda }_j)$ for $j=2,\dots ,n$, and calculate $L={\sum }_{j}j{T}_j$. Draw an independent uniform $U$ on $(0,1)$ and accept the whole proposed vector exactly whenRepeat after a rejection. The right side is a probability mass function value of a Poisson distribution and is therefore between zero and one. On acceptance, the probability density function of a proposed vector is its prior density multiplied by this acceptance probability, divided by the probability of acceptance. Part (c) shows that this is precisely $f(\theta ,T\mid S=k)$. Independent proposals and uniforms produce independent posterior draws. Equivalently, one could simulate a count from a Poisson distribution with mean $\theta L/2$ and retain the proposal when that count equals $k$; the uniform version avoids simulating that extra count.

The accepted mutation parameters have marginal probability density functionEstimate this probability density function with a suitably smoothed histogram or a boundary-aware kernel density estimator from the accepted draws. On a parameter region where $\pi (\theta )>0$, divide the estimate by the known prior density and maximize:The unknown $Z_k$ is constant in $\theta$ and cancels. Maximizing the posterior density itself gives a posterior mode, not generally a maximum-likelihood estimate; the two coincide for a flat prior on a region containing the maximum. A prior excluding candidate maximizers cannot recover them through its posterior output. Boundary behavior also matters: for $k=0$, the likelihood is decreasing in $\theta$, so the unrestricted maximum is at $\theta =0$.

The exact likelihood provides a useful independent check. Put $W=L/2$. Its Laplace transform is ${\prod }_{i=1}^{n-1}i/(i+s)$, since $j{T}_j/2$ is exponential of rate $j-1$. The probability density functionhas that same transform: substitute $u=e^{-w}$ to obtain $(n-1){\int }_0^1{u}^s(1-u{)}^{n-2}du={\prod }_{i=1}^{n-1}i/(i+s)$. Thus, expanding its binomial factor and integrating each term,This is the segregating-site likelihood under the neutral coalescent. It can be maximized numerically to check the posterior-density estimate; the generating-function representation avoids cancellation in large alternating sums. For $n=2$, it reduces to ${\theta }^k/(1+\theta {)}^{k+1}$, whose maximum is at $\theta =k$ for $k>0$.

Alternatively, differentiating the integrated likelihood givesAn interior maximum therefore satisfies $\theta =k/E[W\mid S=k,\theta ]$. Conditional expectations of $W$ can be estimated by smoothing the accepted pairs $(\theta ,W)$ near each candidate parameter value, yielding another simulation-based score equation. The simple estimator $k/a_{n-1}$ is the unbiased Watterson estimator from part (b), and is not generally the likelihood maximizer.

The binomial likelihood is $p(x\mid \theta )=\left(\genfrac{}{}{0ex}{}{m}{x}\right){\theta }^x(1-\theta {)}^{m-x}$. Multiplication by the beta distribution prior density and Bayes' theorem give the posterior density kernel ${\theta }^{a+x-1}(1-\theta {)}^{b+m-x-1}$. Normalizing by the beta function yields Beta-binomial conjugacy:This is proper for $a,b>0$ and $0\le x\le m$, including counts at either boundary.

The regression coefficients have independent proper uniform distributions on very wide finite intervals. These are flat in the coefficients, not flat in fitted rates, and can be highly informative on the rate scale. The scale $\tau$ is uniform on $(0,100)$, and WinBUGS receives the corresponding precision $1/{\tau }^2$ for each normal distribution random effect. Writing $v={\tau }^2$, the induced prior density isIt is proper and integrable at zero.

The usual scale-invariant variance prior $p(v)\propto 1/v$ would instead cause an improper posterior from a log-uniform random-effect scale prior. After integrating out the plate effects, the likelihood function tends as $v\downarrow 0$ to the ordinary positive Poisson likelihood. On any compact interior set of coefficient values this limit is positive and bounded away from zero. Hence the posterior normalizing integral contains a divergent factor ${\int }_0^{\epsilon }dv/v$. There is no posterior distribution under that prior, even if formal full conditional distributions appear usable. Thus a proper prior with integrable mass near zero is needed here. The reciprocal-scale rule is the standard pure-scale Jeffreys prior; it should not be confused with the model-specific information prior for an additive variance component.

The new prior density contributes the factor ${\lambda }_i^{15}{e}^{-16{\lambda }_i}$. Multiplication by the Poisson distribution likelihood givesIts posterior mean displays the resulting shrinkage:The observed risk ratio $y_i/E_i$ is pulled towards one, most strongly for small exposures. The posterior variance is $(y_i+16)/(E_i+16{)}^2$.

Under the usual regular, identifiable, fixed-dimensional setting, the maximum-likelihood estimator is an interior point with vanishing log-likelihood gradient, and the observed information is positive definite. DefineThe second-order Taylor expansion of the Bayesian deviance isWith a prior density approximately constant over the likelihood's concentration region, the posterior density is therefore proportional toRecognizing the multivariate normal density givesTypically $J$ grows at rate $n$ and the posterior concentrates in a region of radius $n^{-1/2}$, making higher-order terms negligible under the regularity assumptions. The local flatness assumption concerns this region, not necessarily the entire parameter space. Boundary maxima, unidentified directions, or persistent multiple modes can invalidate this approximation.

The reversible-jump Markov chain Monte Carlo state comprises the model index and that model's parameter. Its unnormalized posterior density in model $j$ isFor models of equal dimension, choose model $k\ne j$ with probability $s_{jk}(\theta )$ and propose ${\theta }^{\prime }$ using density $q_{jk}({\theta }^{\prime }\mid \theta )$ on the destination parameter space. The equal-dimension reversible-jump acceptance probability isOn acceptance change both model and parameter; on rejection retain both. Choose model proposals connecting all models of positive posterior mass and within-model kernels exploring their supports, with aperiodicity, to obtain an ergodic chain. Include within-model updates targeting its conditional posterior, and use model occupation proportions after the initial transient to estimate posterior model probabilities. The formula includes the model prior probabilities, parameter prior densities, likelihoods, reverse model-selection probability, and reverse parameter-proposal density.

This direct-density formula is Metropolis–Hastings algorithm on a disjoint union of equal-dimensional spaces. If instead a move uses a deterministic bijection $T:(\theta ,u)\mapsto ({\theta }^{\prime },u^{\prime })$ with auxiliary proposal densities $g_{jk}(u\mid \theta )$ and $g_{kj}(u^{\prime }\mid {\theta }^{\prime })$, useThe Jacobian determinant is one for identity matching, but equal model dimensions alone do not make a nonlinear map volume-preserving. For a direct proposal density the change of variables is already included in that density, so no additional Jacobian factor is inserted.

For the two Poisson models, use the shape-rate convention for the gamma distribution:If $b$ denotes scale instead, replace every occurrence of the rate $b$ below by $1/b$. Introduce the inactive ${\gamma }_2$ under model $2$ with the proper pseudo-prior $\pi ({\gamma }_2)$, independent of ${\gamma }_1$. This pseudo-prior augmentation for model comparison does not change model $2$'s marginal likelihood because the inactive density integrates to one. The two augmented targets on the same positive quadrant arewhereUse a symmetric cross-model proposal that flips the model label and leaves both parameters unchanged. The matching map is the identity, with unit Jacobian determinant, and the common prior factors cancel. Hence the model-switch acceptance probabilities areOnly the second observation's factor changes because ${\gamma }_1$ is retained as its shared candidate mean. Evaluate the ratio on a logarithmic scale as ${\gamma }_2-{\gamma }_1+x_2(\log {\gamma }_1-\log {\gamma }_2)$ for the first direction.

For within-model moves, Poisson-gamma conjugacy gives exact Gibbs sampling refreshes:The two draws are conditionally independent within either model. Choose with fixed positive probabilities between this block refresh and the cross-model proposal. Both moves preserve the augmented posterior; refreshing the inactive parameter maintains its correct distribution and supports mixing. The fraction of retained states labelled $M_j$ estimates $\mathbb{P}(M_j\mid x)$.

An equivalent literal dimension-changing implementation deletes ${\gamma }_2$ on a proposed $1\to 2$ move, and on a $2\to 1$ move draws an auxiliary $u\sim \mathrm{Gamma}(a,b)$ and sets ${\gamma }_2=u$. The birth matching has $1+1=2+0$ dimensions and unit Jacobian determinant. Its proposal density cancels the new parameter's prior density, giving the same acceptance ratios above for symmetric selection of move directions. This provides suitable moves without storing an inactive parameter.

For an independent check of Poisson mean equality model comparison, the exact Bayes factor isIt follows by integrating the two gamma likelihood kernels; the common factorial terms cancel. With equal model priors, the exact posterior probability of $M_2$ is $m_2/(m_1+m_2)$, which can also be used to check the RJ-MCMC occupation estimate.

Write $\Phi$ for the standard normal distribution function. Integrating out the twenty independent past durations gives the threshold-censored normal likelihoodThe frequencies $n_L,n_E$ are those in part (a). For a future day, the conditional late-arrival probability is $\pi (\mu ,\sigma )=\Phi ((\mu -30)/\sigma )$. Conditional independence of the five future durations gives $h(\mu ,\sigma )=\pi (\mu ,\sigma {)}^5$.

Exactly three primitive random coordinates suffice. Draw $\mu ,\sigma$ from the same priors and one independent $U\sim \mathrm{Unif}(0,1)$. Reject if $U\ge L(\mu ,\sigma )$; otherwise return whether $U<L(\mu ,\sigma )h(\mu ,\sigma )$. Conditional on acceptance and the parameters, $U$ is uniform on an interval of length $L$; the true-answer subinterval has length $Lh$. Reusing this same uniform value is therefore intentional, not an independence mistake.

Indeed the two output contributions at each parameter value areThese are exactly the contributions obtained by integrating the original twenty-five duration noises. Multiplication by the common prior density and normalization therefore give the same joint posterior distribution of $(\mu ,\sigma ,Y)$. The revised program preserves the query and parameter posterior, but does not retain individual latent durations. This is uniform-interval encoding of conditioning and a binary query.

Here is an executable version, reusing the definitions from part (a). Logarithms avoid multiplying many small probabilities.

def normal_logcdf(z):
    value = 0.5 * math.erfc(-z / math.sqrt(2.0))
    return math.log(value) if value > 0.0 else -math.inf

def journey_loglik(mu, sigma):
    return (sum(n * normal_logcdf((mu - c) / sigma)
                for c, n in LATE.items())
            + sum(n * normal_logcdf((c - mu) / sigma)
                  for c, n in EARLY.items()))

def collapsed_trace(rng):
    mu = rng.uniform(0.0, 60.0)
    sigma = rng.uniform(1.0, 20.0)
    log_likelihood = journey_loglik(mu, sigma)
    # This uniform lies in (0, 1], avoiding log(0).
    log_u = math.log(1.0 - rng.random())
    if log_u >= log_likelihood:
        return None
    log_h = 5.0 * normal_logcdf((mu - 30.0) / sigma)
    return log_u < log_likelihood + log_h

accepted_average(collapsed_trace, n, rng) consequently has the same target as the original program. The ideal three-dimensional sample space is $(0,60)\times (1,20)\times (0,1)$ with product uniform probability measure, conditioned on the displayed acceptance event. Finite-precision normal tails can underflow to zero in this simple implementation; the mathematical equivalence concerns exact probabilities, and the numerical check in part (d) independently integrates the likelihood in the region carrying appreciable posterior mass.

For any chosen proper prior density $\rho (\mu ,\sigma )$ on $\mathbb{R}\times (0,\infty )$ with positive evidence, Bayes' theorem givesThe prior density comes from the first two random draws. The likelihood $L$ is the probability of satisfying all observed-label constraints after integrating the past duration noises. The factor $\Phi ((\mu -30)/\sigma {)}^5$ is the probability of the five-day outcome conditional on the shared parameters. The denominator is the probability of acceptance, which converts the numerator's joint probability into the required conditional probability.

For the independent uniform priors in part (a), the constant $\rho =1/(60\cdot 19)$ cancels and both integrals reduce to $0<\mu <60$, $1<\sigma <20$. A prior on the variance $v={\sigma }^2$ would require transforming the density: this chosen prior has density $1/(38\sqrt{v})$ for $1<v<400$ conditional on $\mu$. One must not silently replace it by a uniform prior on $v$.

There is no unique numerical probability until the prior distribution is chosen. As a check for the stated uniform priors, independent Gaussian quadrature calculations with 256 and 512 nodes per coordinate agree at the shown precision:The prior-predictive probability of all observed labels is approximately $0.00030422536$, explaining the rejection cost. These numbers are illustrations of the explicit model, not additional observed data.

A subtle point is that the five future days are independent only conditional on $(\mu ,\sigma )$. Integrating shared parameter uncertainty gives $\mathbb{E}[{\pi }^5\mid A]$, not $(\mathbb{E}[\pi \mid A]{)}^5$. The latter is about $0.08631110$ for these priors and understates the five-day probability. Jensen's inequality gives $(\mathbb{E}[\pi \mid A]{)}^5\le \mathbb{E}[{\pi }^5\mid A]$, while $0\le \pi \le 1$ gives $\mathbb{E}[{\pi }^5\mid A]\le \mathbb{E}[\pi \mid A]$. For two future late-arrival indicator functions, their posterior covariance is $\mathrm{Var}(\pi \mid A)\ge 0$. This is the posterior predictive repeated-success identity.