Past exam of the mathematics course of the University of Cambridge / 2005 / iii / Paper 47 (Time Series and Monte Carlo Inference)

A real weakly stationary process has finite second moments, a time-independent mean $\mu$, and autocovariance function ${\gamma }_k=\mathrm{Cov}(X_t,X_{t+k})$ depending only on the lag. Use the angular-frequency conventionThe time-series spectral density is nonnegative and integrates to ${\gamma }_0$. The displayed Fourier series is ordinarily justified by ${\sum }_k|{\gamma }_k|<\infty$. More generally, weak stationarity gives a spectral measure of a stationary time series, not necessarily a spectral density. If a spectral density exists without this summability, its Fourier inversion can be understood using Fejér sums, which converge in $L^1$ and at its Lebesgue points. Thus spectral density existence should not be inferred from weak stationarity alone.

A white noise process $WN(0,{\sigma }^2)$ has mean zero and $\mathrm{Cov}({\epsilon }_t,{\epsilon }_s)={\sigma }^2{1}_{\{t=s\}}$. Neither independence nor a normal distribution is needed for this second-order definition. Its autocovariance function is zero at nonzero lags, soFor a real absolutely summable filter $(c_r)$, the series defining $Y_t$ has mean-square convergence: its tails have $L^2$ norm at most $\Vert X_0{\Vert }_2{\sum }_{\text{tail}}|c_r|$. Its mean is $\mu {\sum }_r{c}_r$, andThe sum is absolutely convergent, since $|{\gamma }_X(j)|\le {\gamma }_X(0)$ by the Cauchy-Schwarz inequality. It depends only on $k$, proving that the output is a weakly stationary process. Substitute the spectral integral for ${\gamma }_X$ and interchange the sums and integral, justified by $f_X\in L^1$ and ${\sum }_{r,s}|c_r{c}_s|<\infty$. This givesBecause the coefficients are real, the two factors are complex conjugates. Hence the spectral density transformation under a linear filter isFor complex coefficients the same formula holds with the conjugate-covariance convention.

Applying this result to the order-two moving-average model yieldsEquivalently, its nonzero autocovariances are ${\gamma }_0={\sigma }^2(1+{\theta }_1^2+{\theta }_2^2)$, ${\gamma }_1={\sigma }^2{\theta }_1(1+{\theta }_2)$ and ${\gamma }_2={\sigma }^2{\theta }_2$, together with ${\gamma }_{-k}={\gamma }_k$.

For the order-one moving-average process, invert its filter by a geometric series. The required coefficients areThey are absolutely summable for $|\theta |<1$. Direct cancellation giveswith mean-square convergence. Thus the inverse filter recovers exactly the original white noise process, including its variance; it is an invertible time-series representation.

Let $\Phi (z)=1-{\varphi }_{1}z-\cdots -{\varphi }_p{z}^p$. The causality root criterion for an autoregressive model,is a sufficient condition for a weakly stationary process solution, and gives the causal solution $X_t={\sum }_{j\ge 0}{\psi }_j{\epsilon }_{t-j}$. Indeed $1/\Phi$ is analytic in a disk of radius greater than one, so its Taylor coefficients ${\psi }_j$ decay geometrically. The series therefore has mean-square convergence and defines a linear filter of a stationary time series; multiplying by $\Phi$ verifies the equation.

For the specified order-two model, put $\phi =(1+\sqrt{5})/2$. Since $\phi -{\phi }^{-1}=1$,When $\alpha \ne 0$, its roots are $1/(\alpha \phi )$ and $-\phi /\alpha$. Both have modulus greater than one exactly when $|\alpha |<1/\phi$. The case $\alpha =0$ is white noise and satisfies the same inequality. ConsequentlyThis is the causal stationarity condition selected above. Without causality, a two-sided autoregressive model can also have a noncausal stationary autoregression when its polynomial has no unit-circle zeros. Thus the inequality is not a necessary condition for all conceivable weak stationarity solutions; the causal interpretation is the one required by the subsequent innovations and forecasting questions.

To find the Wold representation, expand $1/\Phi (z)={\sum }_{j\ge 0}{\psi }_j{z}^j$. Matching coefficients gives ${\psi }_0=1$, ${\psi }_1=\alpha$ and ${\psi }_j=\alpha {\psi }_{j-1}+{\alpha }^2{\psi }_{j-2}$. HenceHere $F_1=F_2=1$ are the Fibonacci numbers. The coefficient formula also follows by partial fractions in the displayed factorization. Absolute summability follows from $|\alpha \phi |<1$ and $|\alpha /\phi |<1$. The past observation span equals the past noise span, because the causal series gives one inclusion and ${\epsilon }_t=X_t-\alpha X_{t-1}-{\alpha }^2{X}_{t-2}$ gives the reverse inclusion. Thus these noises are the linear innovations, and there is no additional deterministic component. This is the Fibonacci-coefficient autoregression.

