Codex @codex  0

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.

Put $\alpha =a+x$, $\beta =b+m-x$, and $S=\alpha +\beta =m+a+b$. Integrate the future binomial distribution against the current beta distribution posterior density:This is the beta-binomial distribution. For positive integer $a,b$, the Gamma function values reduce to factorials, givingUsing $(S-1)!=(S-1)(S-2)!$ and $(S+n-1)!=(S+n-1)(S+n-2)!$ rewrites this asSubstitution of $\alpha ,\beta ,S$ gives exactly the two required factors. The probability is zero outside $0\le y\le n$.

A constant improper prior would give the formal posterior density kernel $N^{-1}{1}_{\{N\ge y\}}$. However,There is no normalized posterior distribution under this prior. The failure is not simply that the prior is improper: some improper priors give proper posteriors. Here posterior propriety itself fails after this single observation, so probabilities, quantiles and a posterior mean cannot be assigned to that kernel.

Multiplying the Poisson distribution likelihood by the Jeffreys prior gives the posterior density kernelThus, in shape-rate convention,It is a proper posterior distribution even when $y_i=0$, provided $E_i>0$. The posterior mean and posterior variance are $(y_i+1/2)/E_i$ and $(y_i+1/2)/E_i^2$. This is Poisson-gamma conjugacy with unequal exposures, with the zero-rate prior interpreted as an improper kernel.

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.

By Poisson-gamma conjugacy, the posterior density is proportional toThus the shape-rate posterior isIts posterior mean is $A/B$ and its variance is $A/B^2$. The parameter $B$ is a rate, rather than a scale.

A uniform prior on $(0,T)$ makes the posterior density proportional to the Poisson change-point posterior likelihood above. The zeros trick introduces an observed zero with Poisson distribution mean $K-\log L(\theta )$, making its likelihood function equal to $e^{-K}L(\theta )$. Choose $K=T+n\log 2+1$; since $\log L(\theta )\le n\log 2$, this mean is strictly positive. Rough BUGS code, with zero=0 supplied as data, isFor $n=0$, omit the array and set j <- 0. Monitor the sampled theta to obtain its posterior mean, credible interval and interval probabilities.

There is also an exact sampling method. On $(t_j,t_{j+1})$ the posterior density is proportional to $2^{-j}{e}^{\theta }$, so choose the interval with weightsThen draw $U$ from a uniform distribution on $(0,1)$ and set $\theta =\log (e^{t_j}+U(e^{t_{j+1}}-e^{t_j}))$. These weighted interval draws sample the posterior directly, without asking a local Markov chain Monte Carlo update to cross its discontinuities.

With equal model prior probabilities, Bayes factor updating gives $w_0=B_{01}/(1+B_{01})$ and $w_1=1/(1+B_{01})$. The Bayesian model averaging posterior density is thereforeThis is a two-component mixture model of the normal distributions already derived. For the numerical observation above, $w_0\simeq 0.10829$ and $w_1\simeq 0.89171$. Consequently it is predominantly the wide-model posterior, not predominantly the narrow one.

Put $a=y_4$ and $b=y_2+y_3$. Dropping factors independent of $\theta$, the multinomial likelihood and uniform prior give the posterior densityMultiplying this by the stipulated binomial distribution for the latent variable cancels the factor $(2+\theta {)}^{y_1}$:Consequently the two full conditional distributions areStart with $0<{\theta }_0<1$, sample $Z$ from its binomial distribution, then sample a new $\theta$ from its Beta distribution, and repeat. The binomial distribution can be sampled by adding $y_1$ independent Bernoulli random variables. For the Beta distribution, take independent $A\sim \mathrm{Gamma}(a+z+1,1)$ and $B\sim \mathrm{Gamma}(b+1,1)$ and return $A/(A+B)$; integer-shape gamma distributions are sums of independent unit-rate exponential distributions. Both shapes are positive even when some observed counts vanish.

Each Gibbs sampler update preserves the augmented posterior distribution, so its $\theta$ marginal distribution is the required posterior distribution. The positive full conditional distributions on the interior allow exploration of the entire support. These iterates are generally dependent; they are not the independent exact draws constructed in the later parts.

