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

Write $H={\int }_{\mathbb{R}}h(x)dx$, and assume $0<H<\infty$; a uniform distribution on the region requires positive area. Set $x=V/U$, $s=U$, so the inverse transformation is $U=s$, $V=sx$. Its Jacobian determinant has absolute value $s$ for $s>0$. The region becomes $x\in \mathbb{R}$, $0<s\le \sqrt{h(x)}$, and its area isThe original joint probability density function is therefore $2/H$ on $C_h$. The change of variables formula givesThe boundary $U=0$ has probability zero, so division by $U$ causes no ambiguity for the resulting random variable.

Integrating out the second coordinate gives the probability density function of the marginal distributionThus the ratio has the normalized target probability density function. The factor $s$ from the Jacobian determinant is essential; merely transforming the rectangular coordinates as if they remained uniform would give a different probability density function. This is the basis of the ratio-of-uniforms method.

Here $h(x)$ attains its maximum one at $x=1$, so $a=1$. For $x\le 0$, put $t=-x\ge 0$. Then $x^{2}h(x)=e^{-1/2}{t}^2{e}^{-t/2}$. Its derivative is $e^{-1/2}{e}^{-t/2}t(2-t/2)$, so its maximum is at $t=4$, with value $16e^{-5/2}$.

For $0\le x\le 1$, $x^{2}h(x)=x^2{e}^{(x-1)/2}$ is increasing and has maximum one. For $x\ge 1$, $x^{2}h(x)=x^2{e}^{-(x-1)/2}$ has logarithmic derivative $2/x-1/2$, so its maximum is at $x=4$, with value $16e^{-3/2}>1$. Therefore the Laplace ratio-of-uniforms envelope has the finite boundsThese calculations establish finiteness and identify the tight axis-aligned bounds, rather than relying only on exponential tail decay.

Independently generate $U$ uniform on $(0,1)$ and $V$ uniform on $(b_1,b_2)$, using the bounds above. Accept the proposal ifotherwise generate a fresh pair. On acceptance return $X=V/U$. The bounding rectangle contains the entire ratio region: $U\le \sqrt{h(X)}\le a$, and $V=XU$ lies between the extrema of $x\sqrt{h(x)}$. Conditioning a uniform rectangle point on acceptance makes it uniform on $C_h$, so the preceding marginal distribution calculation proves that the returned value has the desired Laplace distribution.

For this target $H=4$, so the region has area two. The rejection sampling acceptance probability isEndpoints have zero probability and can be excluded to avoid divisions by zero in implementation.

Use the quantile function $F^{-1}(u)=\inf \{x:F(x)\ge u\}$ for $0<u<1$, so a strictly increasing CDF is not required. For any real $x$, the quantile event $F^{-1}(U)\le x$ agrees, up to endpoint null sets, with $U\le F(x)$. HenceThus $F^{-1}(U)$ has distribution function $F$ and probability density function $f$. This proves inverse transform sampling, including a probability density function with intervals on which the CDF is constant.

Integration of the Laplace distribution probability density function givesSolving these two equations for $x$ givesGenerate one uniform $U\in (0,1)$ and return this piecewise expression. The branches agree at $U=1/2$ and the inverse transform sampling proof establishes the target law without any rejection step.

Assume independent samples $X_i\sim g$, $\int f|\theta |<\infty$, and $g>0$ wherever $f\theta$ can contribute. The importance weight is $w(x)=f(x)/g(x)$, and the ordinary importance sampling estimator isIts expectation is $\int \theta (x)f(x)dx=\mu$, by cancellation of $g$ in the expectation integral. The strong law of large numbers gives almost-sure convergence to $\mu$. Common probability support is sufficient; the precise support condition for importance sampling only requires covering the nonzero integrand for this particular expectation.

Write $f=\tilde{f}/Z_f$ and $g=\tilde{g}/Z_g$. If the constants needed to evaluate $f/g$ are unavailable, use the self-normalized importance sampling estimatorThe unknown common factor $Z_g/Z_f$ cancels. Under full target probability support and integrability, the numerator and denominator sample means converge to $\int \tilde{f}\theta /Z_g$ and $\int \tilde{f}/Z_g$, respectively, whose ratio is $\mu$. The same construction applies if only one normalization is unknown by taking the other probability density function already normalized. Sampling from $g$ itself is assumed available, even when its normalization is not explicitly known.

This ratio estimator is generally biased at finite $n$, unlike the ordinary normalized-weight estimator. Also, its variance is not the formula in the next part: that formula concerns the ordinary unbiased estimator. Full probability support for $f$ is needed here because the denominator estimates the entire target normalizing integral, not just the probability support of $f\theta$.