The best linear prediction from an infinite past is the orthogonal projection onto the closed span of $X_T,X_{T-1},\dots$. Future white noise is orthogonal to that span. Substitution in the recursion therefore givesThe respective errors are ${\epsilon }_{T+1}$ and $\alpha {\epsilon }_{T+1}+{\epsilon }_{T+2}$, so their mean squared errors are ${\sigma }^2$ and $(1+{\alpha }^2){\sigma }^2$. No normal distribution assumption is needed for these linear forecasts.

A linear congruential generator starts with an integer seed $S_0\in \{0,\dots ,M-1\}$ and iteratesThe outputs lie in $[0,1)$ and are used as pseudorandom numbers approximating a uniform distribution. They are deterministic once the seed is fixed. The largest possible period is $M$, since there are only $M$ states; good period alone does not establish independence or good multivariate uniformity.

The full-period criterion for a mixed congruential generator says that period $M$ for every seed occurs exactly whenFor $M=12^k=2^{2k}{3}^k$, these conditions reduce toThe multiplier must satisfy both $a\equiv 1(\mathrm{mod}3)$ and $a\equiv 1(\mathrm{mod}4)$. A positive shift alone is insufficient: for example, a shift divisible by three cannot give period $M$ when the multiplier is one modulo three. If the required coprime shift has already been chosen, the condition on the multiplier alone is $a=1+12r$ for an integer $r$.

Use the ratio-of-uniforms method onWriting $x=v/u$ gives $u\le e^{-x^2/(4{\sigma }^2)}$ and $v=ux$. Consequentlybecause the derivative vanishes at $x=\sigma \sqrt{2}$. Thus the normal ratio-of-uniforms envelope is the rectangle $0<u<1$, $|v|<\sigma \sqrt{2/e}$.

Generate independent $U_1,U_2$ with uniform distribution on $(0,1)$, set $u=U_1$ and $v=(2U_2-1)\sigma \sqrt{2/e}$, and accept precisely whenReject and repeat otherwise. On acceptance return $X=v/u$. One can compare logarithms instead of exponentials to avoid underflow.

For a proof, the change of variables $(u,x)\mapsto (u,v=ux)$ has Jacobian determinant $u$. A uniform point in $\mathcal{R}$ therefore gives an $x$-density proportional toNormalizing yields the required normal distribution with mean zero and variance ${\sigma }^2$. The region area is $\frac{1}{2}\int h(x)dx=\sigma \sqrt{\pi /2}$; dividing it by the rectangle area $2\sigma \sqrt{2/e}$ gives acceptance probability $\sqrt{\pi e}/4$.

For $\alpha ,\beta >0$, integrating the Weibull distribution probability density function givesThe inverse transform sampling method takes $U$ uniform on $(0,1)$ and solves $G(X)=U$. HenceSince $1-U$ is also uniform, $\{-\beta \log U{\}}^{1/\alpha }$ is an equivalent implementation. For any $x\ge 0$, the resulting event $X\le x$ has probability $G(x)$, verifying its distribution.

Assume $0\le p\le 1$. Draw an independent component indicator with a Bernoulli distribution of success probability $p$. On success, generate the normal distribution from part (i) and return its absolute value; this has the half-normal distribution probability density function $2f(x)$ on $x\ge 0$. On failure, generate the Weibull distribution from part (ii). Use fresh independent uniforms for the chosen component.

The law of total probability gives the output probability density functionas required. The factor two belongs inside the half-normal component; the probability of choosing that component is $p$, not $2p$.

For independent draws $X_i$ from a normalized proposal $g$, use the importance sampling identityThus the unbiased estimator and its variance areThe probability support condition is $g>0$ wherever $\theta f\ne 0$, and the variance formula requires a finite second moment. Since both probability densities are normalized, no self-normalization of weights is needed.

