Re: [Help-gsl] Does anybody know how to use FFT to compute numerical integration?

[Michael]

Hi Evgeny,

Here is my numerical difficulty I've met painstakingly if I don't
treat the Fourier integral carefully -- it is very tricky and
extremely challenging:

------------------

Is there a way to decompose the characteristic function into discrete part
and the continous part of the distribution and do the inverse transforms
only on the continous part?

I got very strange behavior on my inversions.

Let's fix the following notations:

Probability Density Function: f(x)
Cumulative Distribution Function: P(x)
Complementary Distribution Function: Q(x)=1-P(x)
Characteristic Function: F(v), which is the Fourier Transform of f(x),
defined as F(v)=integrate(f(x)*exp(i * v* x), x from -infinity to
+infinity).

I use the numerical intergral in Matlab:

f_hat(x)=integrate(F(v) * exp(-i * v* x), v from -infinity to +infinity),

And I used truncation: for large L and Matlab function "quadl"

f_hat(x)=integrate(F(v) * exp(-i * v* x), v from -L to +L),

and also I tried transformation:

f_hat(x)=integrate(F(atanh(y)) * exp(-i * atanh(y)* x)/(1-y^2),  y
from -1+epsilon to +1-epsilon),

---------------------------------

In all cases above, I got some of the recovered density function f_hat(x) to
be negative, which is wrong.

And when I sum the f_hat(x) up, and did sum(f_hat(x))*delta_x,

and I got this sum to be more than 1, which is impossible theoretically,
because I am actually approximating the cumulative probability distribution
function P(x), it should be between 0 and 1.

---------------------------------

I just don't understand how can a Fourier inverse transform of a
characteristic function, which is the probability density function, be
negative; and how can the cumulative function be more than 1?

It is worth noting that if I do the inverse Laplace transform to obtain the
cumulative distribution function directly from inverting F(v)/v, I got the
results agreeing with the above results. The cumulative function P(x) is
still non-monotonous and some places can go beyond 1.

It is really weird...

Anybody can help?




-------------------------------------------------------------------------------------------


I think the discrete component in the distribution is the bad apple.

I want to decompose the characteristic functiuon into discrete part and
continous part.

My characteristic function is of the form:

F(v)=exp{-a*i*v/(i*v-b)+ c*log(1-d*i*v)/(i*v-b)}

It looks like when v -> plus/minus infinity,

F(v) -> exp(-a).

I then claim that there is a "delta(x)*exp(-a)" component in the probability
density function f(x), which is a spike at zero.

I then tried to define:

F1(v)=exp{-a*i*v/(i*v-b)+ c*log(1-d*i*v)/(i*v-b)} - exp(-a),

and tried to obtain the continous part of the probability density function
f_c(x) from inverse Fourier transform of the above F1(v),

I thought now I can safely truncate the infinite integral at large
truncating point and do the numerical integral.

But it didn't help at all. Still the same wierd behavior occured.

And it is even worse. The obtained density function f_c(x) is smaller than
f_hat(x), and the cumulative function is also always smaller than the
previous cumulative function, and the cumulative function does not stabalize
to 1 at all. (Please be advised that my previous cumulative function
stabalize to 1 when x grows large, which is correct).

What's wrong?


On 6/28/07, Michael <address@hidden> wrote:
> Thanks! I probably have to do a smarter FFT... otherwise I don't think
> a plain discretization of the function and then take FFT I will be
> able to obtain a good Fourier integral. Fourier integral are hard
> integrals and I have met with enough numerical instability and
> difficulties... I am really looking for smarter FFT to do the
> numerical integration.
>
> Thanks a lot!
>
> On 6/28/07, Evgeny Kurbatov <address@hidden> wrote:
> > The sufficient example of interpolation is on
> > http://www.gnu.org/software/gsl/manual/html_node/Interpolation.html
> >
> > But you can realize your own minimal and efficient linear interpolation
> > on the grid with given mappings v(i) and i(v).  If the mapping i(v) is
> > known, you will not need for lookup of index 'i' by variable 'v', so it
> > will be fast...
> >
> > Michael wrote:
> > > It doesn't seem to be efficient to me? Will the interpolation be fast?
> > > How do you handle it? Thanks!
> >