# Using Convolve with functions that contain error functions, and iterative convolutions

I need to produce convolutions of convolutions to study some probability distributions. My starting distribution is the Rayleigh distribution, which I convolve with itself as:

RayleighRayleighConvolved = Convolve[PDF[RayleighDistribution[σ]][x], PDF[RayleighDistribution[σ]][x], x, X];
RRDistribution[σ_] := Evaluate[ProbabilityDistribution[RayleighRayleighConvolved,{X,-Infinity,Infinity}]];
PDF[RRDistribution[σ]]


Mathematically this would be $$R_{2} = (R * R)$$ where $$R$$ is the Rayleigh distribution. The result contains error functions (Erf[]).

This executes very quickly and behaves as expected:

Plot[PDF[RRDistribution[1]][X], {X, 0, 5}]


I now want to convolve by new distribution with the Rayleigh distribution again so this would be $$R_{3} = ((R * R) * R) = (R_{2} * R)$$

I try this in Mathematica with the same approach as above

Convolve[PDF[RRDistribution[σ]][x], PDF[RayleighDistribution[σ]][x], x, X]


But this fails just returning the command in input form.

I will briefly recap in integral notation what I am trying to achieve: $$\hat{R}_{N} = \int_{-\infty}^{+\infty}\hat{R}(t)_{N-1} R(t - \tau) \ d \tau$$ Where $$R$$ is just the Rayleigh distribution and $$\hat{R}_{N}$$ is the $$N$$'th convolution. So for $$N=1$$ the convolution would be the Rayleigh with itself: $$\hat{R}_{N=1} = \int_{-\infty}^{+\infty}\hat{R}(t)_{N=0} R(t - \tau) \ d \tau \\ \hat{R}_{N=1} = \int_{-\infty}^{+\infty}R(t) R(t - \tau) \ d \tau$$ I can't get any further than $$N = 1$$. I think this is because the result, $$\frac{1}{4 \sigma^{3}}\exp\left(\frac{-x^{2}}{2\sigma^{2}}\right) \left( 2 x \sigma + \exp\left(\frac{x^{2}}{4\sigma^{2}}\right) \sqrt{\pi} \left( x^{2} - 2 \sigma^{2} \right) \rm{erf}(x/2 \sigma) \right)$$, which comes from running

RayleighRayleighConvolved = Convolve[PDF[RayleighDistribution[σ]][x], PDF[RayleighDistribution[σ]][x], x, X];


contains an error function, where there is no closed solution; explaining why I can't go any further.

The papers recommended to me by J.M (thanks for these), show methods of approximating what I want to achieve.

It was also suggested to use compute in Mathematica with

TransformedDistribution[...]


If I do this for $$N = 1$$, that is

TransformedDistribution[u + v , {u \[Distributed] RayleighDistribution[\[Sigma]], v  \[Distributed] RayleighDistribution[\[Sigma]]}]


Which when plotted produces the same result I got using Convolve. When I exetend for

TransformedDistribution[u + v + w, {u \[Distributed] RayleighDistribution[\[Sigma]], v  \[Distributed] RayleighDistribution[\[Sigma]], w \[Distributed] RayleighDistribution[\[Sigma]]}]


Which takes an extremely long time to compute, and in fact I haven't seen it complete.

I want to use the results in MLE type evaluations like in FindDistributionParameters[] so speeding up execution would be extremely useful!

I looked at the papers suggested in the comments. The second one especially very interesting. In the paper they approximated the $$n$$'th convolved Rayleigh distribution as: $$f_{L}(t) = \frac{t^{2 L - 1} \exp\left( - \frac{t^{2}}{2b}\right) }{2^{L - 1}B^{L} (L - 1)!} - \frac{(t - a_{2})^{2L - 2} \exp\left( - \frac{a_{1}(t - a_{2})^{2}}{2b} \right)}{2^{L-1} b \left( \frac{b}{a_1}\right)^{L}(L - 1)!} a_{0} \left( b (2 L t - a_{2}) - a_{1}t(t - a_{2})^{2} \right)$$ where $$b = \frac{\sigma^{2}}{L}((2L - 1)!!)^{1/L}$$, the constants $$a_{i=0,1,2}$$ are dependant on the number of Rayleigh averages/convolutions, $$L$$. Lets put this into MM:

    b[\[Sigma]_, L_] := \[Sigma]^2/L ((2L-1)!!)^(1/L)