In the proposed example, normalize $g$ to $g(x)=3x^2$ on $(0,1)$. It is sampled by $X=U^{1/3}$ for uniform $U$. Taking $\theta (x)=x$ givesAs a check, ${\mu }_g=\frac{4}{\pi }{\int }_0^{1}x\sqrt{1-x^2}dx=4/(3\pi )$. The second moment of one weighted observation isso $\mathrm{Var}({\hat{\mu }}_g)=16/(9{\pi }^{2}n)$. The weight's singularity at zero does not make this variance infinite.

To prove the optimal proposal rule, apply the Cauchy-Schwarz inequality:Equality holds when the two factors are proportional almost everywhere, giving the minimum-variance importance distributionFor the present nonnegative integrand, normalization givesUnder this proposal, $xf(x)/g_0(x)=4/(3\pi )$ is constant, so the variance is exactly zero. If sampling is wanted, its distribution function is $1-(1-x^2{)}^{3/2}$, with inverse $X=\sqrt{1-(1-U{)}^{2/3}}$. In general the optimal proposal's normalizing constant can itself be the unknown integral, limiting this zero-variance construction's practical use.

The plug-in estimate isFor a concrete bootstrap construction, suppose the rows $(x_i,y_i)$ are independent and identically distributed, $\mathbb{E}(e_i\mid x_i)=0$, the design moment matrix is nonsingular, and the required moments exist. Use a paired bootstrap: sample $n$ rows with replacement, form $X^{\ast },Y^{\ast }$, refit by ordinary least squares, and compute ${\hat{\theta }}^{\ast }={\hat{\beta }}_1^{\ast }{\hat{\beta }}_2^{\ast }$. Repeat many times. Singular resampled designs require a suitable rank-preserving implementation or replacement replicate; with a regular continuous design their frequency vanishes asymptotically.

The bootstrap-t confidence interval method also needs a standard error in the original sample and every replicate. One asymptotically valid heteroscedasticity-robust covariance estimate isThe delta method for the coefficient product givesRecompute ${\hat{V}}^{\ast },d^{\ast },{\hat{s}}^{\ast }$ in each paired replicate. This is bootstrap inference for a product of regression coefficients.

With fixed rather than random design, resample centered residuals if errors are identically distributed with homoskedasticity, or use a properly specified wild bootstrap for heteroscedastic errors, keeping $X$ fixed and refitting after generating $Y^{\ast }$. The printed assumption $\mathbb{E}{e}_i=0$ alone does not establish independent observations or validity of a bootstrap confidence interval. The formulas below give asymptotic 95% coverage under the stated regularity conditions, not an exact finite-sample guarantee. In particular the ordinary first-order studentization needs $({\beta }_1,{\beta }_2)\ne (0,0)$ and a nonzero asymptotic variance; when both coefficients vanish the product has nonstandard behavior.

Let $q_{.025}^{\ast }$ and $q_{.975}^{\ast }$ be the empirical quantiles of the paired-bootstrap values ${\hat{\theta }}^{\ast }$. The percentile bootstrap confidence interval isIts endpoints use the bootstrap distribution of the estimate itself. The regularity assumptions and resampling scheme are those in the parent solution.

For each bootstrap replicate form the studentized pivotLet $t_{.025}^{\ast },t_{.975}^{\ast }$ be its empirical quantiles. The bootstrap-t confidence interval approximatesSolving both inequalities for $\theta$ reverses the quantile order, givingUsing the same fixed standard error in every bootstrap numerator would not implement this studentized procedure; each replicate needs its own ${\hat{s}}^{\ast }$.

The Metropolis–Hastings algorithm is useful when a target posterior can be evaluated up to a normalizing constant but direct independent sampling or simple conditional sampling is unavailable. It constructs a Markov chain with that posterior as invariant law. After adequate convergence, averages of functions of the chain estimate posterior expectations, and the retained draws approximate marginal distributions, quantiles and credible intervals.

For the random-walk Metropolis algorithm, at current $\theta$ propose ${\theta }^{\prime }=\theta +Z$, where $Z$ has a symmetric distribution, often $N(0,\Sigma )$. With a symmetric proposal the acceptance probability isAccept the proposal when an independent uniform is below this probability; otherwise retain the current value. For a nonsymmetric proposal $q({\theta }^{\prime }\mid \theta )$, include $q(\theta \mid {\theta }^{\prime })/q({\theta }^{\prime }\mid \theta )$ in the ratio. Detailed balance follows because the accepted flow is the minimum of the two proposed directional flows; rejection supplies the remaining self-transition probability. Irreducibility, appropriate recurrence and aperiodicity are also needed for convergence.

