Past exam of the mathematics course of the University of Cambridge / 2017 / iii / Paper 206 (Applied Statistics)

Write the exponential dispersion family in the conventionAssume that the support does not depend on the parameters, that $\theta$ is in the interior of the natural parameter space, and that differentiation can pass under the normalizing integral or sum. Differentiating $\int f=1$ at fixed dispersion parameter givesso the expectation is $b^{\prime }(\theta )$. A second differentiation givesConsequentlyHere $b$ is the cumulant function of an exponential family; the same identities follow by differentiating $\log \mathbb{E}{e}^{tZ}=[b(\theta +\varphi t)-b(\theta )]/\varphi$. For a known observation weight $w$, replace $\varphi$ by $\varphi /w$. These regularity assumptions for differentiation under the integral sign matter: differentiation of a parameter-dependent support can introduce additional boundary terms.

For fixed known $r>0$, take the natural parameter of an exponential familyThe negative binomial distribution has log-likelihoodThus this is a negative binomial exponential family, with dispersion parameter $\varphi =1$ in the convention of part (a). The carrier contains the Gamma function and depends on the fixed size $r$, which need not be an integer. Differentiating the cumulant function of an exponential family givesThereforeThe size $r$ controls overdispersion relative to a Poisson distribution; it is different from the unit dispersion parameter used by this fixed-size generalized linear model. Another common convention applies to the scaled response $Z=Y/r$: then $b(\theta )=-\log (1-e^{\theta })$ and $\varphi =1/r$. Keeping track of the response scale prevents confusing these two representations. At $\lambda =0$ the law is concentrated at zero and is obtained as a boundary limit, rather than by a finite natural parameter of an exponential family.

The negative binomial regression assumes independent random variables conditional on their recorded days, withThe logarithmic link function applies to the conditional expectation, not to the observed count. The displayed maximum-likelihood estimates areand the dispersion parameter for the fitted family is one. Thus the fitted variance is ${\hat{\lambda }}_i+{\hat{\lambda }}_i^2/\hat{r}$, practically the Poisson distribution variance at these means.

Increasing the day by one multiplies the conditional expectation by $e^{{\hat{\beta }}_1}\simeq 0.996836$, a decrease of about $0.3164\%$. Comparing days 365 and 1 gives $e^{364{\hat{\beta }}_1}\simeq 0.3155$. The negative slope has Wald statistic approximately $-35.79$, so under the stated model assumptions there is very strong evidence of a decreasing conditional mean:This is an association with calendar day, not a causal effect established by the regression. The printed negative-binomial residual deviance is $7544.3$ on $1298$ residual degrees of freedom, which is a reason to investigate the adequacy of the assumed count model and dependence structure. In particular, the small reported slope standard error cannot make a misspecified model reliable.

Introduce a latent Bernoulli distribution indicator $A_i$, with $A_i=1$ denoting a structural zero and $\mathrm{Pr}(A_i=1)=\pi$. Conditional on $A_i=0$, retain the negative binomial regression with ${\lambda }_i=e^{x_i^T\beta }$ and $x_i=(1,d_i{)}^T$. This gives a zero-inflated negative binomial model:A participant who fishes may still catch zero, so an observed zero does not reveal the latent class. This distinguishes zero inflation from a hurdle model, which truncates the count component at zero. The marginal expectation becomes $(1-\pi ){\lambda }_i$, and the law of total variance givesA constant $\pi$ is a parsimonious starting model. If there is evidence that nonparticipation changes with day or other recorded predictors, use logistic regression for ${\pi }_i$ instead. Conditional independence of visitors remains an assumption; extra zeros alone do not establish its validity.

Let $I_0=\{i:y_i=0\}$ and $I_{+}=\{i:y_i>0\}$. Under the proposed zero-inflated negative binomial model, independence gives the likelihood functionThe positive-count factor uses the ordinary, untruncated negative binomial distribution. With $f_i=f_{\mathrm{NB}}(y_i;r,{\lambda }_i)$, the log-likelihood isFor latent indicators $a_i$, the complete-data log-likelihood iswhere $a_i=1$ is allowed only when $y_i=0$. This separation is what makes the expectation-maximization algorithm convenient.