Let $Y_i=\theta (X_i)f(X_i)/g(X_i)$. The independent summands have expectation $\mu$ and second momentWhen this integral is finite, independence givesIf the second-moment integral diverges, the estimator can still be unbiased and consistent, but its variance is infinite; the finite-variance assumption cannot be omitted.

Put $A=\int |\theta (x)|f(x)dx$, with $0<A<\infty$. The Cauchy-Schwarz inequality applied to $|\theta |f/\sqrt{g}$ and $\sqrt{g}$ givesEquality holds exactly when $|\theta |f/\sqrt{g}$ is proportional to $\sqrt{g}$ on the contributing probability support. Therefore the minimum-variance importance distribution isIf $\theta$ has one sign, $A=|\mu |$ and the ideal proposal gives zero variance. Its normalization or exact sampling may nevertheless be impractical. If $A=0$, the integrand vanishes almost everywhere and the estimator is identically zero.

The minimizing probability density function may vanish where $\theta =0$ even though $f>0$ there. This is valid for the ordinary estimator, which needs probability support only for its integrand. If one insists on exactly the same positive probability support as $f$, the bound is instead an infimum whenever this vanishing occurs, approached by $g_{\delta }=(1-\delta )g_0+\delta f$ as $\delta \downarrow 0$: because $g_{\delta }\ge (1-\delta )g_0$, the second moment lies between $A^2$ and $A^2/(1-\delta )$, proving convergence to the bound.

Take the indicator function $\theta (x)=1_{\{x\ge k\}}$. Its expectation under the standard Cauchy distribution is the desired upper-tail probability. To use the proposed tail probability density function in the remaining parts, assume $k>0$: a probability density function proportional to $x^{-2}$ on $[k,\infty )$ is not normalizable when that interval contains zero.

For $k>0$, ${\int }_k^{\infty }{x}^{-2}dx=1/k$. Thus the unnormalized integral is $1/k$ and the multiplicative normalization is $k$:Although this proposal does not cover the whole Cauchy distribution probability support, it covers the entire nonzero integrand $f\theta$, which is the required support condition for importance sampling.

For $x\ge k$, the proposal CDF is $G(x)=1-k/x$. Therefore inverse transform sampling with $V\sim U(0,1)$ givesEquivalently, put $U=1-V$, which is also uniform, and use $X=k/U$. This equivalent parametrization matches the uniforms used in the final estimator.

For $X_i=k/U_i$, the indicator is always one, and the importance weight reduces toAveraging these independent contributions givesAs a direct normalization check, its expectation is ${\int }_0^{1}k/[\pi (u^2+k^2)]du={\pi }^{-1}\arctan (1/k)$. For $k>0$ this equals $1/2-{\pi }^{-1}\arctan k$, the exact standard Cauchy distribution tail probability. The substitution $X=k/U$, rather than $k/(1-U)$ with the same symbol left unchanged, explains why the displayed denominator contains $U^2$.

For a target probability density function $\pi ({\theta }_1,\dots ,{\theta }_p)$, one systematic sweep of the Gibbs sampler updates coordinates successively:Each update leaves the target invariant: keeping the other coordinates with their target marginal and replacing one coordinate by its conditional distribution reconstructs the target joint probability density function. Therefore composition of the coordinate updates also preserves $\pi$. It defines a Markov chain whose consecutive draws are generally dependent. Under irreducibility, aperiodicity and the usual recurrence conditions it converges to the target law from an arbitrary valid start; a stationary start gives exact target marginals immediately. Invariance alone does not prove these ergodicity conditions, and a systematic full sweep need not itself be reversible.

Write $v={\sigma }^2$, $d=k+1$, $X=X_k$, $a=a_k$, $\mu ={\mu }_k$ and $\Sigma ={\Sigma }_k$, with $\Sigma$ positive-definite. For fixed $v$, the terms involving $a$ in the posterior exponent areComplete the square. With $V_k=(X_k^T{X}_k/v+{\Sigma }_k^{-1}{)}^{-1}$ and $m_k=V_k(X_k^{T}y/v+{\Sigma }_k^{-1}{\mu }_k)$, this givesFor fixed $a_k$, collect the powers and exponential terms involving $v$:Hence, in the paper's shape-scale inverse-gamma distribution convention,There is no extra $(k+1)/2$ in this shape, because the normal coefficient prior is independent of $v$. These are the independent normal and inverse-gamma regression priors conditionals.