In the displayed trace plots, the first parameter moves in small increments with pronounced long-term wandering and a shift from its early region to a later region. This is consistent with an undersized proposal or a strongly correlated slow direction. Measure its acceptance rate: if it is very high, increase that coordinate's proposal variance. The second parameter moves rapidly within regions but changes its level around iteration 4500; that is evidence against treating the full run as a settled stationary sample, and it suggests investigating dependence with the first parameter or separated posterior regions. Its local scale may already be reasonable, so a blanket reduction of all proposal variances would be unwarranted. The third parameter has conspicuous long plateaus, indicating many repeated states: if the recorded acceptance is low, reduce its proposal variance; also check constraints and dependence before deciding that scale alone is the cause.

Rerun with coordinate-specific scale changes, inspect joint dependence, consider a covariance-based block proposal, and run longer chains from dispersed starting points. The MCMC trace diagnosis and proposal tuning should include actual acceptance rates and effective sample sizes; the three marginal trace plots do not by themselves identify correlations or prove convergence. Further practical details are in the three subparts.

Convergence means that the chain's distribution approaches the intended posterior, independently of its initial state, under the relevant ergodicity assumptions. It is separate from obtaining enough effectively independent draws for a precise answer. Start several chains at dispersed plausible points, examine trace plots and stability of posterior summaries, compare between-chain and within-chain variation, and inspect autocorrelation and effective sample size of a Markov chain. Discard an initial transient when appropriate and report Monte Carlo uncertainty.

No finite diagnostic proves convergence. Agreement among chains started in the same basin can miss another mode. If different starting points give incompatible summaries, extend runs or improve proposals rather than merely deleting more observations. Extra thinning does not repair poor exploration; retaining all post-transient draws usually preserves more information for estimating uncertainty.

Single-component updates propose one coordinate while holding the others fixed, with acceptance based on the resulting target ratio. They are cheap and permit individual scale tuning, but can crawl along a narrow correlated ridge. Block updates propose several coordinates jointly and can follow posterior covariance directions, potentially improving mixing substantially.

A practical compromise uses small correlated blocks, or a full Gaussian block with covariance estimated in a pilot run and rescaled for acceptance. Blocking independent coordinates need not help, and a badly scaled high-dimensional block may be rejected frequently. Compare computing cost per effective draw, rather than acceptance rate alone. A Metropolis-within-Gibbs sweep of coordinate updates and a properly accepted block update both preserve the target.

Small random-walk proposal variance gives many accepted tiny moves and strong autocorrelation; excessive variance proposes mostly low-density points, giving low acceptance and long rejection plateaus. Tune each coordinate or block during a pilot run to balance movement and acceptance, and inspect effective sample size per unit computation. Rescale or transform parameters with very different units, and account for posterior correlation rather than adjusting only diagonal variances.

There is no universal optimal acceptance percentage for every target or update dimension. Keep the chosen proposal fixed in the production run, or use an adaptive method with a proved diminishing-adaptation and containment justification. The direction of the trace-based changes for this particular run is given in the parent solution.

For each observation time, use only the currently susceptible risk set $S_t$, with new-infection indicators $Y_{it}$. PutThis is the spatial Bernoulli infection model, under conditional independence of the indicators. Distances must be positive. Already infected trees should not be included as new susceptible trials. The minus sign in the power kernel is present in the original PDF, although it is absent in the converted TeX.

The parameter $\alpha$ must be nonnegative: for $\alpha <0$ and positive exposure, $1-e^{-\alpha A_{it}}<0$ is not a probability. Therefore the stated unrestricted normal prior distribution cannot literally define the model. A coherent interpretation is a normal prior distribution restricted to $\alpha >0$,Here $\Phi$ is the standard normal distribution function. Use the same proper prior distribution in both kernel models. The given Gaussian prior distributions for $\beta$ and $\gamma$ are mathematically valid on the real line for finite positive distances, though restricting them to positive values would impose decaying kernels and would require the corresponding normalized prior distributions.

The full conditional kernels areNeither is generally normal or another elementary conjugate family: the infection factors contain $1-e^{-\alpha A}$, and $A$ depends nonlinearly on $\beta$. Thus a simple conjugate Gibbs sampler is unavailable. This does not make Gibbs sampling logically impossible; exact rejection or slice conditional draws could be used if implemented. A straightforward alternative is Metropolis–Hastings algorithm updates within each coordinate.

