Implementation of the FFT for Statistical Functions

In this section we demonstrate how to use the FFT to calculate a density from a characteristic function.  Stable distributions are used in the example, because they are otherwise difficult to calculate in the probability domain.  The techniques, however, can be generalized to any characteristic function.

The functions must be set up in the correct order for everything to work.  For the FFT in the transform domain the list of values is rearranged.  If the characteristic function φ(t) is mapped across an ascending list corresponding to values of t, the list must be rearranged so that it begins at 0.  The left half of the list is placed at the right.  It is critical that the list begin with 0.
Normal ordering:

Map[φ[#] &, Range[-8, 7]]

ReOrder[list] reorders a frequency domain list for the time domain and vice versa.  If Length[list] is odd the output list will be truncated to the next smallest integer.

Rearranged ordering:

Map[φ[#] &, ReOrder[Range[-8, 7]]]

Using ReOrder twice restores the order.

ReOrder[%]

Since the statistical distribution is infinite, only part of it can be used with the FFT; it is important to have a large enough section of the distribution.  A segment which contains nearly all of the probability of the distribution this range is mapped across the N samples; it must also have symmetry on the x-axis.  SR is the sampling rate.  The width of the x-axis is approximately N/SR, which will be used as a scaling factor.   This interval must contain nearly all the probability density.  It is divided into N equally spaced intervals from {-N/(2 SR), N/(2 SR)-1/SR}. The functions are set up only for even numbered samples.  For use with the FFT, the function rearranges the order so that it begins at  t = 0, where  φ(0) = 1.  The amplitude of the whole sample must be rescaled for the probability domain so that it's integral will be 1.  When it is recast in the probability domain, the function XTable creates the proper x values and the order used for the InverseFourier is put back to normal using ReOrder.

DFTCF[cf,SR,N] puts out a discrete FT of a distribution specified by the characteric function cf in the frequency domain compatible with FourierParameters → {-1, 1}.  SR is the sampling rate; SR is greater than 0.  N is the number of points and must be even.  The output is scaled ready for inversion.

CF1c[u, α, β, γ, δ] outputs a compiled version of the S(α, β, γ, δ; 1) characteristic function.  All the inputs must be real numbers.

cf = DFTCF[CF1c[#1, 1.5, 0.5, 1, 0] &, 16, 8192];
ListPlot[Abs[cf], PlotRange -> All, PlotLabel -> "Abs[Characteristic Function]"]
ListPlot[Arg[cf], PlotRange -> All, PlotLabel -> "Arg[Characteristic Function]"]
pdf = Chop[InverseFourier[cf]];
g1 = ListPlot[Transpose[{XTable[16, 8192], ReOrder[pdf]}], PlotRange -> {{-25, 25}, All}, PlotLabel -> "Density Function"]

The plot of the inverse Fourier transform of the sampled characteristic function is the sampled density function.  Interpolation can be used to obtain values between the points.

Everything can also be done from the probability domain to create a discrete characteristic function.  The scaling is not needed in this case, because cf1[[1]] wasn't adjusted to have a value of 1.  This method of deriving the discrete characteristic function is less precise, but ultimately the numerical differences are small.  

cf1 = Quiet[Chop[Fourier[ReOrder[(SPDF[#1, {1.5, 0.5, 1, 0}] &) /@ XTable[16, 8192]]]]];
ListPlot[Abs[cf1], PlotRange -> All, PlotStyle -> Blue, PlotLabel -> "Abs[Characteristic Function]"]
ListPlot[Arg[#1] & /@ cf1, PlotRange -> All, PlotStyle -> Blue, PlotLabel -> "Arg[Characteristic Function]"]
pdf1 = Chop[InverseFourier[cf1]];
g2 = ListPlot[Transpose[{XTable[16, 8192], ReOrder[pdf1]}], PlotRange -> {{-25, 25}, All}, PlotStyle -> Blue, PlotLabel -> "Density Function"]

© Copyright 2008 mathestate    Fri 29 Feb 2008