The Gaussian normalization printed in the PDF is incorrect. The proper $d$-dimensional prior probability density function is $(2\pi {)}^{-d/2}|{\Sigma }_k{|}^{-1/2}\exp [-(a_k-{\mu }_k{)}^T{\Sigma }_k^{-1}(a_k-{\mu }_k)/2]$. For fixed $k$ this correction is constant in $a_k,v$ and leaves the two conditionals unchanged, but it must be retained when comparing different model orders in part (b).

Start with $v^{(0)}>0$ and any coefficient vector. At sweep $m$, draw $a_k^{(m+1)}$ from the multivariate normal distribution above using $v^{(m)}$, then draw $v^{(m+1)}$ from the inverse-gamma distribution using the newly drawn coefficients. For the second step one may generate $1/v$ from a gamma distribution with the same shape and rate equal to the displayed inverse-gamma scale.

Each block update uses an exact conditional distribution, so it preserves the joint posterior distribution. The resulting pairs form a dependent Gibbs sampler sample. With proper priors, positive scale and positive-definite coefficient covariance, these conditionals have full probability support on ${\mathbb{R}}^{k+1}\times (0,\infty )$ and the usual posterior ergodicity applies. A burn-in reduces dependence on initialization; retained sweeps are still serially dependent and Monte Carlo errors should account for that dependence.

The order prior and probabilities of choosing upward or downward moves are not specified, so include them explicitly. Let ${\rho }_k$ be the model-order prior probability, $b_k$ the birth-selection probability at order $k$, and $d_{k+1}$ the reverse death-selection probability. Denote the correctly normalized coefficient prior by $p_k(a)$ and the probability density function of $z\sim N(0,{\sigma }_a^2)$ by $q(z)$. Common variance-prior factors cancel since $v={\sigma }^2$ is unchanged.

Put $m=n^{-1}{\sum }_i{x}_i^{k+1}$, $r=y-X_{k}a$, and $w_i=x_i^{k+1}-m$. The centered polynomial birth move in reversible-jump sampling gives $X_{k+1}{a}^{\prime }=X_{k}a+zw$, so the residual-sum-of-squares change isThe map $(a_0,a_1,\dots ,a_k,z)\mapsto (a_0-mz,a_1,\dots ,a_k,z)$ is triangular with determinant one. Thus the reversible-jump Markov chain Monte Carlo ratio isFor an entirely explicit Gaussian expression, define $Q_k=(a-{\mu }_k{)}^T{\Sigma }_k^{-1}(a-{\mu }_k)$ and $Q_{k+1}=(a^{\prime }-{\mu }_{k+1}{)}^T{\Sigma }_{k+1}^{-1}(a^{\prime }-{\mu }_{k+1})$. Dividing the two normalized Gaussian priors and the proposal probability density function givesEqual model-order priors and equal birth/death selection probabilities remove the first ratio, if that is the intended convention. Dropping the order-dependent normal-prior constants instead would give the wrong acceptance rule.

From order $k+1$, set $z=a_{k+1}^{\prime }$, remove that final coefficient, retain $a_i=a_i^{\prime }$ for $1\le i\le k$, and restore the intercept byKeep $v$ fixed. This deterministically recovers the lower-order state and the auxiliary variable that would reproduce the higher-order state under a birth. The inverse map also has absolute Jacobian determinant one. In particular, simply deleting the highest coefficient without restoring the intercept is not the reverse of the supplied move.

For the recovered $a,z$ of the preceding part, compute the forward ratio $R$ from part (i). The reverse probability density function ratio and reverse Jacobian determinant invert every factor, soEquivalently, with $L_k(a;v)$ denoting the Gaussian regression likelihood,The factor $q(z)$ is present even though no random variable is generated during the death proposal: it is the probability density function of the auxiliary variable required by the reverse birth. These reciprocal rules give detailed balance on the paired dimension-matching moves.

Put $v={\sigma }^2>0$, $\bar{x}=n^{-1}{\sum }_i{x}_i$ and $S={\sum }_i(x_i-\bar{x}{)}^2$. The likelihood function and log-likelihood areFor every fixed $v$, the unique maximum in $\mu$ is $\bar{x}$. At that value, $\partial \ell /\partial v=-n/(2v)+S/(2v^2)$ changes sign from positive to negative at $v=S/n$. Thus, for $S>0$,This is the global maximum-likelihood estimate, since the first maximization was global for each $v$ and the second has a unique peak. If $S=0$, the likelihood function is unbounded at $\mu =\bar{x}$, $v\downarrow 0$; there is no maximizer in the parameter space $v>0$. The freezing argument below uses the nondegenerate case $S>0$.