For example, put $\eta =\log \alpha$ and use symmetric Gaussian random-walk proposals for $(\eta ,\beta )$. The target in these coordinates isThe factor $e^{\eta }$ is the Jacobian determinant and must be retained. Accept according to the ratio of this target, either for a block or one coordinate at a time. Posterior averages and quantiles of $e^{\eta }$ and $\beta$ estimate the requested marginal posterior distributions.

To compare kernels, introduce $M\in \{1,2\}$ with positive model prior distribution probabilities $w_1,w_2$. Under model two replace $A_{it}$ by ${\sum }_{j\in I_t}{e}^{-\gamma d_{ij}}$, defining $L_2(\alpha ,\gamma )$. The joint model/parameter target iswhere $\zeta =\beta$ or $\gamma$ and $p_M$ is its normalized proper prior distribution. Within each model use the preceding random-walk updates. For a concrete model switch, keep $\alpha$, propose the other model's kernel parameter independently from its prior distribution, and use equal forward and reverse switch-selection probabilities. A move from $(1,\alpha ,\beta )$ draws $\gamma \sim p_{\gamma }$ and has ratioAccept with the minimum of this ratio and one. The auxiliary-coordinate swap $(\alpha ,\beta ,\gamma )\mapsto (\alpha ,\gamma ,\beta )$ has absolute Jacobian determinant one, so no additional factor occurs; unequal switch-selection probabilities would require their reverse/forward ratio. This is a valid model-switching Markov chain Monte Carlo construction even though the two active models have equal parameter dimension. Prior distribution proposals may be inefficient, so fitted proposals can be substituted provided their probability densities remain in the acceptance ratio.

After checking mixing, estimate posterior model probabilities by the fractions of retained iterations in each model. The Bayes factor in favor of the first kernel isUse the posterior distributions within each model as well as these odds to assess the competing distance kernels. Proper prior distributions and their model-dependent normalizing constants are essential for this comparison.

For minimizing an energy $F(\theta )$, simulated annealing at temperature $T>0$ uses the Boltzmann distribution ${\pi }_T(\theta )\propto e^{-F(\theta )/T}$. At each temperature run a Metropolis–Hastings algorithm or Gibbs transition preserving that law, then reduce $T$. With a symmetric proposal, accept ${\theta }^{\prime }$ with probability $\min \{1,e^{-[F({\theta }^{\prime })-F(\theta )]/T}\}$: downhill moves are accepted and uphill moves remain possible. Keep the best state encountered. As temperature decreases, equilibrium mass favors global minima, provided the target is proper and the relevant concentration hypotheses hold. An arbitrary finite cooling run does not by itself prove global convergence.

Let $s={\sum }_{i=1}^m{x}_i$, $\bar{x}=s/m$ and $S={\sum }_i(x_i-\bar{x}{)}^2$. For the Poisson distribution model the observations must be nonnegative integers. Its log-likelihood isFor $s>0$, differentiation gives ${\ell }_P^{\prime }=-m+s/\lambda$ and ${\ell }_P^{\prime \prime }=-s/{\lambda }^2<0$. If $s=0$, the likelihood is decreasing and the maximizer is the boundary value zero. For the normal distribution model write $v={\sigma }^2>0$, givingFor each $v$, this is maximized at $\mu =\bar{x}$, and differentiating there gives $v=S/m$. HenceThe normal variance estimate has denominator $m$, not $m-1$. A finite positive normal maximum requires $S>0$; if all observations coincide, the likelihood is unbounded as $\mu =\bar{x},v\downarrow 0$. The normal annealing calculation below assumes $S>0$.

To maximize likelihood, choose the energy $-\ell$, so the annealing law is $L^{1/T}$, not $L^{-1/T}$. This sign is important in the Poisson request: literally inserting a positive log-likelihood into a minimizing Boltzmann law would give a nonnormalizable probability density function growing like $e^{m\lambda /T}$ at infinity. We use likelihood-power annealing with the maximizing sign.

For the normal distribution model, take base measure $d\mu dv$. The joint annealing probability density function isCompleting the square gives the Gibbs sampling conditional for $\mu$. Matching the conditional $v$ probability density function to the supplied inverse-gamma distribution convention givesIntegrating out $\mu$ contributes a factor $\sqrt{v}$, so the marginal law isThus this joint Boltzmann law is proper for $T<m/3$. Its variance marginal and the normal conditional both concentrate at the stated estimates as $T\downarrow 0$. For example, for $T<m/7$ its marginal mean and variance areAlso ${\mathbb{E}}_T(\mu -\bar{x}{)}^2=(T/m){\mathbb{E}}_{T}v\to 0$.