For EM for zero-inflated negative binomial regression, at iteration $k$ apply Bayes' theorem to each zero in the zero-inflated negative binomial model. The E step computes the conditional expectation of its latent indicator:Treating these values as fixed, the M step maximizes the expected complete-data log-likelihood. Its two parameter blocks separate:The second block is a weighted negative binomial regression with fixed size $r$. Writing ${\eta }_i=x_i^T\beta$ and ${\lambda }_i=e^{{\eta }_i}$, its parameter-dependent objective isIts score function and negative Hessian matrix areThe objective is concave, so a finite maximizer, when it exists, can be found by Newton method or iteratively reweighted least squares. The latter uses working response ${\eta }_i+(y_i-{\lambda }_i)/{\lambda }_i$ and weights $w_{i}r{\lambda }_i/(r+{\lambda }_i)$. This is not an unweighted regression on only the positive counts: zeros retain fractional weight in the count component.

Start with $0<\pi <1$ and finite coefficients, alternate the two steps, and monitor the observed log-likelihood. Exact maximization makes the expectation-maximization algorithm nondecreasing in that likelihood, although the mixture may have multiple stationary points. All-zero data can lead to boundary fits and poor identification; multiple initializations and explicit checks for boundary parameters are useful.

For independent random variables $Y_i\sim \mathrm{Bernoulli}(p_i)$, the full logistic regression has linear predictorLayout A is the reference level in a regression factor. The smaller model sets ${\beta }_B={\beta }_C=0$, leaving the same logit link and the price predictor. In either caseThe saturated statistical model allows a separate probability per observation; its maximizing probabilities are $p_i=y_i$, with log-likelihood zero under $0\log 0=0$. Hence the binomial deviance for either fitted model isEquivalently, write the summands as $2[y_i\log (y_i/{\hat{p}}_i)+(1-y_i)\log ((1-y_i)/(1-{\hat{p}}_i))]$. The fitted probabilities differ between the two models, but the definition of binomial deviance is the same. A residual deviance for individual binary data need not itself have an accurate chi-squared distribution approximation; the nested-model difference used next has a different asymptotic justification.

There are $120$ observations, since the intercept-only model has $119$ residual degrees of freedom. Adding price introduces one coefficient and decreases the binomial deviance by $148.263-101.920=46.343$. Adding a three-level factor introduces two coefficients and decreases it by $101.920-97.571=4.349$. Thus the missing entries areThe likelihood-ratio test comparing the smaller and full logistic regression tests $H_0:{\beta }_B={\beta }_C=0$ against at least one nonzero layout contrast. Its statistic is $2({\hat{\ell }}_1-{\hat{\ell }}_2)=4.349$. Under ordinary interior-parameter, full-rank and finite-estimate regularity conditions, Wilks theorem gives an asymptotic chi-squared distribution with two degrees of freedom. The reported p-value $0.1136$ exceeds $0.05$, so at that level retain the simpler price-only model.

The table is an analysis of deviance for nested generalized linear models: price is entered before layout. Thus its price row tests the price-only model against an intercept, whereas the layout row tests layout after allowing for price. The latter is the relevant comparison of model1 and model2. Failing to reject a layout effect is not proof that layouts have identical effects; it means these data do not justify the extra two parameters under this test.

Using the smaller logistic regression, approximate pointwise Wald confidence intervals are an estimate plus or minus $1.96$ times its standard error:Both intervals rely on the large-sample normal approximation and correct model assumptions. The negative slope indicates decreasing purchase probability as price rises. A one-unit price increase multiplies the odds by $e^{-0.08207}\simeq 0.9212$; it does not subtract a fixed amount from the probability. There is no significant additional layout effect in the previous likelihood-ratio test, rather than a demonstration of no effect.

In the selected model layout B has the same prediction rule as the other layouts. At price 100 the linear predictor is $5.16258-0.08207(100)=-3.04442$. Applying the inverse logit link givesThe intercept refers to price zero; if that lies outside the observed price range, its direct substantive interpretation would require extrapolation.

