Dieter Kaiser | 5 Sep 23:56

Implementation of Gamma functions

I have started to implement support for the Incomplete Gamma function. Maxima
only knows the symbols gammaincomplete and gammagreek.

I have already finished the routines for the numerical evaluation of the
Incomplete Gamma function for Float, Complex float and Bigfloat. Now I am
working on the routines for Complex Bigfloats. I have used again an expansion in
a series or continued fractions for the Regularized Incomplete Gamma function. 

In principle it is not necessary to implemenent a new algorithm. We can use the
routines for the Exponential Integral E because gamma(a,z) =
z^a*expintegral_e(1-a,z). But it is interessting to see that both algorithm
converge to every desired accurracy and give correct and equivalent results.

Because the different Gamma functions are interconnected it would be possible to
implement the following functions:

Incomplete Gamma function (The upper tail of the Gamma Incomplete function.)
Complement of the Incomplete Gamma function
Generalized Incomplete Gamma function
Regularized Incomplete Gamma function
Generalized Regularized Incomplete Gamma function

Perhaps the names could be:

gamma_incomplete(a,z)
gamma_greek(a,z)
gen_gamma_incomplete(a,z1,z2)
reg_gamma_incomplete(a,z)
gen_reg_gamma_incomplete(a,z1,z2)

(Continue reading)

Raymond Toy | 6 Sep 06:55

Re: Implementation of Gamma functions



On Fri, Sep 5, 2008 at 5:56 PM, Dieter Kaiser <drdieterkaiser <at> web.de> wrote:
I have started to implement support for the Incomplete Gamma function. Maxima
only knows the symbols gammaincomplete and gammagreek.

I have already finished the routines for the numerical evaluation of the
Incomplete Gamma function for Float, Complex float and Bigfloat. Now I am
working on the routines for Complex Bigfloats. I have used again an expansion in
a series or continued fractions for the Regularized Incomplete Gamma function.
Very cool!
 

Or we could use names like

gamma_incomplete(a,z)
gamma_greek(a,z)                 = 1 - gamma_incomplete(a,z)
gamma_incomplete_gen(a,z1,z2)    = gamma_incomplete(a,z1)-gamma_incomplete(a,z2)
gamma_incomplete_reg(a,z)        = gamma_incomplete(a,z)/gamma(a)
gamma_incomplete_gen_reg(a,z1,z2)= gamma_incomplete_gen(a,z1,z2)/gamma(a)


I prefer these names.  I find incomple and greek rather confusing because I can never remember which one is the tail.

Maybe gamma_incomplete_tail for the integral from z to inf and gamma_incomplete for the integral from 0 to z?



Perhaps more interessting would be to implement a Log Gamma function and the
Inverse Gamma functions.

Well, gamma_incomplete gives us a way to evaluate erf, which we don't have today for complex or bigfloat numbers.

Since we have a gamma function, the log gamma isn't too hard for numerical evaluation.  Inverse gamma is probably not too hard if a Newton iteration can be applied.

Ray
 

_______________________________________________
Maxima mailing list
Maxima <at> math.utexas.edu
http://www.math.utexas.edu/mailman/listinfo/maxima
Mario Rodriguez | 6 Sep 08:25
Favicon

Re: Implementation of Gamma functions

Dieter Kaiser escribió:
> I have started to implement support for the Incomplete Gamma function. Maxima
> only knows the symbols gammaincomplete and gammagreek.
> 
> I have already finished the routines for the numerical evaluation of the
> Incomplete Gamma function for Float, Complex float and Bigfloat. Now I am
> working on the routines for Complex Bigfloats. I have used again an expansion in
> a series or continued fractions for the Regularized Incomplete Gamma function. 
> 

Thanks for working on this. Particularly, I find these functions very 
useful for the distrib package. In contrib/numdistrib.lisp I wrote time 
ago some float based functions for the incomplete gamma and beta, as 
well as for their inverses.

Mario
Robert Dodier | 6 Sep 17:40