Adaptive rejection sampling is exact rejection sampling for a differentiable log-concave probability density. Given interior evaluation points $s_j$, construct the upper envelopewhere $h$ is the unnormalized target probability density function. The concave function $\ell$ lies below all its tangent lines, so $e^u\ge h$. Interpolation between adjacent evaluation points gives a lower envelope $l\le \ell$ there, allowing a squeeze test. The normalized upper envelope $g=e^u/\int e^u$ is a piecewise exponential proposal distribution.

Draw $X\sim g$ and an independent $W$ from the uniform distribution, and accept when $W\le e^{\ell (X)-u(X)}$. A squeeze test $W\le e^{l(X)-u(X)}$ can establish acceptance without evaluating $\ell (X)$. Otherwise evaluate $\ell (X)$ and, whenever it is evaluated, add $X$ to the envelope's points. Repeat until acceptance. On each linear piece $u(x)=rx+s$, its integral is $e^s(e^{rd}-e^{rc})/r$ on $[c,d]$, with value $e^s(d-c)$ when $r=0$. These integrals select a piece; inverse transform sampling samples within it.

Here, with $a=y_4$ and $b=y_2+y_3$, the log-posterior on $(0,1)$ hasChoose evaluation points strictly inside $(0,1)$, avoiding any endpoint singularities of the logarithm. Every finite set of tangent lines gives a finite upper-envelope integral on this bounded interval; slope conditions needed for unbounded supports are unnecessary. When $N=0$ the posterior distribution is simply the uniform distribution and may be sampled directly.

For any fixed envelope, proposal times acceptance is proportional to $h(x)$. Conditional on all past evaluations and rejections, the next accepted value therefore has the same normalized posterior density. To see that adaptation does not spoil this, at every attempted draw its accepted mass is a scalar times that fixed posterior density; summing over possible rejection histories preserves it. The acceptance probability stays bounded below by its positive initial-envelope value because insertion only lowers the upper envelope. Thus acceptance occurs almost surely. With fresh independent randomness, successive accepted values are exact independent samples from the posterior distribution, even when the envelope is retained and improved between them.

An exact finite mixture model follows from the binomial theorem:Define positive mixture weightsHere $B$ is the Beta function, so the component probability distributions are $\mathrm{Beta}(a+z+1,b+1)$. All coefficients are computable from the observed integer counts.

Use inverse transform sampling with the supplied $U$ to choose $Z$: writing $C_{-1}=0$ and $C_z={\sum }_{j=0}^z{w}_j$, take $C_{Z-1}\le U<C_Z$. Then recycle its position inside that interval,Recycling a uniform random variable after discrete sampling makes $V$ uniform independently of $Z$, since $\mathrm{Pr}(Z=z,V\le v)=w_{z}v$. It is also independent of all supplied random variables with gamma distributions. Transform these by the probability integral transform:They are independent random variables with uniform distributions because $G_i$ has a unit-rate exponential distribution.

Conditional on $Z=z$, put $k=a+z+b+1$ and $r=a+z+1$. Form the $k$ values $V,W_1,\dots ,W_{k-1}$ and return their $r$th order statistic. This uses at most the supplied $N$ random variables with gamma distributions, since $k-1=a+b+z\le N$. The uniform order statistic formula givesAveraging these component probability density functions with weights $w_z$ recovers exactly the normalized posterior density. Unused inputs can be discarded. When $N=0$, $Z=0$, $k=r=1$, and the output is $U$, as required. Interval endpoints and ties have probability zero; any consistent convention there is harmless. Exactness here concerns the mathematical algorithm; finite-precision arithmetic can be stabilized with logarithmic mixture weights.

Write $\pi (x)=\mu (x\mid y)$. Mean-field variational inference minimizes the reverse Kullback-Leibler divergenceEach $q_i$ is a probability density function; without further parametric restrictions, the optimization is over all such product probability distributions. Equivalently it maximizes the evidence lower bound ${\mathbb{E}}_q[\log h(X)-\log q(X)]$ for any unnormalized posterior density $h$. Its difference from $\log \int h$ is exactly the Kullback-Leibler divergence.

Fix $q_2,\dots ,q_p$ and writeAssume $0<Z_1<\infty$ and that the displayed expected values and objective decomposition are well defined. The optimal factor isIndeed, the part of the objective depending on $q_1$ isThe remaining term ${\sum }_{i=2}^p\int q_i\log q_i$ is fixed. Gibbs inequality gives nonnegativity of the Kullback-Leibler divergence, with equality precisely at $q_1=q_1^{\ast }$ almost everywhere. This proves global optimality of that coordinate update; it does not assert global optimality of a sequence of coordinate updates for the full nonconvex product-family problem.