The quasibinomial regression retains the logit link and mean function $p_i=(1+e^{-x_i^T\beta }{)}^{-1}$ but uses the working variance function $\mathrm{Var}(Y_i)=\varphi p_i(1-p_i)$. Its quasi-score equation isso multiplying by the common dispersion parameter does not change the coefficient estimates. Iteratively reweighted least squares therefore gives the same fitted means and coefficients as the binomial model. The usual Pearson dispersion estimator isThe estimated coefficient covariance matrix is $\hat{\varphi }(X^{T}WX{)}^{-1}$, with $W_{ii}={\hat{p}}_i(1-{\hat{p}}_i)$. From the intercept standard errors, $(1.08921/1.13031{)}^2\simeq 0.92860$; the rounded slope errors give approximately the same factor. The display uses approximate Student's t-tests with 118 residual degrees of freedom instead of fixed-dispersion normal tests.

For a genuinely individual binary response, the Bernoulli distribution identity $Y_i^2=Y_i$ forces $\mathrm{Var}(Y_i)=p_i(1-p_i)$. Thus $\varphi \ne 1$ is a working quasi-likelihood specification, not a different independent binary distribution with that mean. The mild estimated underdispersion does not establish an improved model, and the mean predictions are unchanged. In particular, a scalar rescaling does not model correlation among users, and usual likelihood-based Akaike information criterion comparisons are unavailable for a family without a specified probability likelihood. The output gives no convincing reason to prefer model3 to model2.

With $K\ge 2$ distinct knots ${\kappa }_1<\cdots <{\kappa }_K$, a natural cubic spline is a $C^2$ function on the real line, cubic on each interval between consecutive knots and affine outside the two outer knots. The latter condition gives $m^{\prime \prime }({\kappa }_1)=m^{\prime \prime }({\kappa }_K)=0$.

There are $4(K-1)$ coefficients for the interior cubic pieces. Matching values and the first two derivatives at the $K-2$ interior knots imposes $3(K-2)$ independent constraints, and the two natural endpoint conditions impose two more. The affine continuations are then determined by endpoint values and slopes. The dimension of a vector space is thereforeOne explicit basis uses $1,x$ and $d_j(x)-d_{K-1}(x)$ for $j=1,\dots ,K-2$, whereThe cubic and quadratic terms cancel in each difference above ${\kappa }_K$, making it affine there; below ${\kappa }_1$ the truncated powers vanish. The changes in third derivative at individual knots show independence of these basis functions. This count uses the linear independence of the spline constraints and counts the listed outer knots as boundary knots. If a convention instead lists $K$ interior knots and two additional boundaries, the dimension is $K+2$. With only one listed knot and linearity on both sides, the $C^2$ condition forces a globally affine function, so the exceptional dimension is two.

Choose knots independently of the eventual validation responses, and let $b(x)=(b_1(x),\dots ,b_K(x){)}^T$ be a natural cubic spline basis as in part (a). A regression spline writes $m(x)=b(x{)}^T\beta$. Define the design matrix $B_{ij}=b_j(X_i)$. Under the stated normal linear model, maximum likelihood estimation for $\beta$ is the same as ordinary least squares:Differentiation gives the normal equations $B^{T}B\hat{\beta }=B^{T}Y$. If the design has full column rank,If it does not, a Moore-Penrose inverse gives the fitted values, but not every coefficient or off-sample prediction is identifiable. Merely choosing distinct knots does not guarantee full rank for an arbitrary set of observed predictor values.

More knots allow greater curvature and reduce approximation bias at the cost of larger estimation variance. Too many or poorly located knots can cause unstable fits. Equally spaced knots are simple; empirical predictor quantiles put more flexibility where observations are dense; subject knowledge can concentrate knots near likely changes. Both the number and positions can be compared using cross-validation or an appropriate penalized model-selection criterion. Adaptive knot selection must be included within validation, since treating a data-selected basis as prespecified understates its flexibility. Boundary knots should cover the region where inference is intended, and natural linear tails do not by themselves justify far-out extrapolation.