fL[b_, L_, a0_, a1_, a2_, t_] :=( t^(2L-1) Exp[-(t^2/(2b))])/(2^(L-1) b^L (L - 1)!) - ((t - a2)^(2L-2) Exp[-a1 (t-a2)^2/(2b)])/(2^(L-1) b (b/a1)^L (L - 1)!) a0 (b (2 L t - a2) - a1 t (t - a2)^2)

La0a1a2 = {{"L","a0","a1","a2"},{3,0.0164,0.306,0.9928},{4,0.0198,0.2413,0.976},{5,0.0221,0.1972,0.9654},{6,0.0236,0.1645,0.9583},{7,0.0248,0.1386,0.9531},{8,0.0257,0.1172,0.9491},{9,0.0264,0.0989,0.946},{10,0.027,0.0829,0.9434},{11,0.0275,0.0686,0.9412},{12,0.0279,0.0557,0.9393},{13,0.0283,0.044,0.9377},{14,0.0286,0.033,0.9363},{15,0.0288,0.0229,0.935},{16,0.0291,0.0133,0.9338}};


If we plot for $$L = 4$$

L = 4;
Show[
Histogram[Mean[Table[RandomVariate[RayleighDistribution[1], 100000], {i, 1, L}]], "FreedmanDiaconis", "PDF"],
Plot[
fL[b[0.52, La0a1a2[[L - 1]][[1]]],La0a1a2[[L - 1]][[1]], La0a1a2[[L - 1]][[2]], La0a1a2[[L - 1]][[3]], La0a1a2[[L - 1]][[4]], t +0.01],
{t, 0 , 8}, PlotRange->All
], PlotRange->{{0, 4}, All}
]


One can see that this is a pretty good approximation however after the main crest of the distribution we see this smaller bump.

• This is a well-studied problem; see e.g. this paper or this paper. I don't think there is a known closed form, but you can at least express the distribution as TransformedDistribution[u + v + w, {u \[Distributed] RayleighDistribution[σ], v \[Distributed] RayleighDistribution[σ], w \[Distributed] RayleighDistribution[σ]}]. – J. M.'s ennui Apr 26 '20 at 12:08
• Perhaps you should edit your question to include what you discovered, before anything else. ;) – J. M.'s ennui Apr 26 '20 at 13:01
• OK. Why not use the papers that @J.M. mentioned to find an approximation to $(R*R)*R$ and compare that of $(R*R)*R$ using NIntegrate rather than Integrate for a few values of $\sigma$? If it close enough for $N=2$, then I would bet money the approximations will be more than close enough for $N>2$. But this assumes that you know how you're characterizing "closeness" and what value of closeness you need. Or is it "I'll known it when I see it." ? (I'm intentionally playing Devil's advocate here.) – JimB Apr 28 '20 at 23:26
• With fL and $L=3$ can you duplicate the points in Figure 5 of that second article? I'm getting negative values for values of $t>4$ while the figure shows those as positive. Also, that article doesn't seem to mention how they got the exact densities. – JimB Apr 29 '20 at 21:54
• I think I figured it out. In that second article, Table 1 has the numbers in columns $a_1$ and $a_2$ reversed. That's why using the coefficients as is results in negative densities in the right tail. – JimB May 2 '20 at 5:07