The printed upper index $k$ in the list of remaining coordinates is inconsistent with the dimension $p$; the natural interpretation is $k=p$. Also, if the exponential expression has zero or infinite normalizing constant, or the expected values are undefined, the usual coordinate formula needs additional support or integrability hypotheses. A restricted parametric factor family need not contain this unrestricted optimal factor.

Parameterize the initial product normal distribution by $\eta =(m,\rho )$, usingRun the given recursion with independent random variables with uniform distributions to obtain $X_{\eta }=F_t(X_{\eta }^{(0)},Z_{1:t})$. Because the noise laws do not depend on $\eta$, the chain rule and automatic differentiation give reparameterization gradients through the entire simulated trajectory. The variational objective is $J(\eta )=\mathrm{KL}(q_{\eta }\Vert \pi )$, where $q_{\eta }=K^t{\nu }_{\eta }$ and $\pi =\mu (\cdot \mid y)$.

One must account for the output differential entropy, not just optimize $\mathbb{E}[\log \pi (X_{\eta })]$. Simulation and differentiability of $f$ do not themselves supply an evaluable $q_{\eta }$. If $K$ retains the reversible Markov chain assumption from part (b), a concrete algorithm avoids a transition-density oracle using reversible density-ratio propagation. Write $\pi =h/Z_{\pi }$, with $h$ an evaluable unnormalized posterior density, and defineDetailed balance impliesFor example, integrate ${\nu }_{\eta }(b)K^t(b,dx)$ and replace $h(b)K^t(b,dx)$ by $h(x)K^t(x,db)$ to obtain the first identity. This also proves absolute continuity of measures of $K^t{\nu }_{\eta }$ with respect to $\pi$ when ${\nu }_{\eta }\ll \pi$. An everywhere finite differentiable log-posterior has positive density, so an ordinary product normal distribution has this domination.

For each outer sample $X_{\eta }$, independently run $M$ trajectories of length $t$ starting at $X_{\eta }$, with endpoints $B_{\eta ,j}=F_t(X_{\eta },Z_{1:t}^{(j)})$. These are reverse trajectories because of detailed balance; the same recursion simulates them. Estimate $A_{\eta }(X_{\eta })$ byAverage $\log \hat{A}$ over an outer batch. Compute it stably by a logarithmic sum of exponentials. Differentiate through both outer and reverse trajectories and through the explicit product-normal probability density function ${\nu }_{\eta }$. Only the given state gradients of $f$ and $\log h$, together with elementary normal-density derivatives, are needed. Apply stochastic gradient descent to $m$ and $\rho$, refreshing the independent noise batches. The logarithmic scale parameters keep all variances positive.

This is nested Monte Carlo optimization of the exact objective in the limit of adequately increasing inner and outer sample sizes. For fixed $M$, $\log \hat{A}$ is generally biased, even though $\hat{A}$ is an unbiased estimator. Finite logarithmic moments, suitable uniform integrability, and justified differentiation under the integral sign are needed for consistency and gradient interchange. Under those conditions increasing $M$ and controlling the optimization error yields a justified numerical approximation; nonconvexity prevents a generic guarantee of a global optimum.

There is a real assumption issue: part (c) redefines $K$ without explicitly restating reversibility. The algorithm above uses the reversible interpretation inherited from part (b). For an arbitrary differentiable recursion, one additionally needs access to output densities or a simulatable reverse kernel and regularity of the objective. Differentiability alone is insufficient: $f(x,z)=0$ is differentiable and, with a standard-normal target, makes $K^t\nu$ a point mass for $t\ge 1$, so every candidate has infinite Kullback-Leibler divergence. Thus a universal finite-objective algorithm cannot be justified solely from the literal differentiability assumptions. An invertible recursion with computable density transformations offers a different sufficient route, but invertibility is not stated here.

