Convolution of two probability distributions using FFT

In summary, the conversation is about coding an algorithm to compute the convolution of two probability distributions using the FFT. The user is trying to understand why their results do not correspond to any definition of convolution. Another user suggests that it may be a case of circular convolution, where the indices of the distributions are shifted. This explanation is consistent with the results and other tests, providing clarification for the original user.
  • #1
jamie_m
14
0
I've been trying to code an algorithm to compute the convolution of two probability distributions. using the FFT. This relies on the "convolution theorem":

(p*q)[z] = FFT^{-1}(FFT(p) \cdot FFT(q))

However, when I test it using the distributions

p={0.1, 0.2, 0.3, 0.4}
q={0.4, 0.3, 0.2, 0.1}

the result that's output - {0.24, 0.22, 0.24, 0.3} - doesn't correspond to any definition or generalised definition of convolution that I can find.

(I'm fairly sure that one of the 0.24s, for instance, results from p[0]q[2] + p[1]q[1] + p[2]q[0] + p[3]q[3] being computed and divided by n)

Can someone well-versed in FFTs (and/or convolutions) shed any light on what might be going on? I enclose the source code, which uses C++ and the FFTW package, below:

Code:
void convolution_of_two_probability_distributions_with_fft(long double* p, long double* q, long double* p_star_q, unsigned long int n)
{
	fftwl_plan fpcp;
	fftwl_plan fpcq;
	fftwl_plan f_reverse;
	unsigned long int ulx;

	fftwl_complex* transformed_p = (fftwl_complex*)fftwl_malloc(sizeof(fftwl_complex) * n);
	fftwl_complex* transformed_q = (fftwl_complex*)fftwl_malloc(sizeof(fftwl_complex) * n);
	fftwl_complex* transformed_p_dot_transformed_q = (fftwl_complex*)fftwl_malloc(sizeof(fftwl_complex) * n);

	fpcp = fftwl_plan_dft_r2c_1d(n, p, transformed_p, FFTW_ESTIMATE);
	fpcq = fftwl_plan_dft_r2c_1d(n, q, transformed_q, FFTW_ESTIMATE);
	f_reverse = fftwl_plan_dft_c2r_1d(n, transformed_p_dot_transformed_q, p_star_q, FFTW_ESTIMATE);

	fftwl_execute(fpcp);
	fftwl_execute(fpcq);

	for (unsigned long int ulx=0; ulx < n; ulx++)
	{
		//transformed_p_dot_transformed_q[ulx] = transformed_p[ulx] * transformed_q[ulx];
		//There's no overloaded * for fftwl_complex.
		transformed_p_dot_transformed_q[ulx][0] = (transformed_p[ulx][0] * transformed_q[ulx][0]) - (transformed_p[ulx][1] * transformed_q[ulx][1]);
		transformed_p_dot_transformed_q[ulx][1] = (transformed_p[ulx][0] * transformed_q[ulx][1]) + (transformed_p[ulx][1] * transformed_q[ulx][0]);
	}

	fftwl_execute(f_reverse);

	for (ulx=0; ulx < n; ulx++)
	{
		p_star_q[ulx] /= n;
	}

	fftwl_destroy_plan(fpcp);
	fftwl_destroy_plan(fpcq);
	fftwl_destroy_plan(f_reverse);

	fftwl_free(transformed_p);
	fftwl_free(transformed_q);
	fftwl_free(transformed_p_dot_transformed_q);
}
 
Last edited:
Mathematics news on Phys.org
  • #2
Perhaps it's "circular convolution"

[itex] (p*q)_n = \sum_{m=0}^3 { p_m q_{n-m} } [/itex]

e.g.
[itex] (p*q)_0 = \sum_{m = 0}^3 {p_m} q_{0-m} [/itex]
[itex] = p_0 q_0 + p_1 q_{-1} + p_2 q_{-2} + p_3 q_{-3} [/itex]

Where we regard [itex] q_{-1} = q_3 , q_{-2}= q_2 [/itex] etc

[itex] = (0.1)(0.4) + (0.2)(0.1) + (0.3)(0.2) + (0.4)(0.3) = 0.24 [/itex]
 
  • #3
Stephen Tashi said:
Perhaps it's "circular convolution"

[itex] (p*q)_n = \sum_{m=0}^3 { p_m q_{n-m} } [/itex]

I think you're right - this is consistent with the existing results and with another pair of distributions I've tested it with. Thanks for that!
 

1. What is the purpose of convolving two probability distributions using FFT?

The purpose of convolving two probability distributions using FFT (Fast Fourier Transform) is to find the probability distribution of the sum of two independent random variables. This is useful in many areas of science and engineering, such as signal processing, image processing, and probability theory.

2. How does FFT help in computing the convolution of two probability distributions?

FFT is an efficient algorithm for computing the discrete Fourier transform, which is used to represent a signal or function as a sum of complex sinusoidal functions. By taking the FFT of the two probability distributions, multiplying them together, and then taking the inverse FFT, we can compute the convolution of the two distributions much faster than using traditional methods.

3. What are the advantages of using FFT for convolution of probability distributions?

The main advantage of using FFT for convolution of probability distributions is speed. FFT can greatly reduce the computation time compared to traditional methods, especially when dealing with large datasets. Additionally, FFT is a stable and accurate algorithm, providing reliable results.

4. Can FFT be used for any type of probability distribution?

Yes, FFT can be used for any type of probability distribution as long as it is discrete and has a finite range. It is commonly used for Gaussian distributions, but it can also be applied to other distributions such as Poisson, binomial, and exponential.

5. Are there any limitations or drawbacks to using FFT for convolution of probability distributions?

One limitation of using FFT for convolution of probability distributions is that it assumes the distributions are independent. If the distributions are not independent, then FFT will not produce accurate results. Additionally, FFT is only suitable for discrete distributions, so it cannot be used for continuous distributions without discretization.

Similar threads

Replies
1
Views
1K
  • General Math
Replies
14
Views
4K
  • Set Theory, Logic, Probability, Statistics
Replies
0
Views
963
  • Set Theory, Logic, Probability, Statistics
Replies
3
Views
1K
Replies
2
Views
2K
  • Set Theory, Logic, Probability, Statistics
Replies
2
Views
2K
Replies
3
Views
2K
  • Electrical Engineering
Replies
1
Views
873
  • Differential Equations
Replies
2
Views
1K
  • Calculus and Beyond Homework Help
Replies
4
Views
984
Back
Top