Next: A Yet Another Shot Up: The Goodness of Fit Previous: Definition of

### Definition of

The big_gamma function, , is given in terms of a continued fraction expansion:

 = = (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) = = fn(x) = = (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:

Now we'll do likewise with f2:

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:

Let us list our obtained expressions for Ai and Bi in order, and attempt some substitutions:

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:

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:

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:

and performing appropriate substitutions:

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:

Matching this against equation 2.31 yields:

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

hence

 (2.33)

Similarly for B1:

Matching this against equation 2.32 yields:

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

hence

 (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 :

we find that

and

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 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:

 assume DO now an becomes a2 now bn becomes b2 now A0 becomes A2 now B0 becomes B2 now an becomes a3 now bn becomes b3 now A1 becomes A3 now B1 becomes B3 UNTIL

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 now A0 becomes A2 now B0 becomes B2 now A1 becomes A3 now B1 becomes B3

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.

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 2.15 shows the Gnuplot command file used to convert try.out into a graph shown in Figure 2.16.

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