next up previous index
Next: A Yet Another Shot Up: The Goodness of Fit Previous: Definition of

   
Definition of $\Gamma (a,x)$

The big_gamma function, $\Gamma (a,x)$, is given in terms of a continued fraction expansion:

$\displaystyle \Gamma(a, x)$ = $\displaystyle e^{-x}x^a
\frac{1}
{x + \frac{1 - a}
{1 + \frac{1}
{x + \frac{2-a...
... + \frac{2}
{x + \frac{3-a}
{1 + \frac{3}
{x + \frac{4 - a}
{1 + \ldots}}}}}}}}$  
  = $\displaystyle e^{-x}x^a
\left(\frac{1}{x+}\, \frac{1-a}{1+} \, \frac{1}{x+} \, ...
...-a}{1+}
\frac{2}{x+}\,
\cdots \, \frac{n-a}{1+} \, \frac{n}{x+}\, \cdots\right)$ (2.27)

Continued fraction expansions may look dreadful but they converge extremely rapidly compared to power series expansions. For this reason they've been of great interest to numerical analysts and, of course, to computational scientists too. In case of the incomplete gamma function the continued fraction expansion converges best where the series expansion does not, so there is also a very nice complementarity here.

Unlike a power series expansion the continued fraction expansion seemingly cannot be evaluated in steps from ``left to right'', i.e., from the lowest terms, gradually adding higher and higher terms to the sum. Here it looks like we have to start from the end, sic! That is, we have to guess the last term, and then gradually build up the whole fraction starting from the last term in the expansion and moving towards the beginning.

But it can be done the right way, i.e., from ``left to right''. Consider the following generic continued fraction expansions:

fn-1(x) = $\displaystyle b_0 + \frac{a_1}{b_1 +}\,\frac{a_2}{b_2 +}\,
\cdots \, \frac{a_{n-2}}{b_{n-2}+}\,
\frac{a_{n-1}}{b_{n-1}}$  
  = $\displaystyle \frac{A_{n-1}}{B_{n-1}}$  
fn(x) = $\displaystyle b_0 + \frac{a_1}{b_1 +}\,\frac{a_2}{b_2 +}\,
\cdots \, \frac{a_{n-1}}{b_{n-1}+}\,
\frac{a_{n}}{b_{n}}$  
  = $\displaystyle \frac{A_n}{B_n}$ (2.28)

The question is, how An and Bn relate to An-1, Bn-1, An-2, Bn-2, etc.

We can help ourselves with Maxima   in this task. First we define f1 and extract its numerator and denominator:

\begin{figure}
\begin{footnotesize}
\begin{tex2html_preform}\begin{verbatim}(C12...
.../R/ B
1
(C15)\end{verbatim}\end{tex2html_preform}\end{footnotesize}\end{figure}

Now we'll do likewise with f2:
\begin{figure}
\begin{footnotesize}
\begin{tex2html_preform}\begin{verbatim}(C15...
...A
1 2 2
(C18)\end{verbatim}\end{tex2html_preform}\end{footnotesize}\end{figure}

Observe that without asking Maxima has grouped biga[2] in a way that shows the presence of biga[1] in it! Let us now evaluate f3:
\begin{figure}
\begin{footnotesize}
\begin{tex2html_preform}\begin{verbatim}(C18...
... 2 3 1 3
(C21)\end{verbatim}\end{tex2html_preform}\end{footnotesize}\end{figure}

Let us list our obtained expressions for Ai and Bi in order, and attempt some substitutions:
\begin{figure}
\begin{footnotesize}
\begin{tex2html_preform}\begin{verbatim}(C30...
... 2 3 1 3
(C36)\end{verbatim}\end{tex2html_preform}\end{footnotesize}\end{figure}

In order to carry out substitutions we use Maxima  function ratsubst, which takes three arguments: the name of a new variable, its value, and the expression into which the variable should be substituted:
\begin{figure}
\begin{footnotesize}
\begin{tex2html_preform}\begin{verbatim}(C36...
...1 A
3 3
(C40)\end{verbatim}\end{tex2html_preform}\end{footnotesize}\end{figure}

What we seem to be getting is the following relation:

A3 = b3 A2 + a3 A1 (2.29)

Now let us have a closer look at Bi:
\begin{figure}
\begin{footnotesize}
\begin{tex2html_preform}\begin{verbatim}(C40...
...1 A
3 3
(C42)\end{verbatim}\end{tex2html_preform}\end{footnotesize}\end{figure}

Here we got that

B3 = b3 B2 + a3 B1 (2.30)

Could it be that the general formula is:
  
An = bn An-1 + an An-2 (2.31)
Bn = bn Bn-1 + an Bn-2 (2.32)

We can check it by evaluating f4:
\begin{figure}
\begin{footnotesize}
\begin{tex2html_preform}\begin{verbatim}(C42...
... 1 2 2 4
(C45)\end{verbatim}\end{tex2html_preform}\end{footnotesize}\end{figure}

and performing appropriate substitutions:
\begin{figure}
\begin{footnotesize}
\begin{tex2html_preform}\begin{verbatim}(C45...
...2 A
4 4
(C49)\end{verbatim}\end{tex2html_preform}\end{footnotesize}\end{figure}

Well, formulas 2.31 and 2.32 clearly hold. But the recurrence relationships for An and Bn are identical, so clearly, if they are not to be the same their start-up values must be different. What are they? Let us have a look at A1 first:
\begin{figure}
\begin{footnotesize}
\begin{tex2html_preform}\begin{verbatim}(C52...
...A
0 1 1
(C53)\end{verbatim}\end{tex2html_preform}\end{footnotesize}\end{figure}

Matching this against equation 2.31 yields:

A1 = b1 A0 + a1 A-1 = b1 b0 + a1 1

hence

 \begin{displaymath}
A_0 = b_0\qquad\textrm{and}\qquad A_{-1} = 1
\end{displaymath} (2.33)

Similarly for B1:
\begin{figure}
\begin{footnotesize}
\begin{tex2html_preform}\begin{verbatim}(C53...
.../R/ B
1
(C54)\end{verbatim}\end{tex2html_preform}\end{footnotesize}\end{figure}

Matching this against equation 2.32 yields:

B1 = b1 B0 + a1 B-1 = b1 1 + a1 0

hence

 \begin{displaymath}
B_0 = 1\qquad\textrm{and}\qquad B_{-1} = 0
\end{displaymath} (2.34)

The evaluation of a continued fraction expansion now looks much easier. All that we have to do is to

1.
assign values to A-1, A0, B-1, and B0, using equations 2.33 and 2.34
2.
evaluate A1 and B1 from the recurrence relations 2.31 and 2.32 and find A1/B1,
3.
evaluate A2 and B2 likewise, find A2/B2, and check for convergence, if none
4.
evaluate A3 and B3, as before, find A3/B3, and check for convergence,
5.
and so on until the result converges, i.e., until the result no longer changes, within a prescribed accuracy, as we embark on a yet another iteration.

Looking once again at the continued fraction expansion for $\Gamma (a,x)$:

\begin{displaymath}e^{-x}x^a
\left(\frac{1}{x+}\, \frac{1-a}{1+} \, \frac{1}{x+...
...}\,
\cdots \, \frac{n-a}{1+} \, \frac{n}{x+}\, \cdots\right)
\end{displaymath}

we find that

\begin{eqnarray*}A_{-1} &=& 1 \\
A_0 &=& b_0 = 0 \\
B_{-1} &=& 0 \\
B_0 &=& 1
\end{eqnarray*}


and

\begin{eqnarray*}a_1 &= 1 &\Rightarrow A_1 = b_1 A_0 + a_1 A_{-1} = x\cdot 0 + 1...
... n - a &\\
b_{2n} &= 1 &\\
a_{2n+1} &= n &\\
b_{2n+1} &= x
\end{eqnarray*}


Observe that a regular sequencing begins only with a2 and b2, not with a1 and b1. In fact the first term in the continued expansion for $\Gamma (a,x)$ is the odd one. Consequently we shall begin our iterations from evaluating A2 and B2, and we'll use A0, B0, A1, and B1 as our startup values.

Using mathematicised notation our algorithm is going to be as follows:

$ A_0 \leftarrow 0 $  
$ B_0 \leftarrow 1 $  
$ A_1 \leftarrow 1 $  
$ B_1 \leftarrow x $  
$ R_{\textit{\scriptsize old}} \leftarrow A_1 / B_1 $ assume $x \neq 0$
$ n \leftarrow 1 $  
DO  
$ \qquad a_n \leftarrow n - a $ now an becomes a2
$ \qquad b_n \leftarrow 1 $ now bn becomes b2
$ \qquad A_0 \leftarrow b_n A_1 + a_n A_0 $ now A0 becomes A2
$ \qquad B_0 \leftarrow b_n B_1 + a_n B_0 $ now B0 becomes B2
$ \qquad a_n \leftarrow n $ now an becomes a3
$ \qquad b_n \leftarrow x $ now bn becomes b3
$ \qquad A_1 \leftarrow b_n A_0 + a_n A_1 $ now A1 becomes A3
$ \qquad B_1 \leftarrow b_n B_0 + a_n B_1 $ now B1 becomes B3
$ \qquad R_{\textit{\scriptsize new}} \leftarrow A_1/B_1 $  
$ \qquad \delta \leftarrow R_{\textit{\scriptsize new}} - R_{\textit{\scriptsize old}} $  
$ \qquad R_{\textit{\scriptsize old}} \leftarrow R_{\textit{\scriptsize new}} $  
$ \qquad n \leftarrow n+1 $  
UNTIL $ \left\vert\delta/R_{\textit{\scriptsize new}}\right\vert < \epsilon $  

At this stage we have carefully mapped the algorithm, without attempting any economies. But having a closer look at the DO loop as it's been outlined above shows that we can make some savings. First, we don't really need to maintain the variables an and bn. Instead we can use n-a, 1, n and x directly:

DO  
$ \qquad A_0 \leftarrow A_1 + (n - a) A_0 $ now A0 becomes A2
$ \qquad B_0 \leftarrow B_1 + (n - a) B_0 $ now B0 becomes B2
$ \qquad A_1 \leftarrow x A_0 + n A_1 $ now A1 becomes A3
$ \qquad B_1 \leftarrow x B_0 + n B_1 $ now B1 becomes B3
$ \qquad \cdots $  

This saves us 4 assignments within the loop. Assignments are not as cheap as they seem to be, because data movement is involved

Figure 2.13 shows subprogram big_gamma. You should be able to identify clearly all elements. The algorithm is exactly as specified above.


  
Figure: Subprogram big_gamma, which implements function $\Gamma (a,x)$
\begin{figure}\begin{footnotesize}
\begin{tex2html_preform}\begin{verbatim}real(...
...tion big_gamma\end{verbatim}\end{tex2html_preform}\end{footnotesize}\end{figure}

Having now implemented all required functions we can not only calculate the goodness of fit in our problem, which happens to be Q=0.053, but we can test the implementation of function goodness itself.

Figure 2.14 shows a short program, which calculates 1 - Q(a, x) = P(a, x), and writes the results on a data file, try.out. This file can then be displayed using Gnuplot and the graph compared against similar graphs found in the literature.


  
Figure: Program try, which evaluates P(a, x) = 1 - Q(a, x) for four different values of a: 0.5, 1.0, 3.0, and 10.0, and for $x \in [0, 15.0]$.
\begin{figure}\begin{footnotesize}
\begin{tex2html_preform}\begin{verbatim}progr...
...nd program try\end{verbatim}\end{tex2html_preform}\end{footnotesize}\end{figure}

Figure 2.15 shows the Gnuplot command file used to convert try.out into a graph shown in Figure 2.16.


  
Figure 2.15: Gnuplot command file used to plot graph shown in Figure 2.16.
\begin{figure}\begin{footnotesize}
\begin{tex2html_preform}\begin{verbatim}set x...
...''
replot
quit\end{verbatim}\end{tex2html_preform}\end{footnotesize}\end{figure}


  
Figure 2.16: Function P(a, x) = 1 - Q(a, x) for a = 0.5, 1.0, 3.0, 10.0.
\includegraphics[width=10cm]{try.eps}


next up previous index
Next: A Yet Another Shot Up: The Goodness of Fit Previous: Definition of
Zdzislaw Meglicki
2001-02-26