P-star density for two exponential rates (source code)

= P-star density for two exponential rates
{c}
{title2=$p^*(\widehat\psi,\widehat\lambda)$}

For independent samples of size $n$ from <exponential distributions> with rates $\lambda$ and $\psi\lambda$, the <maximum-likelihood estimators> are $\widehat\lambda=n/S$ and $\widehat\psi=S/T$, where $S$ and $T$ are the respective sample sums. Their normalized <P-star approximation> is exact:
$$
p(\widehat\psi,\widehat\lambda;\psi,\lambda)=\frac{n^{2n}\psi^n\lambda^{2n}}{\Gamma(n)^2\widehat\psi^{n+1}\widehat\lambda^{2n+1}}\exp\left\{-\frac{n\lambda}{\widehat\lambda}\left(1+\frac\psi{\widehat\psi}\right)\right\},\qquad \widehat\psi,\widehat\lambda>0.
$$
Indeed, $S$ and $T$ have independent <gamma distributions> with shapes $n$ and respective rates $\lambda,\psi\lambda$. The <Jacobian determinant> of their inverse transformation is $n^2/(\widehat\psi^2\widehat\lambda^3)$, giving the displayed <density>. Alternatively the square root of the determinant of the fitted <observed information> is $n/(\widehat\psi\widehat\lambda)$; multiplying by the exponential of the fitted <log-likelihood> difference gives the same shape. Its normalizing constant in the <P-star approximation> is $c_n=n^{2n-1}e^{-2n}/\Gamma(n)^2$. Integrating over $\widehat\lambda$ gives
$$
p(\widehat\psi;\psi)=\frac{\Gamma(2n)}{\Gamma(n)^2\psi}\left(\frac{\widehat\psi}\psi\right)^{n-1}\left(1+\frac{\widehat\psi}\psi\right)^{-2n},
$$
so $\widehat\psi/\psi$ has the <F-distribution> with degrees of freedom $(2n,2n)$.