For $\lambda >0$, minimize the smoothing spline objective over the Sobolev space of functions with square-integrable second derivative on the observed interval:The roughness penalty measures curvature; its null space consists of affine functions. To see why the minimizer is a natural cubic spline, let $g$ be the natural spline interpolating the values of any competitor $m$ at its distinct observation locations, and put $h=m-g$. Integration by parts twice on each knot interval gives $\int g^{\prime \prime }{h}^{\prime \prime }=0$: $g^{\prime \prime \prime \prime }=0$ on each interval, $h$ vanishes at the knots, $g^{\prime \prime }$ is continuous, and $g^{\prime \prime }$ vanishes at the outer endpoints. Hencewhile both functions have the same residual sum of squares. This proves the reduction to a finite-dimensional spline space.

For distinct $X_i$ and $n\ge 2$, use its $n$-dimensional basis, write $B_{ij}=b_j(X_i)$, and define the positive semidefinite matrixThe penalized least squares problem is $\Vert Y-B\beta {\Vert }^2+\lambda {\beta }^T\Omega \beta$, soUse linear continuation outside the observation interval. The fitted-value matrix is $S_{\lambda }=B(B^{T}B+\lambda \Omega {)}^{-1}{B}^T$; its trace gives the effective degrees of freedom. As $\lambda$ decreases to zero, the natural spline approaches interpolation, with trace $n$; as it grows without bound, the fit approaches ordinary least squares on $1,x$, with trace two. For $\lambda =0$ the criterion alone does not determine values between observations; the natural interpolant is the limiting convention. Repeated observation locations are combined by using distinct knots and their multiplicities as weights, rather than pretending the spline design has $n$ independent interpolation values.

Because no family is specified, the command fits a Gaussian identity-link generalized additive model:The $s_j$ are centered penalized regression splines with cubic regression bases; centering separates them from the intercept. The age curve in the original figure falls to a minimum near age 25, then rises markedly, particularly from about 40 onwards. Its effective degrees of freedom are $4.05$, reflecting curvature. The white-cell curve has effective degrees of freedom one and is nearly flat with a slight positive slope; the pointwise uncertainty bands are consistent with no substantial effect. These plots display centered partial mean contributions, not numbers of prescriptions by themselves.

A simpler mean model keeps the age regression spline and replaces the white-cell smooth by a linear term. A separate issue is that the response is a nonnegative count: a Poisson regression with logarithmic link function respects that support and positivity of the mean. A reasonable candidate isIf diagnostics show overdispersion, one can fitusing a negative binomial regression. Alternatively, if a Gaussian approximation is adequate, the minimal simplification is gam(npres ~ s(age, bs="cr") + wbc). These are candidates to compare by diagnostics and validation, rather than guaranteed improvements from the partial plots alone. Retain the curved age effect; consider a linear or removable white-cell effect, and assess an appropriate count family. Refit before interpreting the old curves on a new link scale.

There is a genuine syntax error in the printed command: random=list(...) must be a separate argument, following a comma rather than a plus. The intended call isWith the default Gaussian family, this generalized additive mixed model isHere $(b_{0j},b_{1j}{)}^T\stackrel{\mathrm{iid}}{\sim }{N}_2(0,D)$ across doctors, $D$ is an estimated covariance matrix, and the Gaussian errors have variance ${\sigma }^2$ and are independent of the random effects. The formula ~age includes a random intercept as well as a random slope; it does not force these two effects to be uncorrelated.

The random intercept allows different prescribing baselines, and the random slope allows doctor-specific linear modifications to the population age curve. Conditional on these effects the observations are independent under the model, whereas two patients of the same doctor have additional covariance $(1,a_{ij})D(1,a_{kj}{)}^T$. This accounts for clustered data and avoids treating doctor variation as independent patient-level noise. The population smooth and doctor-specific predictions are different targets. With only four doctors, estimates of the random-effect covariance can be imprecise; the Gaussian support issue for counts also remains. The new model accounts for doctor-specific baselines and age slopes.

The random-slope linear mixed model iswith independence between the random effects and errors. The estimates areIncubators are treated as exchangeable representatives of possible growth environments, so a random slope permits their growth rates to differ while estimating a population mean growth rate. The formula 0+hours deliberately removes a random intercept: larvae are assigned only after hatching, so it is reasonable to assume no incubator-specific baseline at time zero. Random assignment supports this common-baseline assumption in expectation; it does not prove that observed initial sizes or other baseline differences are exactly identical.

