| Multiplying very high with very low numbers: erfc * exp [message #88263] |
Thu, 03 April 2014 02:35  |
tho.siegert
Messages: 6
Registered: April 2014
|
Junior Member |
|
|
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
|
|
|
|
| Re: Multiplying very high with very low numbers: erfc * exp [message #88264 is a reply to message #88263] |
Thu, 03 April 2014 06:21  |
lecacheux.alain
Messages: 325
Registered: January 2008
|
Senior Member |
|
|
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.
|
|
|
|
| Re: Multiplying very high with very low numbers: erfc * exp [message #88352 is a reply to message #88264] |
Wed, 16 April 2014 05:09 |
tho.siegert
Messages: 6
Registered: April 2014
|
Junior Member |
|
|
On Thursday, April 3, 2014 3:21:31 PM UTC+2, alx wrote:
> 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.
On Thursday, April 3, 2014 3:21:31 PM UTC+2, alx wrote:
> 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.
Okay, since I already actually had a lnerfc function, but was too silly to make it work properly, i post my solution:
function lnerfc, x, y
a = [-1.26551223d, 1.00002368d, 0.37409196d, 0.09678418d, -0.18628806d, 0.27886807d, -1.13520398d, 1.48851587d, -0.82215223d, 0.17087277d]
t = 1d / (1d + 0.5d * abs(x))
tau = t * exp( -x*x + (a[0] + t * (a[1] + t * (a[2] + t * (a[3] + t * (a[4] + t * (a[5] + t * (a[6] + t * (a[7] + t * (a[8] + t * a[9]))))))))))
y = alog(t) + ( -x*x + (a[0] + t * (a[1] + t * (a[2] + t * (a[3] + t * (a[4] + t * (a[5] + t * (a[6] + t * (a[7] + t * (a[8] + t * a[9]))))))))))
lt0 = where(x lt 0d,/null)
y[lt0] = y[lt0] + alog(2d / tau - 1d)
return, y
end
It is again taken from Numerical recipes, Chapter 6.2, Special Functions, just translated to logarithm space. This is indeed based on Chebyshev fitting. Thanks alx!
Regards,
Thomas
|
|
|
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
