Solution (source code)

= Solution

Consider a regular $d$-parameter <statistical model> for which the data can locally be expressed as $(v,a)$, where $v=\widehat\theta$ is the <maximum-likelihood estimator> and $a$ is a suitable exact or higher-order <ancillary statistic>. Write the <log-likelihood> as $\ell(\theta;v,a)$, and let
$$
j(v;v,a)=-\left.\partial_\theta\partial_\theta^{\mathsf T}\ell(\theta;v,a)\right|_{\theta=v}
$$
be its fitted <observed information>. The <P-star approximation> to the conditional sampling density of $v$, given $a$, is
$$
\boxed{p^*(v\mid a;\theta)=c(\theta,a)
|j(v;v,a)|^{1/2}
\exp\{\ell(\theta;v,a)-\ell(v;v,a)\}.}
$$
The normalizing factor is independent of $v$ and is chosen so the density integrates to one over the feasible estimator coordinates. Its leading Gaussian value is $(2\pi)^{-d/2}$, but exact normalization can change it. Both the fitted likelihood and the fitted information vary with $v$; the observed dataset is not held fixed while integrating over this sampling coordinate.

To see the structure, expand the likelihood difference for $v$ close to $\theta$. Its leading term is minus one half of the information-weighted quadratic displacement, while $|j(v)|^{1/2}$ supplies the local density scale. This recovers the first-order <multivariate normal distribution> of an efficient <maximum-likelihood estimator>. The formula is invariant under smooth one-to-one parameter changes: at a likelihood maximum the score vanishes, so the Hessian transforms as a quadratic form. Its square-root determinant supplies exactly the <Jacobian determinant> required for the estimator density to transform.

There is a direct derivation in a full regular <exponential family>. Let $S$ be the <sufficient statistic> for an independent sample, with likelihood
$$
\ell(\theta;s)=\theta^{\mathsf T}s-n\kappa(\theta)+\text{data term}.
$$
The maximum-likelihood equation is $s=n\nabla\kappa(v)$, and its derivative is $j(v)=n\kappa''(v)$. Applying a multivariate <saddlepoint density approximation> to $S$ gives
$$
f_S(s;\theta)\simeq(2\pi)^{-d/2}|j(v)|^{-1/2}
\exp\{\ell(\theta;s)-\ell(v;s)\}.
$$
Changing variables from $s$ to $v$ multiplies by $|ds/dv|=|j(v)|$, producing the displayed P-star form. In more general models, ancillary conditioning supplies the relevant local sample coordinates; an ancillary is not an arbitrary extra statistic.

For an illustrative exact case, take an independent sample from an <exponential distribution> with mean $\theta$, so $v=\bar Y$. The fitted <observed information> is $j(v)=n/v^2$, and
$$
\ell(\theta;v)-\ell(v;v)=n\log(v/\theta)-nv/\theta+n.
$$
Thus P-star has density shape $v^{n-1}\theta^{-n}e^{-nv/\theta}$. Normalizing gives
$$
p^*(v;\theta)=\frac{n^n}{\Gamma(n)\theta^n}v^{n-1}e^{-nv/\theta},\qquad v>0,
$$
which is exactly the <gamma distribution> sampling density of the mean. Here sample ratios provide ancillary coordinates and are independent of the mean.

In suitable regular models the unnormalized leading construction has relative error of order $n^{-1}$, and higher-order ancillary constructions with normalization can attain relative order $n^{-3/2}$. Such accuracy needs the stated smoothness and ancillary hypotheses; it is not a universal consequence of asymptotic normality. The exact scale-model example is likewise a special property, not a general identity. The formula is useful for refined conditional inference, score distributions, and nuisance-likelihood adjustments, including <modified profile likelihood>. It describes a sampling density of the estimator, rather than a posterior density of the parameter.