Re: Implementation of Gamma functions

On 9/5/08, Dieter Kaiser <drdieterkaiser <at> web.de> wrote:

>  Or we could use names like
>
>  gamma_incomplete(a,z)
>  gamma_greek(a,z)                 = 1 - gamma_incomplete(a,z)
>  gamma_incomplete_gen(a,z1,z2)    = gamma_incomplete(a,z1)-gamma_incomplete(a,z2)
>  gamma_incomplete_reg(a,z)        = gamma_incomplete(a,z)/gamma(a)
>  gamma_incomplete_gen_reg(a,z1,z2)= gamma_incomplete_gen(a,z1,z2)/gamma(a)

Dieter, as usual I have only 1/2 clue about what's going on technically
so I'll restrict myself to some comments about the names.

I'm in favor of big-endian names so the above with gamma_adjective
makes sense to me.

The "greek" bit seems obscure --- surely there is some more widely
recognized word we could substitute for "greek".

I'm usually opposed to cryptic abbreviations in seldom-used functions,
but spelling out "generalized" and "regularized" make for very long names ...
Maybe that's how it should be, I'm undecided. Be that as it may, if an
abbreviation is used, I would suggest "genl" for "generalized".

Log-gamma and gamma inverse functions would be great too,
but there's only so much time in a day, so whenever you want to work
on those, that's terrific.

Thanks for your help,

Robert Dodier
Kostas Oikonomou | 6 Sep 18:59
Favicon

Re: Implementation of Gamma functions

Dieter,

Thanks for the great work on gamma functions.

I agree with Ray and Robert that the gamma_* names are 
preferable.

And I also find gamma_greek hard to remember/understand. 
How about gamma_incomplete_compl()?

				Kostas

Robert Dodier wrote:
> On 9/5/08, Dieter Kaiser <drdieterkaiser <at> web.de> wrote:
> 
>>  Or we could use names like
>>
>>  gamma_incomplete(a,z)
>>  gamma_greek(a,z)                 = 1 - gamma_incomplete(a,z)
>>  gamma_incomplete_gen(a,z1,z2)    = gamma_incomplete(a,z1)-gamma_incomplete(a,z2)
>>  gamma_incomplete_reg(a,z)        = gamma_incomplete(a,z)/gamma(a)
>>  gamma_incomplete_gen_reg(a,z1,z2)= gamma_incomplete_gen(a,z1,z2)/gamma(a)
> 
> Dieter, as usual I have only 1/2 clue about what's going on technically
> so I'll restrict myself to some comments about the names.
> 
> I'm in favor of big-endian names so the above with gamma_adjective
> makes sense to me.
> 
> The "greek" bit seems obscure --- surely there is some more widely
> recognized word we could substitute for "greek".
> 
> I'm usually opposed to cryptic abbreviations in seldom-used functions,
> but spelling out "generalized" and "regularized" make for very long names ...
> Maybe that's how it should be, I'm undecided. Be that as it may, if an
> abbreviation is used, I would suggest "genl" for "generalized".
> 
> Log-gamma and gamma inverse functions would be great too,
> but there's only so much time in a day, so whenever you want to work
> on those, that's terrific.
> 
> Thanks for your help,
> 
> Robert Dodier
> _______________________________________________
> Maxima mailing list
> Maxima <at> math.utexas.edu
> http://www.math.utexas.edu/mailman/listinfo/maxima
Dieter Kaiser | 6 Sep 21:55

Re: Implementation of Gamma functions

Thank very much for the feedback.

I like long names too. For use it is possible to define a short name via a
function definition or an alias. Second, I would like the convention that the
main name of the function is always at the beginning e.g. laguerre_generalized
and not gen_laguerre. The long names would look like:

gamma_incomplete(a,z)
gamma_incomplete_generalized(a,z1,z2)
gamma_incomplete_regularized(a,z)    
gamma_incomplete_generalized_regularized(a,z1,z2)