For a new incubator, no observations are available to estimate its realized random slope, and its mean effect is zero. Under squared-error loss the best prediction integrates over that effect and the residual error:This is the plug-in conditional expectation given the population parameters, rather than a prediction using a fitted effect from one of the four existing incubators. Ignoring parameter-estimation uncertainty, the prediction variance for a new individual is $10^2{\hat{\tau }}^2+{\hat{\sigma }}^2=375.29023$; it is not just the uncertainty in the population mean.

Let $X$ have columns $1,t$, and let $Z$ have one column per incubator, with entry $t_i$ in the column for observation $i$'s incubator and zeros elsewhere. Integrating the Gaussian random effects givesConsequently its marginal log-likelihood isThe likelihood-ratio test of $H_0:{\tau }^2=0$ against ${\tau }^2>0$ usesFor each covariance, the optimizing coefficients are the generalized least squares estimator $(X^T{V}^{-1}X{)}^{-1}{X}^T{V}^{-1}Y$. Both fits should use the same marginal likelihood convention, for example fitting the alternative with update(fly.model, REML=FALSE) and the null as a Gaussian linear model. The displayed alternative REML criterion alone cannot supply $T$.

The variance component is constrained to be nonnegative, and zero lies on the boundary. Thus the interior-parameter hypothesis of Wilks theorem fails; the usual ${\chi }_1^2$ approximation is inappropriate. A $\frac{1}{2}{\delta }_0+\frac{1}{2}{\chi }_1^2$ limit sometimes applies to a single variance component under additional asymptotic conditions, but four incubators do not give a persuasive large-number-of-groups approximation.