There is a direct convergence proof for the actual changing-temperature Gibbs algorithm, beyond convergence of its equilibrium probability densities. Start at a finite deterministic $v_0>0$, choose deterministic $0<T_n\le m/10$ with $T_n\to 0$, draw ${\mu }_n$ conditional on $v_{n-1}$, then draw $v_n$ conditional on ${\mu }_n$ using the two boxed distributions. Conditional normal moments giveThe inverse-gamma distribution first and second moments consequently yieldandIn the first recursion the coefficient of the preceding mean is at most $1/6$; in the second the coefficient of the preceding second moment is at most $1/8$. The means therefore stay bounded, and then so do the second moments. Letting $T_n\to 0$ in these formulas gives $\mathbb{E}{v}_n\to S/m$ and $\mathbb{E}{v}_n^2\to (S/m{)}^2$. Together with $\mathbb{E}({\mu }_n-\bar{x}{)}^2=T_n\mathbb{E}{v}_{n-1}/m\to 0$, this provesThis Gaussian likelihood Gibbs annealing convergence uses the explicit conditional structure and does not assume that one Gibbs sampling sweep is an exact equilibrium draw.

For the Poisson distribution model, the likelihood-power law with respect to $d\lambda$ has the normalized probability density functionsoThe mean tends to $\bar{x}$ and the variance to zero, proving convergence to a point mass at the Poisson maximum-likelihood estimate, including the zero-boundary case.

The Akaike information criterion is $2k-2\ell (\hat{\theta })$, counting one parameter for Poisson and two for normal. With $0\log 0=0$ and $S>0$, the two forms areandChoose the smaller criterion. Likelihood constants that differ between models, including the factorial and $2\pi$ terms, cannot be discarded.

For a fully specified trans-dimensional annealing for penalized likelihood, let $M=P$ or $N$ and defineUse the combined target, relative to counting measure on models and Lebesgue measure on their parameters,where $w_P,w_N>0$ are fixed reference weights. For $S>0$ and $0<T<m/3$ both model components are integrable. The factor $2T$ is a convenient temperature convention keeping the within-model laws exactly those already derived. The extra factor $e^{-k_M/T}$ makes the zero-temperature objective AIC rather than unpenalized likelihood.

Use, for example, probability $1/2$ for a model-change attempt in either model and probability $1/2$ for a within-model update. Within the Poisson distribution model draw a fresh gamma value from its boxed law; within the normal distribution model perform the boxed Gibbs sampling sweep. For a Poisson-to-normal proposal, draw $u\sim q(u)=N(0,{\tau }^2)$ with a fixed ${\tau }^2>0$, and setThe inverse normal-to-Poisson proposal is $\lambda =v$ with reverse auxiliary $u=\mu -v$. This dimension matching uses one auxiliary variable to augment the one-dimensional Poisson parameter space. The absolute Jacobian determinant isIf the forward model-jump selection probability is $b$ and the reverse probability is $d$, the reversible-jump Markov chain Monte Carlo acceptance ratio isAccept with probability $\min (1,R_{P\to N})$. For the reverse move compute $\lambda =v,u=\mu -v$ and useFor the suggested equal selection probabilities, $b=d=1/2$ and those factors cancel. The probability density function $q$ must still be included, despite symmetry of its normal distribution: it is an auxiliary draw in only the forward dimension-increasing proposal. Pairing a move with its inverse makes the accepted forward and reverse probability flows equal after the change of variables, proving detailed balance for each fixed-temperature model-switch kernel.

Cool through proper temperatures, allow enough model and parameter exploration at each stage, and retain the smallest observed $E_M$. The minima of the two energies are exactly the two AIC values above. Low-temperature equilibrium mass concentrates on the better criterion and its fitted parameters; if they tie, the minimizing set contains both models. Unlike the direct normal Gibbs argument, this multimodel exploration claim needs an appropriate cooling/mixing justification and is not guaranteed by an arbitrary short run.

Finally, an exact statistical comparison must put both likelihoods on the same observable and reference-measure convention. A Poisson distribution model is a probability mass function on counts, whereas the displayed normal distribution model is a continuous probability density function. The formulas above are the question's formal Gaussian working-likelihood comparison. A coherent exact normal approximation for count observations can instead use probabilities of rounding intervals; that changes its likelihood and generally loses the simple analytic normal estimates. This distinction matters when interpreting the numerical AIC difference, while leaving the dimension-matching and annealing construction applicable to properly specified likelihoods.