At temperature $T>0$, target a probability density function proportional to $\exp (f(\theta )/T)$ on a parameter region where it is normalizable, possibly with a specified proper base measure. A Metropolis–Hastings algorithm proposal $q(\theta ,{\theta }^{\prime })$ has acceptance probabilityFor a symmetric proposal, improving moves are accepted and a loss $\Delta$ is accepted with probability $e^{-\Delta /T}$. Run updates while gradually lowering $T$ toward zero. This simulated annealing scheme initially permits exploration and increasingly favors high objective values. Appropriate cooling and mixing conditions are needed for a general global convergence theorem; a finite arbitrary schedule does not guarantee finding a global maximum.

To maximize a positive likelihood function, take $f=\log L$, which has the same maximizers. The resulting likelihood-power annealing target is $L^{1/T}$, and available Gibbs sampler conditional distributions can replace Metropolis proposals.

Use the Gaussian likelihood Gibbs annealing target ${\pi }_T(\mu ,v)\propto L(\mu ,v{)}^{1/T}$ relative to $d\mu dv$, with $S>0$. Completing the normal square and recognizing an inverse-gamma kernel giveThe shape is $n/(2T)-1$, because the probability density function in $v$ is $v^{-n/(2T)}{e}^{-[S+n(\mu -\bar{x}{)}^2]/(2Tv)}$. The base measure matters: inserting a prior factor $1/v$ would change this shape.

Integrating out $\mu$ contributes a factor proportional to $v^{1/2}$, so the marginal is $v\sim \mathrm{IG}(n/(2T)-3/2,S/(2T))$. Thus the joint target is proper for $T<n/3$. Choose a decreasing sequence $0<T_m\le n/10$ tending to zero and a finite starting value $v_0>0$. At iteration $m$, sample ${\mu }_m$ from the first conditional using $v_{m-1}$, then sample $v_m$ from the second using ${\mu }_m$. This is an explicit annealing algorithm entirely using Gibbs sampling, with all its inverse-gamma moments finite for the convergence proof.

We prove convergence of the actual sequential updates above, not just concentration of ideal equilibrium distributions. Write $M_m=\mathbb{E}{v}_m$ and $N_m=\mathbb{E}{v}_m^2$. Given $v_{m-1}$, the centered normal draw satisfiesThe conditional inverse-gamma distribution moments therefore giveFor $T_m\le n/10$, the coefficient multiplying $M_{m-1}$ is at most $1/6$, so $M_m$ stays bounded by a geometric-recursion argument. The coefficient multiplying $N_{m-1}$ is at most $1/8$, and the other terms are bounded using the first-moment bound, so $N_m$ also stays bounded. These bounds start from the finite deterministic $v_0$.

Since $T_m\to 0$, the boundedness makes every term containing $T_m{M}_{m-1}$ or $T_m^2{N}_{m-1}$ vanish. Hence $M_m\to S/n$ and $N_m\to S^2/n^2$. ConsequentlyThusThis establishes mean-square convergence for the specified Gaussian annealer. It does not assert that every implementation or cooling rule for every objective has the same guarantee, nor that mean-square convergence alone is almost-sure convergence.

For observed data $x$ and latent variables $Z$, let ${\ell }_c(\theta ;x,Z)=\log p_{\theta }(x,Z)$ be the complete-data log-likelihood. From a current parameter ${\theta }^{(m)}$, the expectation-maximization algorithm computesand then chooses ${\theta }^{(m+1)}$ maximizing this function in its first argument. Thus the E-step averages the missing information conditional on the current model and observed data; the M-step treats those expectations as fixed while maximizing.

For a short proof of likelihood function monotonicity, let $q_m(Z)=p_{{\theta }^{(m)}}(Z\mid x)$. ThenThe Jensen inequality for the logarithm gives $\ell (\theta ;x)-\ell ({\theta }^{(m)};x)\ge Q(\theta \mid {\theta }^{(m)})-Q({\theta }^{(m)}\mid {\theta }^{(m)})$. The M-step makes the right side nonnegative. Thus observed likelihood function cannot decrease. This proves the usual EM likelihood monotonicity, while allowing convergence to local stationary points rather than guaranteeing a global maximum.

Put ${\alpha }_1=\alpha$, ${\alpha }_2=1-\alpha$, and introduce a categorical class label for each observation, encoded by $Z_{i1},Z_{i2}\in \{0,1\}$ with $Z_{i1}+Z_{i2}=1$. Its prior class probability is ${\alpha }_j$, and conditional on class $j$ the observation probability density function is $f_j(x_i)$. Thus its joint observation-label probability density function isIndependence across observations yields the finite Gaussian mixture with a common variance complete-data likelihood functionSumming over the two possible labels at each observation recovers ${\prod }_i[\alpha f_1(x_i)+(1-\alpha )f_2(x_i)]$, the observed-data likelihood function. The product requested in the PDF is therefore a joint complete-data likelihood function. Strictly conditioning on $Z$ would remove the factors ${\alpha }_j$; the PDF's conditional notation for this joint product is imprecise.