Use a parametric bootstrap under the null: estimate $\beta ,{\sigma }^2$ there; simulate independent Gaussian errors with that fitted variance at exactly the observed times and grouping labels; refit both models by maximum likelihood for each simulated sample; and calculate $T^{\ast }$. The upper-tail proportion, conventionally $(1+\#\{T^{\ast }\ge T\})/(B+1)$, estimates the p-value. Boundary fits with estimated ${\tau }^2=0$ must remain in the simulation distribution. A restricted-likelihood ratio and a correspondingly calibrated bootstrap are another option when the fixed-effect design is identical, but one must not mix ordinary and restricted likelihoods in a single statistic.

The marginal Akaike information criterion is $\mathrm{AIC}=-2\ell ({\hat{\vartheta }}_{\mathrm{ML}})+2k$, where $k$ counts estimated population parameters. Here they are ${\beta }_0,{\beta }_1,{\sigma }^2,{\tau }^2$, so $k=4$. The realized incubator slopes have been integrated out of the likelihood function; they are not four additional free population coefficients. From the marginal Gaussian log-likelihood, the expression isAll hats here denote the ordinary maximum likelihood estimation fit.

The displayed value $14041.1$ is instead a restricted maximum likelihood criterion. Thus it cannot simply be substituted for the ordinary maximum-likelihood deviance. The logLik method defaults to the fitting convention. If one applies AIC(fly.model) directly to the displayed REML object, the default restricted log-likelihood and the same parameter count give the software's restricted-likelihood-based value $14041.1+8=14049.1$. This number should be labelled as such; it is not the ordinary marginal AIC derived above. To obtain the latter, refit and calculateThe output does not contain that numerical answer. In particular, restricted criteria must not be used to compare models with different fixed-effect designs. A conditional AIC targeting predictions for existing groups is yet another criterion and does not follow by counting each latent slope as an ordinary parameter.

Add a location-specific random intercept to the random-slope linear mixed model:Take the location effects, incubator slopes and individual errors to be mutually independent, with the original variances ${\tau }^2$ and ${\sigma }^2$. The random intercept models variation already present at hatching because of the origin of the eggs. Because larvae are subsequently assigned to incubators, location and incubator are crossed random effects, not automatically nested groups. The corresponding code isThe resulting covariance between two observations contains ${\omega }^2$ if they share a location and ${\tau }^2{t}_i{t}_k$ if they share an incubator. If the three named locations themselves are the only populations of interest, use a fixed location factor instead:Choosing between these interpretations depends on the inferential target; with only three locations, a location variance component will be weakly estimated.

The command fits a Gaussian generalized additive mixed modelThe population growth curve $s$ is a centered penalized regression spline with a cubic regression basis. The incubator term remains a random slope without a random intercept, because ~0+hours removes the intercept. A population linear component is allowed within the smooth's unpenalized part; curvature is penalized rather than being required.

A gamm fit returns a list with gam and lme components. Plot the population smooth and pointwise uncertainty bands withInspect whether the fitted curve departs materially from an affine function, particularly in regions supported by observations. The effective degrees of freedom near one suggest that the penalization has selected an approximately linear centered smooth. Bands and partial residuals help distinguish convincing curvature from noisy deviations, although pointwise bands are not a simultaneous test of linearity. A formal comparison would need to account for both smoothing selection and the incubator dependence. The smooth plot tests the visual plausibility of a common straight growth curve while retaining incubator-specific slope variation.

For a weakly stationary process, every observation has the same finite expectation $\mu$. Linearity of expectation, without any assumption of independence, givesThe variance is a sum of all pairwise covariances:For a real-valued weakly stationary process, the autocovariance is even. There are $n$ pairs at lag zero and $n-h$ pairs in each direction at positive lag $h$. ThereforePositive autocovariances inflate the uncertainty relative to the independent-observation formula; negative terms may reduce it. An unbiased estimator need not have variance tending to zero: a process consisting of the same random variable at every time is weakly stationary, but averaging cannot reduce its variance.

For a known mean zero, a common sample autocovariance function convention isThe weakly stationary process property gives $\mathbb{E}[Y_t{Y}_{t+h}]=\gamma (h)$, soThe lag-zero estimate is unbiased, but at nonzero lags the divisor $n$ causes a finite-sample multiplicative bias. An alternative convention uses $(n-h{)}^{-1}$ instead; for that convention the expectation is exactly $\gamma (h)$. These are two normalizations of the sample autocovariance function, so the divisor must be specified. Here no estimated-mean centering term occurs because zero is known. Replacing zero by a sample mean would introduce additional bias terms. Extend to negative lags by symmetry; there are no observed lagged pairs for $h\ge n$.

A zero-mean autoregressive process of order one satisfieswhere ${\epsilon }_t$ is white noise with variance ${\sigma }_{\epsilon }^2>0$. A causal time series uses present and past innovations only. For $|\alpha |<1$ the serieshas mean-square convergence because ${\sum }_j{\alpha }^{2j}<\infty$, and it supplies the causal weakly stationary process. Its variance is ${\sigma }_{\epsilon }^2/(1-{\alpha }^2)$. For $h\ge 0$, only the common innovations contribute to the autocovariance:Evenness then givesWhen $\alpha <0$, the signs alternate; when $\alpha =0$, the process is white noise. With nondegenerate innovations, $|\alpha |\ge 1$ does not give a causal finite-variance stationary solution. These covariance calculations need uncorrelated innovations, not Gaussianity; Gaussian innovations additionally yield a Gaussian process.

The selected model is a moving-average process of order one, with the plus-sign convention used by R:The Gaussian likelihood function models innovations as independent $N(0,{\sigma }_{\epsilon }^2)$. Its autocovariance is $\gamma (0)={\sigma }_{\epsilon }^2(1+{\theta }^2)$, $\gamma (1)={\sigma }_{\epsilon }^2\theta$, and zero at larger lags, with autocorrelation $\rho (1)=\theta /(1+{\theta }^2)\simeq 0.4631$. The estimated coefficient has modulus less than one, giving invertibility of a moving-average model; a finite moving-average process is stationary regardless of this invertibility restriction.

The original plots show fluctuations around a roughly constant level, no persuasive deterministic trend, and a conspicuous positive lag-one correlation with the other plotted correlations within the approximate sampling bands. They make a stationary nonseasonal working model reasonable for this short series. stationary=TRUE excludes differenced alternatives, and seasonal=FALSE excludes seasonal ARIMA terms. Neither option is proved by a 60-day plot: annual temperature seasonality cannot be excluded from two months, and a longer record or residual diagnostics might require a different model. The options are reasonable local simplifications, not universal statements about temperature dynamics.

Let ${\mathcal{F}}_t$ contain observations through time $t$. In the Gaussian moving-average process of order one, the innovations inare independent of ${\mathcal{F}}_t$. Hence the conditional expectation is the constant mean. For known parameters, the forecast and its mean squared prediction error areThe mean prediction error itself is zero; the requested uncertainty is naturally interpreted as expected squared error. Plugging in the fit gives forecast $15.0133$, mean squared error $6.731002$, and root mean squared error about $2.5944$.

The usual naive time-series forecast repeats $Y_t$. Its unconditional expected squared error issince $\gamma (7)=0$. ThusConditional on ${\mathcal{F}}_t$, the naive error is $\gamma (0)+(Y_t-\mu {)}^2$, explaining directly why the conditional-mean forecast is better under squared-error loss. These are fitted-model calculations that ignore parameter-estimation uncertainty. Estimating the unknown mean from the past adds its estimation mean squared error to the seven-step forecast risk; the printed summary does not provide a complete assessment of that extra uncertainty. This is distinct from claiming a signed mean error as a measure of forecasting accuracy.

There is a literal domain error in the printed formula. A semivariogram indexed by a signed displacement must satisfy $\gamma (h)=\gamma (-h)$, since reversing the two observations leaves the squared increment unchanged. For positive ${\sigma }^2$ or ${\tau }^2$, the displayed positive-lag value is nonzero whereas the displayed negative-lag value is zero. Thus the claim on all of $\mathbb{R}$ is false as written.

Interpret the argument as a nonnegative distance, or use the even extension of the Gaussian semivariogram. For ${\sigma }^2,{\tau }^2\ge 0$ and $\varphi >0$, the intended model isIt admits the covariogramThe Gaussian kernel is a positive-definite kernel, as can be seen from the Fourier transform of a Gaussian: in one dimensionand the nonnegative spectral density makes every finite covariance quadratic form nonnegative. In higher dimensions use the product Gaussian density. Adding independent Gaussian white noise of variance ${\tau }^2$ adds the diagonal nugget covariance. Thus a zero-mean stationary Gaussian random field exists with this covariance, and $C(0)-C(h)$ gives the corrected semivariogram.

For nonzero structured variance, the nugget effect is the right limit ${\tau }^2$ at zero, the sill of a semivariogram is ${\tau }^2+{\sigma }^2$, and the range of a semivariogram is infinite if defined as the distance at which the sill is exactly reached. A 95-percent practical range solves $1-e^{-{\varphi }^2{h}^2}=0.95$:The scale $1/\varphi$ is sometimes called the range parameter; it is neither the exact range nor the 95-percent practical range. If ${\sigma }^2=0$, there is no nontrivial structured range.

The fact that a semivariogram does not determine stationarity is a further distinction: a semivariogram specifies increment variation, not a unique process or absolute covariance. Adding an independent random constant adds a constant to $C$ without changing $\gamma$. More strongly, subtracting $W(0)$ from a stationary field $W(s)$ preserves its semivariogram but generally makes its variance depend on $s$. The corrected statement is that this semivariogram admits a second-order stationary model, with the displayed covariance using the usual convention $C(h)\to 0$ at large distance; it does not force every compatible process to be stationary.

Write $z=(Z_{s_1},\dots ,Z_{s_n}{)}^T$, let ${\Sigma }_{ij}=C(\Vert s_i-s_j\Vert )$, let $c_i=C(\Vert s_i-s_0\Vert )$, and put $v_0=C(0)$. Assume the covariance matrix $\Sigma$ is invertible. Because the mean is known to be zero, every linear predictor $w^{T}z$ is unbiased. Its mean squared prediction error isCompleting the square yieldsThe positive-definite matrix property makes the last term nonnegative, proving the simple kriging formulasThe residual is uncorrelated with $z$. Since the joint law is multivariate normal, it is independent of $z$, so the same predictor and variance describe the conditional distribution. Consequently a 95-percent prediction interval isThere is no constraint that the weights sum to one: that is required for an unknown constant mean in ordinary kriging, not for the stated known-zero-mean problem. If $\Sigma$ is singular, use a Moore-Penrose inverse with the corresponding covariance compatibility condition. At an already observed location prediction of that very same field value has zero error; predicting an independent noisy replicate is a different target. Estimated covariance parameters make this a plug-in interval and generally add uncertainty that the formula does not include.