# Limits for Triple Integration

I have 4 random variables: $p_1, p_2, p_3, p_4$

The joint probability distribution function of $p_1, p_2, p_3$ is:

$f(p1,p2,p3) = p_1^{b_1 + x_1 - 1} p_2^{b_2 + x_2 - 1} p_3^{b_3 + x_3 - 1} (1-p_1-p_2-p_3)^{x4+b4-1}$

where $x_1,x_2,x_3,x_4,b_1,b_2,b_3,b_4$ are positive real numbers.

I want to triple integrate f($p_1, p_2, p_3$) over the following criteria:

\begin{align*} p_1 + p_2 + p_3 + p_4 &= 1\\ p_1 + p_2 + p_3 &\le 1\\ p_1 &> p_2\\ p_1 &> p_3\\ p_1 &\ge 0\\ p_2 &\ge 0\\ p_3 &\ge 0\\ p_4 &\ge 0 \end{align*}

Can someone tell me what should be the limits of my triple integral and in which order f($p_1, p_2, p_3$) should be integrated?

• Generally, you can do integrals like these in two ways: you can define a region reg = ImplicitRegion[...] (or other valid region object) that satisfies the constraints and then call Integrate using this region as Integrate[f[p1,p2,p3,p4], {p1, p2, p3, p4} \[Element] reg]. The other way is to include the constraints into the integrant by multiplying f with statements like Boole[p1 > p3] to make it zero in the places where you don't want to have a contribution. For the constraint p1+p2+p3+p4 == 1, you probably need a DiracDelta Jun 19 '17 at 9:08
• Altho mef has already informed you about DirichletDistribution[], you might nevertheless be interested in looking up CylindricalDecomposition[] for determining your integration bounds. Aug 14 '17 at 14:44

Therefore, writing:

Reduce[p1 + p2 + p3 < 1 && p1 > p2 && p1 > p3 && p1 > 0 && p2 > 0 && p3 > 0]

you can also get the integration intervals.

Unfortunately, a simple analytical result doesn't seem to me to be able to get it.

Then, for example, writing:

f[b1_, b2_, b3_, b4_, p1_, p2_, p3_, p4_, x1_, x2_, x3_, x4_] :=
p1^(b1 + x1 - 1) p2^(b2 + x2 - 1) p3^(b3 + x3 - 1) p4^(b4 + x4 - 1);

A := ImplicitRegion[p1 + p2 + p3 < 1 && p1 > p2 && p1 > p3 &&
p1 > 0 && p2 > 0 && p3 > 0, {p1, p2, p3}];

Integrate[f[1, 1, 1, 1, p1, p2, p3, 1 - p1 - p2 - p3,
1, 1, 1, 1], {p1, p2, p3} \[Element] A]

I get:

1/15120

which is what you want.

• You missed the condition $p_1 + p_2 + p_3 + p_4 = 1$. And how can you get the integration result without the values of x1, x2, b1, b2, etc.? Jun 19 '17 at 10:42
• I did not forget it, that condition simply implies that p4=1-p1-p2-p3 and is taken into account in the integer. Unfortunately, that integrator does not allow a primitive expressible in terms of elementary functions, so you have to replace numeric values. If you explicitly explain the context, you can probably help it more effectively.
– TeM
Jun 19 '17 at 10:46
• $p_i$s are vote shares for party i. The joint probability distribution function is available. I want to calculate the probability when party 1 beats party 2 and party 3. The criteria there is $p_1 > p_2$ and $p_1 > p_3$. Also, sum of vote shares for these 3 parties will be less than equal to 1, i.e. $p_1 + p_2 + p_3 \leq 1$. And the sum of all vote shares should be 1, i.e. $p_1 + p_2 + p_3 + p_4 =1$ Jun 19 '17 at 11:01
• Unfortunately, I do not see an alternative to numerical calculation, i.e. the replacement of numeric values to the parameters. Naturally, this can be automated, for example through the Table[] function, but the analytical path does not seem to me to be feasible, MMA raises the white flag.
– TeM
Jun 19 '17 at 12:02

It looks like you are using the Dirichlet distribution, but you have omitted the constant of integration. For simplicity, let $a_i = b_i + x_i$. The density function for this distribution is given by

PDF[DirichletDistribution[{a1, a2, a3, a4}], {p1, p2, p3}]

In order to find the limits of integration that respect your additional constraints, it is convenient to work with a specialized parameter vector where $a_i = 1$. In this case the pdf is flat over the simplex. Before finding the proper limits, first note

Integrate[
PDF[DirichletDistribution[{1, 1, 1, 1}], {p1, p2, p3}],
{p1, 0, 1}, {p2, 0, 1}, {p3, 0, 1}
]
(* 1 *)

Now for the limits that respect the additional constraints (i.e., $p_1 > p_2$ and $p_1 > p_3$):

Integrate[
PDF[DirichletDistribution[{1, 1, 1, 1}], {p1, p2, p3}],
{p1, 0, 1}, {p2, 0, p1}, {p3, 0, p1}
]
(* 1/3 *)

We may run a simulation to check the reasonableness of this result:

Count[
RandomVariate[DirichletDistribution[{1, 1, 1, 1}], 10^5],
{p1_, p2_, p3_} /; p1 > p2 && p1 > p3
]/10^5 // N
(* 0.3307 *)

The analytical result seems reasonable.

This

RegionPlot3D[p1+p2+p3<=1 && p1>p2 && p1>p3 && p1>=0 && p2>=0 && p3>=0,
{p1,0,1}, {p2,0,1}, {p3,0,1}, PlotPoints->100, AxesLabel->{p1,p2,p3}]

shows the five sided prism that your viable points live inside. From that you can see the viable ranges for p2 and p3. I believe that p4 is implicitly represented in that diagram.

Assuming[x1>0 && x2>0 && x3>0 && x4>0 && b1>0 && b2>0 && b3>0 && b4>0,
Simplify[Integrate[Boole[p1+p2+p3<=1 && p1>p2 && p1>p3 && p1>=0 && p2>=0 && p3>=0]*
p1^(b1+x1-1),
{p1,0,1}, {p2,0,1/2}, {p3,0,1/2}]]]

As you incorporate your additional factors the result is bigger and much slower until it finally can't give you a simple closed form without additional information about your parameters.

Assuming[x1>0 && x2>0 && x3>0 && x4>0 && b1>0 && b2>0 && b3>0 && b4>0,
Simplify[Integrate[Boole[p1+p2+p3<=1 && p1>p2 && p1>p3 && p1>=0 && p2>=0 && p3>=0]*
p1^(b1+x1-1) p2^(b2+x2-1) p3^(b3+x3-1)(1-p1-p2-p3)^(x4+b4-1),
{p1,0,1}, {p2,0,1/2}, {p3,0,1/2}]]

If you incorporate the additional information about your parameters, perhaps even up to some or all of the numeric values, then I think this gives you the result that you are looking for.