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.
TransformedDistribution[u + v + w, {u \[Distributed] RayleighDistribution[σ], v \[Distributed] RayleighDistribution[σ], w \[Distributed] RayleighDistribution[σ]}]. $\endgroup$NIntegraterather thanIntegratefor 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.) $\endgroup$fLand $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. $\endgroup$