This is only a partial answer. (I cannot duplicate the results from the approximation found in this article. That's likely my fault.)

The problem can be described as finding the pdf for the sum of $$L$$ independent and identically distributed Rayleigh random variables:

$$Z_L=\sum_{i=1}^L X_i$$

where $$X_i \sim \text{Rayleigh}(\sigma)$$. All of the moments can be determined for all values of $$L$$ but there doesn't appear to be a closed form for the pdf with $$L>3$$. The mean and standard deviation of $$Z_L$$ is given by

mean[L_, σ_] := L σ Sqrt[π/2];
sd[L_, σ_] := σ Sqrt[L (2 - π/2)];


The exact symbolic equations for the pdf when $$L=1$$ and $$L=2$$ are found with

pdf[1] = PDF[RayleighDistribution[σ], x][[1, 1, 1]];
pdf[2] = Integrate[pdf[1]*(pdf[1] /. x -> z - x), {x, 0, z},
Assumptions -> z > 0] /. z -> x;


(The inclusion of /. z -> x in pdf[2] will hopefully become apparent below.)

For $$L\geq 3$$ numerical integration is used and a "dense" table of pdf values are created. A reasonable value for the upper limit of integration is

zmax[L_, σ_] := mean[L, σ] + 5 sd[L, σ];


Below we set $$\sigma=1$$ to obtain a specific example:

σ = 1;
n = 100;  (* Number of values for each pdf *)
pdf[3] = Table[{z, NIntegrate[pdf[2]*(pdf[1] /. x -> z - x), {x, 0, z}]},
{z, 0, zmax[3, σ], zmax[3, σ]/(n - 1)}];
pdf[4] = Table[{z, NIntegrate[pdf[2]*(pdf[2] /. x -> z - x), {x, 0, z}]},
{z, 0, zmax[4, σ], zmax[4, σ]/(n - 1)}];
pdf[5] = Table[{z, NIntegrate[(pdf[2] /. x -> x1)*(pdf[2] /. x -> x2)*(pdf[1] /. x -> z - x1 - x2),
{x1, 0, z}, {x2, 0, z - x1}]}, {z, 0, zmax[5, σ], zmax[5, σ]/(n - 1)}];
pdf[6] = Table[{z, NIntegrate[(pdf[2] /. x -> x1)*(pdf[2] /. x -> x2)*(pdf[2] /.
x -> z - x1 - x2), {x1, 0, z}, {x2, 0, z - x1}]}, {z, 0, zmax[6, σ], zmax[6, σ]/(n - 1)}];
pdf[7] = Table[{z, NIntegrate[(pdf[2] /. x -> x1)*(pdf[2] /. x -> x2)*(pdf[2] /.
x -> x3)*(pdf[1] /. x -> z - x1 - x2 - x3), {x1, 0, z}, {x2, 0, z - x1}, {x3, 0, z - x1 - x2}]},
{z, 0, zmax[7, σ], zmax[7, σ]/(n - 1)}];
pdf[8] = Table[{z, NIntegrate[(pdf[2] /. x -> x1)*(pdf[2] /. x -> x2)*(pdf[2] /.
x -> x3)*(pdf[2] /. x -> z - x1 - x2 - x3), {x1, 0, z}, {x2, 0, z - x1},
{x3, 0, z - x1 - x2}]}, {z, 0, zmax[8, σ], zmax[8, σ]/(n - 1)}];


A plot of all 8 pdf's:

Show[Plot[{pdf[1], pdf[2]}, {x, 0, zmax[8, σ]}, PlotRange -> All, PlotStyle -> Black],
ListPlot[pdf[#] & /@ Range[3, 8], Joined -> True, PlotStyle -> Black]]


The article mentioned above provides an excellent approximation for $$3\leq L \leq 16$$ but it seems to have a minor typo but with serious consequences. Fortunately it is easily fixed. In Table 1 the values for columns $$a_1$$ and $$a_2$$ are switched. Below is shown the justification for that statement followed by a Mathematica function to implement the excellent approximation in the article.

In short, the article shows a simple approximation to the cdf and enhances that simple approximation by fitting a function to the deviations from the simple approximation. The estimated deviation function is then subtracted from the simple approximation and the derivative is taken to obtain an approximation to the pdf.

(* Set parameters *)
σ = 1;
L = 3;

(* Some useful functions from above *)
mean[L_, σ_] := L σ Sqrt[π/2];
sd[L_, σ_] := σ Sqrt[L (2 - π/2)];
zmax[L_, σ_] := mean[L, σ] + 5 sd[L, σ];

(* Determine "true" cdf using numerical integration *)
n = 100;  (* Number of points to evaluate cdf *)
pdf[1] = PDF[RayleighDistribution[σ], x][[1, 1, 1]];
pdf[2] = Integrate[pdf[1]*(pdf[1] /. x -> z - x), {x, 0, z},
Assumptions -> z > 0] /. z -> x;
pdfL = Table[{z, NIntegrate[pdf[2]*(pdf[1] /. x -> z - x), {x, 0, z}]},
{z, zmax[L, σ]/(n - 1), zmax[L, σ],
zmax[L, σ]/(n - 1.)}];
cdfL = Table[
NIntegrate[pdf[2]*(pdf[1] /. x -> z - x), {z, 0, zz}, {x, 0, z}],
{zz, zmax[L, σ]/(n - 1), zmax[L, σ],
zmax[L, σ]/(n - 1.)}];

(* Determine SAA approxmation to cdf of t=Z/Sqrt[L] *)
b = σ^2 ((2 L - 1)!!)^(1/L) /L;
FSAA = Table[1 - Exp[-(zz/Sqrt[L])^2/(2 b)] Sum[((zz/Sqrt[L])^2/(2 b))^k / k!,
{k, 0, L - 1}], {zz, zmax[L, σ]/(n - 1), zmax[L, σ], zmax[L, σ]/(n - 1.)}];


Now find the difference between the approximation and the true value and perform the regression:

err = FSAA - cdfL;

(* Create a dataset and estimate coefficients of function that describes the error *)
tt = Table[zz/Sqrt[L], {zz, zmax[L, σ]/(n - 1), zmax[L, σ], zmax[L, σ]/(n - 1.)}];
data = Transpose[{tt, err}];
nlm = NonlinearModelFit[data, t a0 (t - a2)^(2 L - 1) Exp[-a1 (t - a2)^2/(2 b)]/(2^(L - 1) (b/a1)^L (L - 1)!),
{{a0, 0.0164}, {a1, 0.3060}, {a2, 0.9928}}, t];
nlm["BestFitParameters"]
(* {a0 -> 0.0163999, a1 -> 0.992766, a2 -> 0.305966} *)


Ones sees that the values of a1 and a2 are reversed in Table 1 of the article.

Finally below is a function that will approximate the pdf of a sum of $$L$$ Rayleigh random variables:

(* Define constants *)
aa0 = {0, 0, 0.0164, 0.0198, 0.0221, 0.0236, 0.0248, 0.0257, 0.0264,
0.027, 0.0275, 0.0279, 0.0283, 0.0286, 0.0288, 0.0291};
aa1 = {0, 0, 0.9928, 0.976, 0.9654, 0.9583, 0.9531, 0.9491, 0.946,
0.9434, 0.9412, 0.9393, 0.9377, 0.9363, 0.935, 0.9338};
aa2 = {0, 0, 0.306, 0.2413, 0.1972, 0.1645, 0.1386, 0.1172, 0.0989, 0.0829,
0.0686, 0.0557, 0.044, 0.033, 0.0229, 0.0133};

(* Set parameters *)
σ = 1;
L = 3;

(* Some useful functions from above *)
mean[L_, σ_] := L σ Sqrt[π/2];
sd[L_, σ_] := σ Sqrt[L (2 - π/2)];
zmax[L_, σ_] := mean[L, σ] + 5 sd[L, σ];

pdfApprox[z_, L_, σ_, a0_, a1_, a2_] := Module[{b, t},
b = σ^2 ((2 L - 1)!!)^(1/L) /L;
t = z/Sqrt[L];
t^(2 L - 1) Exp[-t^2/(2 b)]/(2^(L - 1) b^L (L - 1)!)/Sqrt[L] -
(t - a2)^(2 L - 2) Exp[-a1 (t - a2)^2/(2 b)] a0*
(b (2 L t - a2) - a1 t (t - a2)^2)/(2^(L - 1) b (b/a1)^L (L - 1)! Sqrt[L])]

Show[ListPlot[pdfL, Joined -> True,
PlotStyle -> {{Red, Thickness[0.01]}}],
Plot[pdfApprox[z, L, σ, aa0[[L]], aa1[[L]], aa2[[L]]], {z, 0,
zmax[L, σ]}, PlotStyle -> Blue]]


No real need for a legend because the approximation is so good. The pdf's are right on top of each other.

• JimB this is fantastic, thanks! I was previously using the $f_{SAA}$ as, until your spot of the swapped $a_{1}$ and $a_{2}$ columns, this gave the best estimate. I now have three models to play with!! Seeing as you provided an answer and a solution I think you win this one (again). – Q.P. May 2 '20 at 9:31
• With getting negative estimates of density in the right tail I checked all of the equations starting with differentiating the cdf to obtain the pdf. That all checked out. Then for $L=3$ I created a table of points for the cdf and used NonlinearModelFit to check on the constants in Table 1. I got the exact same values in the table but only if $a1$ and $a2$ were switched. I'll add those steps sometime this weekend. The Journal that published the article doesn't seem to have errata for that article. – JimB May 2 '20 at 13:58
• I have sent an e-mail to the authors of the article about this. – JimB May 3 '20 at 0:40
• Thanks again Jim! You've really gone above and beyond with this question! – Q.P. May 3 '20 at 9:33
• I found a way to get a straightforward and accurate approximation of the pdf for values of $L$ up to 107. (And I'm sure that the numerical instability issues that start with $L=108$ can be fixed with better programming.) Is that of interest? – JimB May 5 '20 at 2:15

Disclaimer: I am not up-to-speed on any current literature on this particular approximation. What I'm about to present has likely been done many years ago.

The problem can be described as the need to find the probability density function (or a good approximation of it) for the sum of $$L$$ independent and identically distributed Rayleigh random variables with parameter $$\lambda$$.

$$Z_L=\sum_{i=1}^L X_i$$

with $$X_i\sim \text{Rayleigh}(\sigma)$$.

The following article described a simple way to approximate the pdf of $$Z$$ for $$L=1, 2,\ldots,16$$ (given that the true pdf was known). (Now that isn't a logical contradiction. The actual distribution for $$L>2$$ can almost certainly be computed using numerical integration and only of that is programmed very carefully. What is presented is an approximation that doesn't need numerical integration.) The article is

Hu, J., Beaulieu, N. C. (2005) Accurate simple closed-form approximations to Rayleigh sum distributions and densities. IEEE Communications Letters 9: 109–111.

The authors discovered that the difference between the cumulative distribution functions (cdf's) of the true distribution and a crude approximation followed a common functional form with 3 parameters that could be estimated using regression for each value of $$L$$. That would provide the adjustment to get a good estimate of the cdf. Then that cdf could be differentiated to obtain the pdf.

But for large values of $$L$$ we don't have the true cdf. Fortunately the next best thing to having the true cdf is being able to generate a gazillion samples from that distribution. (This forms the basis for current widespread use of Bayesian statistics.)

So we generate a large number of samples from the distribution of $$Z_L$$, find the empirical cdf and then proceed as the article did. This allows us to get a good approximation to the true pdf. (However, for values much above 100 the following code has numerical instability problems. Those are relatively easy to fix but someone else can do that if there is interest. And, again, it's very likely that someone has already published about this.)

The code follows what was done in the article and much of the code deals with $$t=Z_/\sqrt{L}$$ rather than $$Z_L$$ directly. The function pdfAdjusted does provide the adjusted pdf of $$Z_L$$.

(* Define some functions to determine what values of the random variable to be concerned about *)
mean[L_, σ_] := L σ Sqrt[π/2];
sd[L_, σ_] := σ Sqrt[L (2 - π/2)];
zmax[L_, σ_] := mean[L, σ] + 5 sd[L, σ];
zmin[L_, σ_] := Max[0.1, mean[L, σ] - 5 sd[L, σ]]

(* Set parameters *)
L = 64;
σ = 1;

(* Generate a large number of random samples of t = X/Sqrt[L] and
construct the empirical cdf function at a fixed set of n values *)
sampleSize = 100000;
data = Total[#]/Sqrt[L] & /@
RandomVariate[RayleighDistribution[σ], {sampleSize, L}];
dist = EmpiricalDistribution[data];
n = 100;
tt = Range[zmin[L, σ], zmax[L, σ], (zmax[L, σ] - zmin[L, σ])/(n - 1.)]/Sqrt[L];
empiricalCDF = CDF[dist, #] & /@ tt;

(* Function that crudely approximates the cdf of t *)
cdfSAA[t_, L_, σ_] := Module[{b},
b = σ^2 ((2 L - 1)!!)^(1/L) /L;
1 - Exp[-t^2/(2 b)] Sum[(t^2/(2 b))^k/k!, {k, 0, L - 1}]]

(* Generate cdfSAA values for each value in tt *)
approxCDF = Max[0, cdfSAA[#, L, σ]] & /@ tt // N;

(* Construct a dataset with t values, the empirical CDF, the approximate CDF
and the difference between the two. *)
(* Only keep those rows with 10^-7 < empiricalCDF < 1 - 10^-7 *)
d = Transpose[{tt, empiricalCDF, approxCDF, approxCDF - empiricalCDF}];
d = Select[d, 10^-7 < #[[2]] < 1 - 10^-7 &];

ListPlot[{d[[All, {1, 2}]], d[[All, {1, 3}]]},
PlotLegends -> {"Empirical CDF", "SAA CDF"},
Frame -> True, FrameLabel -> {"t", "CDF"},
PlotLabel -> "L = " <> ToString[L]]


(* Predict the difference between the approximate cdf and the \
empirical cdf *)
b = σ^2 ((2 L - 1)!!)^(1/L) /L;
nlm = NonlinearModelFit[d[[All, {1, 4}]], {t a0 (t - a2)^(2 L - 1)*
Exp[-a1 (t - a2)^2/(2 b)]/(2^(L - 1) (b/a1)^L (L - 1)!),
a0 > 0 && a1 > 0 && 0 <= a2 < Min[d[[All, 1]]]},
{{a0, 0.04}, {a1, 0.9}, {a2, 0}}, t];

(* Display the fit of the difference in the cdf's *)
Show[ListPlot[d[[All, {1, 4}]], PlotRangeClipping -> False,
PlotRange -> {{zmin[L, σ]/Sqrt[L], zmax[L, σ]/Sqrt[L]}, {0, Automatic}},
PlotLabel -> "L = " <> ToString[L], Frame -> True, FrameLabel -> {"t", "CDF error"}],
Plot[nlm[t], {t, zmin[L, σ]/Sqrt[L], zmax[L, σ]/Sqrt[L]}]]


(* Function that approximates the pdf of the sum of L Rayleigh random variables *)
pdfAdjusted[z_, L_, σ_, a0_, a1_, a2_] := Module[{b, t},
b = σ^2 Exp[Log[(2 L - 1)!!]/L]/L;
t = z/Sqrt[L];
t^(2 L - 1) Exp[-t^2/(2 b)]/(2^(L - 1) b^L (L - 1)!)/Sqrt[L] -
(t - a2)^(2 L - 2) Exp[-a1 (t - a2)^2/(2 b)] a0 (b (2 L t - a2) - a1 t (t - a2)^2)/
(2^(L - 1) b (b/a1)^L (L - 1)! Sqrt[L])]

(* Crude approximation of pdf *)
pdfSAA[z_, L_, σ_] :=
Module[{b, t}, b = σ^2 ((2 L - 1)!!)^(1/L) /L;
t = z/Sqrt[L];
t^(2 L - 1) Exp[-t^2/(2 b)]/(2^(L - 1) b^L (L - 1)!)/Sqrt[L]]

(* Show results *)
{aa0, aa1, aa2} = {a0, a1, a2} /. nlm["BestFitParameters"]
Show[Histogram[Sqrt[L] data, "FreedmanDiaconis", "PDF", Frame -> True,
FrameLabel -> {"Z", "PDF"}, PlotLabel -> "L = " <> ToString[L]],
Plot[{pdfSAA[z, L, σ], pdfAdjusted[z, L, σ, aa0, aa1, aa2]},
{z, zmin[L, σ], zmax[L, σ]},