On Thursday, April 3, 2014 11:35:10 AM UTC+2, tho.s...@gmail.com wrote:
> Hello,
>
> for my MCMC fitting program, I need to evaluate functions of the form (Gaussian with a one sided exponential tail towards lower x-values):
>
>
>
> f(a,b,c,d) * erfc( g(a,b,c,d) ) * exp( h(a,b,c,d) ) := X * Y * Z = F
>
>
>
> where f,g and h are certain functions of the parameters a,b,c and d.
>
>
>
> It almost always happens that the numbers of these three factors are like:
>
>
>
> F = X * Y * Z = 1e2 * 1e-999 * 1e1000 = 1e3
>
>
>
> Which is a big problem since 1e-999 is represented as 0 and 1e1000 is represented as infinity, thus the result being 0, infinity or nan, but definetly not 1e3.
>
> As a work-around, I went to log-space such that:
>
>
>
> F = exp( ln(F) ) = exp( ln(X * Y * Z) ) = exp( ln(X) + ln(Y) + ln(Z) ) =
>
> = exp( ln(f(a,b,c,d)) + ln(erfc(g(a,b,c,d))) + ln(exp(h(a,b,c,d))) ) :=
>
> := exp( Q + W + E )
>
>
>
> Q and E are no problem to evaluate since f() is just a rational function and ln(exp(h())) is just h().
>
> However, W = ln(erfc(g())) contains the same problem as above. If g() is far negative from 0, erfc(g()) is just 2 (and not e.g. 2 - 1e-99). If g() is far positive from 0, erfc(g()) is just 0, returning W as -Inf (as erfc(g()) should actually be something like 1e-99).
>
>
>
> Now, I looked up several representations of the erfc() function in order to build something like a lnerfc - function. I have chosen the erfcc() function in Numerical recipes, Chapter 6, Special Functions (around page 214) which is also given in Wikipedia at http://en.wikipedia.org/wiki/Error_function#Numerical_approx imation
>
> This approximation has two major advantages:
>
> 1) It is represented as proprotional to an exponential function, for which the ln can easily be calculated.
>
> 2) The fractional error is "everywhere less than 1.2e-7".
>
>
>
> Including all these work-arounds, F = X * Y * Z can be calculated to a good enough precision (for me).
>
>
>
> However (again), as you might already think of, it takes a while to calculate F. In a MCMC run, this function has to be evaluated over and over again. If there is more than one such a function present in my data (say N), I need to fit, i.e. evaluate something like:
>
>
>
> sum(F_i, i=0..N)
>
>
>
> over and over again (typically N = 20..30).
>
>
>
> To put it in a nutshell:
>
> I am looking for a speed-up to calculate W = ln(erfc(g(a,b,c,d))).
>
> I know that I can calculate the erfc - function by:
>
> erfc(x) = 1 - sgn(x) * igamma(0.5,x^2)
>
> where igamma is the incomplete gamma-function.
>
> Unfortunately, there is no LNIGAMMA - function in IDL, as for the complete gamma-function (LNGAMMA). As this does not necessarily have to work good then because of the "1 - ".
>
>
>
> I hope you understand the problem and are not overwhelmed by this wall of text.
>
> I appreciate any suggestions.
>
>
>
> Cheers,
>
> Thomas
I am afraid that IDL will not be able to help you without some reformulation of your problem.
In order to avoid underflow and overflow when computing each of your Y and Z functions, you have to find a derived or approximated expression for their product, which indeed is finite and of order about 10.
You might for instance consider Rational Chebyshev approximations of X*Y, which are often used for computing the "erfcx" function (i.e. exp(x^2)*erfc(x)), whose shape is not far from the one you are dealing with.
Hoping this can help you.
alx.
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