The E-step uses the posterior class probabilitiesLet $n_j^{(m)}={\sum }_i{w}_{ij}^{(m)}$, so $n_1^{(m)}+n_2^{(m)}=n$. Up to constants, the expected complete-data log-likelihood isMaximizing the first term with ${\alpha }_1+{\alpha }_2=1$ gives ${\alpha }^{(m+1)}=n_1^{(m)}/n$. Differentiating in each ${\mu }_j$ gives ${\mu }_j^{(m+1)}={\sum }_i{w}_{ij}^{(m)}{x}_i/n_j^{(m)}$. After these updates, maximizing in the shared variance givesThe residuals use the new component means, while the weights remain from the old parameters. The divisor is $n$, not $n-2$, because this is a maximum-likelihood update. Starting with $0<\alpha <1$ and $v>0$ gives positive responsibilities and positive $n_j$; a zero effective component size at a boundary would leave its mean unidentified. These are the EM for Gaussian mixtures with a common variance transitions.

For a weakly stationary process with finite second moments, let $\mu =\mathbb{E}{X}_t$ and $\gamma (h)=\mathrm{Cov}(X_{t+h},X_t)$, independent of $t$. If $\gamma (0)>0$, its autocorrelation function is $\rho (h)=\gamma (h)/\gamma (0)$. Given observations $x_1,\dots ,x_N$, one usual sample autocovariance function isThe correlogram plots $\hat{\rho }(h)$ against lag $h$. The partial autocorrelation function at lag $h$ is the correlation between $X_t$ and $X_{t-h}$ after linearly projecting each onto the intervening observations. Its sample version can be calculated by solving the Yule-Walker regression equationsthen taking the last coefficient ${\hat{\varphi }}_{h,h}$. Equivalently, use the corresponding residual correlation from the sample regressions, with a consistent covariance convention.

For a minimal causal autoregressive model of order $p$, the population partial autocorrelation function is zero after lag $p$, while the autocorrelation function follows an exponentially decaying or damped oscillatory tail. For a moving-average model of order $q$, its autocovariance and hence autocorrelation are zero after lag $q$, because observations farther apart share no driving-noise terms; its partial correlations usually tail off. These order identification by autocorrelation cutoffs patterns suggest AR or MA candidate orders. Sample cutoffs are approximate and should be judged with sampling uncertainty and residual diagnostics; the familiar white-noise bands near $\pm 1.96/\sqrt{N}$ are a null approximation, not universal confidence limits for a dependent series.

For the trend problem, write a symmetric estimator asReproducing arbitrary constant and quadratic trends requires $2a+2b+c=1$ and $8a+2b=0$; the linear-trend condition holds automatically by symmetry. These constraints follow by expanding $T_{t+j}=T_t+j(b_0+2c_{0}t)+c_0{j}^2$ for a trend $T_t=a_0+b_{0}t+c_0{t}^2$. Thus $b=-4a$, $c=1+6a$. Since the errors are uncorrelated and the estimator is unbiased, its mean squared error isThe strictly convex quadratic is minimized at $a=-3/35$, giving the five-point symmetric quadratic trend smootherThis weighted moving average allows negative coefficients; a convex average could not cancel the quadratic moment.

For three points, symmetry gives weights $(a,b,a)$. Unbiasedness for every quadratic requires $2a+b=1$ and $2a=0$, so the only solution is $(0,1,0)$. Thus there is no nontrivial three-point quadratic-preserving smoother. The literal claim that there is no unbiased estimator at all is false, since ${\hat{T}}_t=X_t$ is unbiased with error variance ${\sigma }^2$. This is the three-point quadratic reproduction forces the identity filter distinction.

For seasonal extraction by a centred moving average, use the additive decomposition $X_t=T_t+S_t+{\epsilon }_t$, where $S_t$ has period $m$ and its values over one period sum to zero. An average covering a full period removes the seasonal component and approximately retains a slowly changing trend. If $m$ is odd, take the centred length-$m$ average; if $m$ is even, centre it by averaging two adjacent length-$m$ averages, giving half weights at the two endpoints. Subtract this trend estimate from the data, average residuals within each seasonal phase over different cycles, and centre the resulting phase estimates to sum to zero. For a multiplicative decomposition, use ratios and normalize the phase estimates to average one.