[R] precision of gamma function

[Mike Marchywka]

----------------------------------------
> Date: Wed, 9 Feb 2011 18:23:21 +0100
> From: savicky at praha1.ff.cuni.cz
> To: r-help at r-project.org
> Subject: Re: [R] precision of gamma function
>
> On Wed, Feb 09, 2011 at 05:23:27PM +0200, Chuse chuse wrote:
> > Dear R users,
> >
> > I have to calculate gamma functions for negative numbers beyond -171.4.
> > e.x. gamma(-500.4)
> > I got following:
> >
> > > gamma(-170.4)
> > [1] -5.824625e-308
> > > gamma(-171.4)
> > [1] 0
> > Warning message:
> > underflow occurred in 'gammafn'
> >
> > I have tried to use a recursion getting values a little futher -180.
> > How could I solve this problem? Thank you beforehand.
>
> Hello.
>
> As others pointed out, the function lgamma() should be used. If
> the accuracy of double precision is not sufficient or in case
> of doubt concerning accuracy in extreme cases, Rmpfr package may
> be used. In this case, we get


Unless you expect your end result to be gazillions, it may be worth
considering algrebra. The real question may be,
"What does R support for simplifaction or symbolic manipulation?"



>
> x <- c(-171.4, -170.4)
> y1 <- lgamma(x)
> y1
>
> [1] -712.5781 -707.4341
>
> library(Rmpfr)
> y2 <- lgamma(mpfr(x, precBits=200))
> y1 - y2
>
> 2 'mpfr' numbers of precision 200 bits
> [1] -8.0100329934347335532420397139357516023505624051256280109870069e-14
> [2] -1.0849741205998334826524274950211713854549088941953001164911718e-13
>
> Hope this helps.
>
> Petr Savicky.
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.