= Solution
The <P-star approximation> describes the sampling <density> of a regular <maximum-likelihood estimator> in sample coordinates consisting of the fitted parameter $\widehat\theta$ and a suitable <ancillary statistic> $a$. For a $d$-dimensional parameter, write $\ell(\theta;\widehat\theta,a)$ for the <log-likelihood> expressed in these coordinates, and $j(\widehat\theta;\widehat\theta,a)$ for the fitted <observed information>. The approximation is
$$
p^*(\widehat\theta\mid a;\theta)=c(\theta,a)|j(\widehat\theta;\widehat\theta,a)|^{1/2}\exp\{\ell(\theta;\widehat\theta,a)-\ell(\widehat\theta;\widehat\theta,a)\}.
$$
The <determinant> factor supplies the local volume scale, while the exponential retains the full <likelihood ratio>, rather than replacing it by a quadratic. The <normalizing constant> is chosen by integrating over the fitted-parameter coordinates; its leading regular large-sample value is $(2\pi)^{-d/2}$. These ingredients arise from a <statistical saddlepoint approximation>. Their use requires an interior, locally unique fit, nonsingular information and suitable smoothness; an <ancillary statistic> specifies the conditioning surface when the fit alone does not exhaust the relevant sample coordinates. In a regular model the unnormalized leading expression has relative error of order $n^{-1}$ locally; normalization and conditioning matter for refined inference. It is not a universal exact <density>. Under a smooth one-to-one parameter transformation, the <Hessian> at the fit transforms without a score term, so its <determinant> supplies precisely the <density> <Jacobian determinant>. This proves the invariance of the construction under one-to-one changes of <statistical parameters>.
Here put $S=\sum_iX_i$ and $T=\sum_iY_i$. The <log-likelihood>, apart from a data-only constant, is
$$
\ell(\psi,\lambda)=n\log\psi+2n\log\lambda-\lambda S-\psi\lambda T.
$$
The two <score equations> give
$$
\boxed{\widehat\lambda=\frac nS,\qquad\widehat\psi=\frac ST.}
$$
The sample sums form a <sufficient statistic>, and there is no need to retain an additional ancillary in calculating their joint <density>. Inverting the fitted coordinates gives $S=n/\widehat\lambda$ and $T=n/(\widehat\psi\widehat\lambda)$. With parameter order $(\psi,\lambda)$, the fitted <observed information> is
$$
j(\widehat\psi,\widehat\lambda)=\begin{pmatrix}n/\widehat\psi^2&n/(\widehat\psi\widehat\lambda)\\n/(\widehat\psi\widehat\lambda)&2n/\widehat\lambda^2\end{pmatrix},\qquad |j|^{1/2}=\frac n{\widehat\psi\widehat\lambda}.
$$
Also,
$$
\ell(\psi,\lambda)-\ell(\widehat\psi,\widehat\lambda)=n\log\frac\psi{\widehat\psi}+2n\log\frac\lambda{\widehat\lambda}-\frac{n\lambda}{\widehat\lambda}\left(1+\frac\psi{\widehat\psi}\right)+2n.
$$
Thus the required joint <P-star approximation> is
$$
\boxed{p^*(\widehat\psi,\widehat\lambda)=\frac{c_n n e^{2n}\psi^n\lambda^{2n}}{\widehat\psi^{n+1}\widehat\lambda^{2n+1}}\exp\left\{-\frac{n\lambda}{\widehat\lambda}\left(1+\frac\psi{\widehat\psi}\right)\right\},\quad \widehat\psi,\widehat\lambda>0.}
$$
To integrate out $\widehat\lambda$, substitute $u=b/\widehat\lambda$ in
$$
\int_0^\infty v^{-2n-1}e^{-b/v}\,dv=\frac{\Gamma(2n)}{b^{2n}}.
$$
Taking $b=n\lambda(1+\psi/\widehat\psi)$ and rearranging the powers yields
$$
p^*(\widehat\psi)=\frac{c_ne^{2n}n^{1-2n}\Gamma(2n)}\psi\left(\frac{\widehat\psi}\psi\right)^{n-1}\left(1+\frac{\widehat\psi}\psi\right)^{-2n}.
$$
The <beta function> identity $\int_0^\infty z^{n-1}(1+z)^{-2n}\,dz=\Gamma(n)^2/\Gamma(2n)$ therefore gives
$$
\boxed{c_n=\frac{n^{2n-1}e^{-2n}}{\Gamma(n)^2},\qquad p^*(\widehat\psi)=\frac{\Gamma(2n)}{\Gamma(n)^2\psi}\left(\frac{\widehat\psi}\psi\right)^{n-1}\left(1+\frac{\widehat\psi}\psi\right)^{-2n}.}
$$
Both normalized approximations are actually exact. To check the joint assertion independently, $S$ and $T$ have independent <gamma distributions> with shape $n$ and rates $\lambda$ and $\psi\lambda$. Their joint <density> is $\psi^n\lambda^{2n}S^{n-1}T^{n-1}e^{-\lambda S-\psi\lambda T}/\Gamma(n)^2$. The inverse-coordinate <Jacobian determinant> is
$$
\left|\frac{\partial(S,T)}{\partial(\widehat\psi,\widehat\lambda)}\right|=\frac{n^2}{\widehat\psi^2\widehat\lambda^3},
$$
which reproduces the normalized joint expression. This is the <P-star density for two exponential rates>; its marginal gives the stated <F-distribution>. <Stirling's approximation> gives $c_n=(2\pi)^{-1}(1-1/(6n)+O(n^{-2}))$, so using only the leading normalizer retains the exact shape but not its exact total mass.
The original PDF's displayed <density> constant is missing a square on $\Gamma(n)$. Its printed expression integrates to $\Gamma(n)$, rather than one in general. The normalized <density> above is the one consistent with the stated <F-distribution>.
Back to article page
