# Distribution over the product of three, or n, independent Beta random variables

I would like to calculate the PDF for the product of three independent Beta random variables. Specifically, I would like to find the distribution of the product of the following: $X_1\sim Beta(1,3/2)$, $X_2\sim Beta(3/2,1)$ and $X_3\sim Beta(2,1/2)$.

I can manage to get Mathematica to output a distribution for the case of the the product $X_1 X_2$ by doing the following:

PDF[TransformedDistribution[
u v w, {u \[Distributed] BetaDistribution[1, 3/2],


This takes a few minutes to evaluate. However, when I move to $X_1 X_2 X_3$ Mathematica seems to be eternally stuck on the calculation:

PDF[TransformedDistribution[
u v w, {u \[Distributed] BetaDistribution[1, 3/2],


I could approximate this by a large number of draws from the distributions of $X_1$, $X_2$ and $X_3$; multiplying each set together. However, I would like the analytic form.

Does anyone have any idea how I can do this? The reason I state $n$ in the question is because I would like eventually to generalise this calculation to the product of more $Beta$ random variables.

Best,

Ben

• Might have better luck getting an answer on current SOTA for this over at Math.StackExchange.Com, then writing appropriate MMA code. I don't expect you'll find a compact/low computational complexity solution.
– ciao
Aug 11, 2015 at 22:09
• It has been a while since I've done something like this but if I recall correctly you may be able to do this quicker using the moment generation function of the Beta distributions, performing the product on these, and then matching that form to the form of a know distribution's moment generating function. Aug 11, 2015 at 22:16
• Here's a useful reference - extends earlier work that required integer parameters to (relatively) arbitrary parameters...
– ciao
Aug 11, 2015 at 22:23
• @Edmund good luck with that. Aug 12, 2015 at 8:06

Here is an answer that does not use a 3rd party package and works for an arbitrary amount of Beta distributions.

You can make use of a closed form for the product of n Beta distributions from the Handbook of Beta Distribution and Its Applications, Products and Linear Combinations, I. Products, B. Exact Distributions as found on page 57.

This expresses a product of n Beta distributions as the product of a constant $K$ (a function of the Beta distributions parameters in $\Gamma$ functions) and a Meijer-G function of the Beta distribution parameters.

ClearAll[pdfProductBeta];
pdfProductBeta[
α_ /; VectorQ[α, NumericQ],
β_ /; VectorQ[β, NumericQ],
x_Symbol] /; Dimensions[α] == Dimensions[β] :=
Module[{k},
k = Times @@ (Gamma[Plus @@ #]/Gamma[#[[1]]] & /@ (Transpose@{α, β}));
k MeijerG[{{}, α + β - 1}, {α - 1, {}}, x]
]


For the three Beta distributions in the question

a = {1, 3/2, 2}; b = {3/2, 1, 1/2};
pdfProductBeta[a, b, x]
(* 27/32 π MeijerG[{{}, {3/2, 3/2, 3/2}}, {{0, 1/2, 1}, {}}, x] *)


This result matches to @wolfies answer above that makes use of the mathStatica 3rd party package.

The two Beta case plots quickly so we can compare this easily to the TransfromedDistribution PDF from the built in Mathematica functions.

tpdf = PDF[
TransformedDistribution[
u v, {u \[Distributed] BetaDistribution[1, 3/2],

a = {1, 3/2}; b = {3/2, 1};
mpdf = pdfProductBeta[a, b, x];

GraphicsRow[
Plot[#, {x, 0, 1}, PlotRangePadding -> None] & /@ {tpdf, mpdf}]


For the three Beta case we need to limit the upper bound in the plot as it takes very long to calculate there.

a = {1, 3/2, 2}; b = {3/2, 1, 1/2};
mpdf3 = pdfProductBeta[a, b, x];
(* evaluate pdf at a point *)
mpdf3 /. x -> 0.5
(* 0.479319 *)
Plot[mpdf3, {x, 0, 0.9}, Exclusions -> {0}, AxesOrigin -> {0, 0},
PlotRangePadding -> None, PlotRange -> Full]


With the pdfProductBeta function you can construct the pdf for the product of an arbitrary number of Beta distributions without the need of a third party package.

Hope this helps.

• Very nice indeed to have a general solution from the Handbook of Beta Distribution -- I would go further and say that your answer not only constructs the pdf of the product of Beta random variables without the need of an add-on package, ... it does so without needing Mathematica either. Aug 12, 2015 at 14:03
• @wolfies Well ... I wouldn't want to try constructing, evaluating, or <cringe> plotting it without Mma. :-) It quickly gets process heavy for small $n$. Aug 12, 2015 at 14:09
• Note that this is no longer needed in version 11 and above as TransformedDistribution directly uses the above method for a product of BetaDistributions. Apr 8, 2017 at 13:18

Let random variable $X_i \sim Beta(a_i,b_i)$, with pdf $f_i(x_i)$. The OP is interested in 3 specific parameter combinations:

The pdf of $Y = X_2 X_3$, say $g(y)$, is:

where I am using the TransformProduct function from the mathStatica package for Mathematica, and where domain[g] = {y,0,1}.

The pdf of $Z = X_1 X_2 X_3 = Y* X_1$, say $h(z)$, is then:

All done.

Quick Monte Carlo check

It is always a good idea to check symbolic solutions with Monte Carlo methods.

The following plot compares:

• the exact symbolic solution pdf obtained $h(z)$ [ red dashed curve ], to

• an empirical Monte Carlo simulation of the pdf [ blue squiggly curve ]

All looks good.

• That Meijer's function had to show up indicates that things are not too simple. Tho, those parameters look familiar… what happens if you apply FunctionExpand[]? Aug 12, 2015 at 9:15
• Applying FunctionExpand returns the same MeijerG function. I might add that it appears tractable for most of the domain ... but trying to plot it for $x$ close to 1 becomes very slow (which is also why that is 'left out' in the diagram above). Aug 12, 2015 at 11:44
• Ah, that part I can explain; the contour integration being done under the hood converges rather slowly when the argument is near $1$. A pity that there does not seem to be a simpler form (and I was so sure that there was)… Aug 12, 2015 at 11:46
• @J. M. What is the reason the contour integration converges so slowly near 1? I am trying to generate points of the pdf near there, and was wondering if you had any ideas? Aug 12, 2015 at 14:19
• @ben, it's an inherent limitation of the contour integral method for the $G$-function. (The situation is similar to the one for generalized hypergeometric functions; arguments near $1$ are often very difficult cases, numerically speaking.) That's why I was hoping there might be expressions in terms of simpler special functions. I do not have Mathematica at hand to investigate further, unfortunately. Aug 12, 2015 at 15:04

The problem can indeed be solved explicitly for the product of n = 3 Beta-distributed variables and the explicit parameters of the OP.

In part 1 I show only the results, and turn later, in part 2, to the details of calculation in Mathematica, part 3 is discussion.

Part 1 Results

The PDF of the Beta distribution is given by

f[x_, a_, b_] = Simplify[PDF[BetaDistribution[a, b], x], 0 < x < 1]

(*
((1 - x)^(-1 + b) x^(-1 + a))/Beta[a, b]
*)


Let the random variables and their rescpective distributions be X1 ~ Beta(1,3/2) , X2 ~ Beta(3/2,1) and X3 ~ Beta(2,1/2), and let

f2n(x) = dist( X1*X2 )
f3n(x) = dist( X1*X2*X3 )

be the requested distributions.

Then we find

f2n[x_] = 9/2 (Sqrt[1 - x] - Sqrt[x] ArcTan[Sqrt[-1 + 1/x]]);

f3n[x_] = 27/4 (EllipticE[1 - x] - x EllipticK[1 - x]) -
27/16 \[Pi] MeijerG[{{}, {3/2, 3/2, 3/2}}, {{1/2, 1, 1}, {}}, x];


The functions are normalized

Integrate[f2n[x], {x, 0, 1}]

(* 1 *)

Integrate[f3n[x], {x, 0, 1}]

(* 1 *)


A plot of the two functions is shown here

Plot[{f2n[x], f3n[x]}, {x, 0.0001, 0.9}, PlotRange -> {{0, 1}, {0, 4}},
PlotLabel ->
"Distributions of the product of\ntwo (yellow) and three (blue)\nbeta \
distributed random variables", AxesLabel -> {"x", "f(x)"}]
(* 150812_Plot_Prod_Beta_dist.jpg *)


Since with the Meijer function Mathematica requires very long calculation times close to x = 1 I have left this region out.

Observations

1) I believe that the case of general n should be tackled using the Mellin transformation, as is natural for products (as is Fourier for sums). Our result for n = 3, the MeijerG function, already exhibits this pattern.

Part 2: Derivation

The distributions were calculated here using the general formula

  fn(x) = Integrate( Prod( du f(u,p)) DiracDelta(x-Prod(u)) )


Details will be given later.

Part 3: Discussion

Let me rather start a brief discussion suggested by a comment of Guess who it is.

2.1) The result of wolfie (which I saw only after having finished my calculations) is

f3wolfie[x_] :=
27/32 \[Pi] MeijerG[{{}, {3/2, 3/2, 3/2}}, {{0, 1/2, 1}, {}}, x]


which is not identical at first sight with my result

f3wolfgang[x_] :=
27/4 (EllipticE[1 - x] - x EllipticK[1 - x]) -
27/16 \[Pi] MeijerG[{{}, {3/2, 3/2, 3/2}}, {{1/2, 1, 1}, {}}, x]


But the functions indeed coincide as can be seen in a plot, or some clever FullSimplify?

For f3n[x], the last step in my calculation was this integral

Integrate[27/8 (
Sqrt[-((w x)/(-1 + w))] (Sqrt[-1 + w/x] - ArcCos[Sqrt[x/w]]))/w, {w, x, 1},
Assumptions - 0 < x < 1]


and the two terms in f3wolfgang correspond to the two terms in the this integrand.

2.2) Behaviour close to x = 1

As I haven't found a quick answer in the literature I tackled the poblem "experimentally":

Numerically my solution can be respresented as

f3nn[x_]:=NIntegrate[
27/8 ( Sqrt[-((
w x)/(-1 + w))] (Sqrt[-1 + w/x] - ArcCos[Sqrt[x/w]]))/w, {w, x,
1}]


Plotting this close to 1 with a negative power of (1-x) attached

Plot[1/(1 - x)^k f3nn[x], {x, 0.5, 1.1}]


and playing with k close to 2 shows that f3nn[x] ~= const (1-x)^2. I suspect even that k = 2 is exact because only for this value the function exhibits a sharp shoulder.

Afterwards I found (after some lengthy study) the exact behaviour of the constant to be 27 Pi/16, i.e. the distribution function has the series Expansion

f3n[x] = 27/64 \[Pi] (1 - x)^2 + O[1 - x]^3


Also, the MeijerG-function has the series expansion

MeijerG[{{}, {3/2, 3/2, 3/2}}, {{1/2, 1, 1}, {}}, x] == (1 - x) -
1/8 (1 - x)^2 + O[1 - x]^3


Remark: the developemnts here show that there is a simple integral representation of the MeijerG function, much simpler than the complex integral version. We have

MeijerG[{{}, {3/2, 3/2, 3/2}}, {{1/2, 1, 1}, {}}, x] = 2/Pi Integrate[
27/8 ( Sqrt[-((w x)/(-1 + w))] (Sqrt[-1 + w/x]
- ArcCos[Sqrt[x/w]]))/w, {w, x, 1}]


The difficulties of Mathematica calculating the numerical values close to 1 have been overcome by this representation and the series expansion.

• "using the Mellin transformation" - with the Meijer result, this is what you're implicitly doing anyway, since the $G$-function is effectively an inverse Mellin transform… Aug 12, 2015 at 12:29
• I have been busy doing the calculations and the editing so I did't see the resuls of wolfie. Sorry for that. Aug 12, 2015 at 12:30
• @J. M.: that's what I was - cautiously - saying ;-) Aug 12, 2015 at 12:34
• The reason why I had mentioned that the Meijer result looked familiar in the other answer is that it resembled some of the Mellin-Barnes representations for elliptic integrals that I remember; your result, which at least involves the complete elliptic integrals, seems to be tantalizing. Aug 12, 2015 at 12:41
• Some nice reference to MeijerG is ams.org/notices/201307/rnoti-p866.pdf pointing out specifically the closure under convolutions. Aug 12, 2015 at 13:20