Mathematica 9 is now available
Services & Resources / Wolfram Forums / MathGroup Archive
-----

MathGroup Archive 2011

[Date Index] [Thread Index] [Author Index]

Search the Archive

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