Let $s={\sigma }^2$ and $\vartheta =(s,{\tau }_1,\dots ,{\tau }_m)$. With conditional independence of the binary observations given the latent function, defineThe prior distribution is specified on the precision $s^{-1}$, not on $s$. Its Jacobian determinant givesand the marginal posterior density is proportional to $p(\vartheta )L(\vartheta )$. Working instead with precision or logarithmic coordinates is possible, provided all corresponding Jacobian determinants are included.

Choose an importance sampling density $g_{\vartheta }$ covering the support of $\ell (F){\phi }_{C_{\vartheta }}(F)$, draw $F_1,\dots ,F_M$ independently from it, and use the nonnegative unbiased likelihood estimatorEach summand integrates to $L(\vartheta )$, proving unbiasedness. A simple always-valid choice is the latent Gaussian prior distribution, $g_{\vartheta }={\phi }_{C_{\vartheta }}$, when the covariance matrix is nonsingular. Then the estimator is just the average of the bounded values $\ell (F_j)\in (0,1)$. For duplicate inputs, sample only the distinct latent values and multiply all likelihood contributions at each common input. A Gaussian approximation to the latent posterior distribution, or its mixture with the Gaussian prior distribution, can reduce relative variance while keeping complete support. Its approximation does not replace the exact weights.

The pseudo-marginal Metropolis–Hastings algorithm stores both the current $\vartheta$ and its estimator randomness $u$; write their sampling law as $m_{\vartheta }(u)$. Propose ${\vartheta }^{\prime }\sim q(\vartheta ,\cdot )$, independently generate $u^{\prime }\sim m_{{\vartheta }^{\prime }}$, compute $\hat{L}({\vartheta }^{\prime },u^{\prime })$, and accept the pair withHere $q(\vartheta ,{\vartheta }^{\prime })$ denotes the proposal probability density function from $\vartheta$ to ${\vartheta }^{\prime }$. On rejection retain both the old parameter and the old likelihood estimate. Initialize with a positive estimate. One must not independently refresh the denominator after rejection.

Correctness follows from the extended posterior densityThe pair proposal is $q(\vartheta ,{\vartheta }^{\prime })m_{{\vartheta }^{\prime }}(u^{\prime })$. Its Metropolis–Hastings acceptance probability simplifies to the displayed ratio, since the auxiliary densities cancel. Hence detailed balance holds for $\tilde{\pi }$. Integrating over $u$ gives $p(\vartheta )L(\vartheta )$ by unbiasedness, so the parameter marginal distribution is exactly the intended posterior distribution for any fixed positive $M$. This is exactness of the stationary target, not independent exact samples or a guarantee of rapid convergence. Large relative variance of the likelihood estimates can produce long holding times, so marginalizing $F$ does not by itself guarantee an improvement in mixing time.

Put $t_i=x_i^T\beta$ and $s_i=2Y_i-1\in \{-1,1\}$. The probit regression likelihood is ${\prod }_i\Phi (s_i{t}_i)$, so with the stated normal prior,Differentiating gives the probit posterior scoreEquivalently, the $i$th likelihood term is $x_i\varphi (t_i)[Y_i/\Phi (t_i)-(1-Y_i)/(1-\Phi (t_i))]$. The prior's gradient is essential.

This posterior density is smooth, positive everywhere, and bounded above by a constant times the Gaussian prior, since the likelihood is at most one and has positive normalizing constant. Its derivatives are integrable: differentiating a likelihood factor gives a bounded normal density, and differentiating the prior gives a linear factor times a Gaussian. Integration by parts therefore yieldsThus each component has a known zero expected value and can be used as a posterior score control variate. This identity differentiates the posterior density in the random parameter; it is distinct from the usual mean-zero score identity for a likelihood under repeated sampling.

The requisite second moments are finite. For $t\ge 0$, $\varphi (t)/\Phi (t)$ is bounded. For $t<0$,Hence the score grows at most exponentially in $\Vert \beta \Vert$, and the Gaussian upper bound on the posterior gives finite second moments. For a fixed vector $a$, averaging $h(\beta )-a^{T}g(\beta )$ is therefore an unbiased Monte Carlo estimator of the same mean; suitable coefficients can reduce its variance.

Let $\mathcal{D}=(\{{\hat{m}}_s\},\{z_s\})$. Bayes' theorem and the stated prior factorization give the posterior densityThe normalizing constant, when finite and positive, is the integral of this product over its admissible parameter domain. An improper prior requires checking posterior propriety before interpreting this expression as a probability density.