1
$\begingroup$

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?

$\endgroup$
2
  • 1
    $\begingroup$ 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 $\endgroup$ Jun 19 '17 at 9:08
  • $\begingroup$ Altho mef has already informed you about DirichletDistribution[], you might nevertheless be interested in looking up CylindricalDecomposition[] for determining your integration bounds. $\endgroup$
    – J. M.'s torpor
    Aug 14 '17 at 14:44
1
$\begingroup$

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.

$\endgroup$
4
  • $\begingroup$ 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.? $\endgroup$
    – ProgSnob
    Jun 19 '17 at 10:42
  • $\begingroup$ 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. $\endgroup$
    – TeM
    Jun 19 '17 at 10:46
  • $\begingroup$ $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$ $\endgroup$
    – ProgSnob
    Jun 19 '17 at 11:01
  • $\begingroup$ 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. $\endgroup$
    – TeM
    Jun 19 '17 at 12:02
1
$\begingroup$

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.

$\endgroup$
0
$\begingroup$

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.

This then does your integral, but only including your first factor.

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.