Yes, gamma_greek is the most confusing name and perhaps we should avoid the use
of this function. The Maxima symbol gammagreek is only used in the code of
specint and perhaps could be completely eliminated.

The biggest problem is that it is not always clear what is meant by the
Incomplete Gamma function. A&S gives three different definitions with the
formulas 6.5.1, 6.5.2 and 6.5.3 but do not distinguish the definitions with
different names. The definitions are:

6.5.1: P(a,x)     = 1/Gamma(a) * integral(exp(-t)*t^(a-1),t,0,x) (lower tail)

6.5.2: gamma(a,x) = P(a,x)*Gamma(a)                  (lower case Greek gamma)

6.5.3: Gamma(a,x) = integral(exp(-t)*t^(a-1),t,x,inf)            (upper tail)
                  = Gamma(a)-gamma(a,x)              (upper case Greek gamma)

The Maxima symbol gammagreek follows the definition in A&S 6.5.2. The Maxima
symbol gammaincomplete means the definition 6.5.3.

The Numerical Recipes use the following conventions:

P(a,x) = Gamma(a,x)/gamma(a)   (Incomplete Gamma function)
Q(a,x) = gamma(a,x)/gamma(a)   (Complement of the Incomplete Gamma function)
       = 1-Gamma(a,x)/Gamma(a)

That is opposite to A&S. Here P(a,x) is the upper tail and not the lower tail of
the integral and the function is divided by Gamma(a). The numerical Recipes do
not give extra names for Gamma(a,x) or gamma(a,x).

The names I have suggested follow the conventions of functions.wolfram.com.
There we have the following definition for the Incomplete Gamma function which
is equal to A&S 6.5.3:

Gamma(a,z) = (exp(-t)*t^(a-1),t,x,inf)   (upper tail, Upper case Greek gamma)

A very important feature of this definition is that the interconnections between
a lot of special functions can be expressed more easy and directly with this
definition of the Incomplete Gamma function. 

functions.wolfram.com do not introduce a name for the complement of the
Incomplete Gamma function. For the Gamma Incomplete function divided by gamma(a)
the name Regularized Incomplete Gamma function is introduced and the symbol
Q(a,x).

Sorry, but I have not seen that we have numerical routines for the Gamma
functions in the code of numdistrib.lisp. I will have a look at the routines and
compare them against the routines I have written. 

Dieter Kaiser
Kostas Oikonomou | 7 Sep 01:22
Favicon

Re: Implementation of Gamma functions

> gamma_incomplete(a,z)
> gamma_incomplete_generalized(a,z1,z2)
> gamma_incomplete_regularized(a,z)    
> gamma_incomplete_generalized_regularized(a,z1,z2)

I think these are good choices.   Having used Mathematica 
for a long time, I have experienced the usefulness of long, 
descriptive names for functions or other entities that are 
not frequently used.
And having a convention for how to order the parts of a long 
name is also very helpful to the user.

> Yes, gamma_greek is the most confusing name and perhaps we should avoid the use
> of this function. The Maxima symbol gammagreek is only used in the code of
> specint and perhaps could be completely eliminated.

That would be ok by me.

> The biggest problem is that it is not always clear what is meant by the
> Incomplete Gamma function. A&S gives three different definitions with the
> formulas 6.5.1, 6.5.2 and 6.5.3 but do not distinguish the definitions with
> different names. 
> The names I have suggested follow the conventions of functions.wolfram.com.
> There we have the following definition for the Incomplete Gamma function which
> is equal to A&S 6.5.3:
> 
> Gamma(a,z) = (exp(-t)*t^(a-1),t,x,inf)   (upper tail, Upper case Greek gamma)
> 
> A very important feature of this definition is that the interconnections between
> a lot of special functions can be expressed more easy and directly with this
> definition of the Incomplete Gamma function. 

I'm used to the Wolfram/Mathematica definitions, even though 
they don't fully agree with A&S.
I have also found the regularized versions quite useful, and 
I'm glad to see you including them in Maxima.

						Kostas

Gmane