Solving for endpoints given by two integral equations

I have been trying to work out an example related to hypothesis testing for the scale parameter in an exponential distribution. By following the statistic theory I have been led to the following mathematics question.

For a fixed natural number $n$ is there a way to calculate constants $a, b >0$ that satisfy the following to equations

$$\int_a^b x^{n-1}e^{-x} dx = 0.95 \times (n-1)! ~~\text{and}~~ \int_a^b x^n e^{-x}dx = 0.95 \times n!.$$

I think I can calculate $a$ and $b$ by sampling from a Gamma distribution with scale $1$ and shape $n$ but is there a way more directly get a handle on $a$ and $b$?

Solution 1:

Using algebra.

Writing the equations as $$f_1(n)=\log \Big[\Gamma (n,a)-\Gamma (n,b)\Big]-\log \left(\frac{19}{20} (n-1)!\right)$$ $$f_2(n)=\log \Big[\Gamma (n+1,a)-\Gamma (n+1,b)\Big]-\log \left(\frac{19}{20} n!\right)$$ and defining the norm $$\Phi(a,b)=f_1^2(n)+f_2^2(n)$$ the minimization of $\Phi(a,b)$ doe not make much problems since can easily be built the analytical Jacobian and Hessian.

Below are given some results

$$\left( \begin{array}{cc} n & a & b \\ 1 & 0.04236 & 4.76517 \\ 2 & 0.30350 & 6.40122 \\ 3 & 0.71250 & 7.94830 \\ 4 & 1.20696 & 9.43022 \\ 5 & 1.75808 & 10.8644 \\ 6 & 2.35023 & 12.2623 \\ 7 & 2.97386 & 13.6315 \\ 8 & 3.62263 & 14.9773 \\ 9 & 4.29208 & 16.3036 \\ 10 & 4.97893 & 17.6134 \\ 11 & 5.68069 & 18.9088 \\ 12 & 6.39540 & 20.1917 \\ 13 & 7.12151 & 21.4636 \\ 14 & 7.85776 & 22.7257 \\ 15 & 8.60308 & 23.9789 \\ 16 & 9.35661 & 25.2240 \\ 17 & 10.1176 & 26.4619 \\ 18 & 10.8854 & 27.6931 \\ 19 & 11.6595 & 28.9181 \\ 20 & 12.4393 & 30.1374 \\ \end{array} \right)$$

The problem could also be reduced to a single variable. Using a now deleted answer from @PierreCarre, the two parameters are related by $$e^{-a}\, a^n = e^{-b}\, b^n \implies b=-n W_{-1}\left(-\frac{a }{n}e^{-\frac{a}{n}}\right)$$ where $W_{-1}$ is the second branch of Lambert function. This reduces the problem to one equation for one unknown; this equation can easily be solve using Newton method.

For $20 \leq n \leq 100$, a quick and dirty nonlinear regression $(R^2 > 0.999999)$ gives for $a$ an estimate $$a =\alpha + \beta \,n- \gamma \log(n)$$

$$\begin{array}{clclclclc} \text{} & \text{Estimate} & \text{Standard Error} & \text{Confidence Interval} \\ \alpha & 4.06273 & 0.07949 & \{3.90441,4.22104\} \\ \beta & 0.93189 & 0.00051 & \{0.93087,0.93290\} \\ \gamma & 3.39745 & 0.02724 & \{3.34319,3.45171\} \\ \end{array}$$

Tested for $n=200$, the estimate is $a=172.440$ while the solution is $a=173.537$.