[Date Index]
[Thread Index]
[Author Index]
Re: LogErfc
*To*: mathgroup at smc.vnet.net
*Subject*: [mg122650] Re: LogErfc
*From*: Peter Falloon <pfalloon at gmail.com>
*Date*: Fri, 4 Nov 2011 06:02:36 -0500 (EST)
*Delivered-to*: l-mathgroup@mail-archive0.wolfram.com
*References*: <j8tkfh$eu0$1@smc.vnet.net>
On Nov 3, 7:47 pm, Shawn Garbett <shawn.garb... at gmail.com> wrote:
> I need to have a LogErfc function to be able to compute another
> function in a numerically stable manner. I found the function in the
> GSL, I found it in R, and even in another system, there's a copy of code to do
> it floating around. However, I need it in Mathematica and Log[Erfc[x]]
> doesn't cut in due to underflow. Cody solved this in 1969, and I could
> code his solution, but I was hoping that someone had already solved
> this in Mathematica. Anyone?
As you probably know, the problem here is that Erfc[x] goes to zero
like Exp[-z^2], and so quickly reaches an underflow, whereas
Log[Erfc[x]] is still a number of tractable magnitude.
One way around this is to find a representation in terms of another
special function in which the exponential dependence is explicit: that
way, with a bit of simple manipulation you can apply the Log function
and extract the troublesome exponential dependence. One such
representation is in terms of the hypergeometric U function (http://
functions.wolfram.com/06.27.26.0002.01). For real z>0 the annoying (z/
Sqrt[z^2]) factor in that equation reduces to 1, and you get:
Erfc[z] == Exp[-z^2] * HypergeometricU[1/2, 1/2, z^2] / Sqrt[Pi]
Unfortunately, in Mathematica HypergeometricU[1/2,1/2,z^2] is
converted to Exp[z^2]*Gamma[1/2, z^2] internally before computing, so
we're back where we started (since Gamma[1/2,z^2] suffers the same
Underflow[] issues as Erfc for large z). However, if you use the trick
of converting one of the 1/2 factors to finite precision (i.e. Real
instead of Integer) then a different algorithm is used, which doesn't
involve the pesky exponential term.
So to cut to the chase, here is a simple definition which should solve
your problem:
Clear[logErfc]
logErfc[z_?NumericQ /; Precision[z]<Infinity] := -z^2 - (1/2)Log[Pi] +
Log[HypergeometricU[SetPrecision[1/2, Precision[z]], 1/2, z^2]]
Cheers,
Peter.
Prev by Date:
**Re: How to eliminate noises? A better way perhaps.**
Next by Date:
**Re: nVidia Optumus prevents using CUDA?**
Previous by thread:
**Re: LogErfc**
Next by thread:
**Sending an email with a condition be verified**
| |