Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
integrate a dirichlet pdf when n is not an integer

Subject: integrate a dirichlet pdf when n is not an integer

From: frankie tsoi

Date: 23 Jun, 2010 07:39:06

Message: 1 of 5

Hi,

I have generated a bunch of data with the dirichlet distribution and then fit by mle. The problem is the resulting estimators of v1, v2, and v3 are not integers. May I ask if anyone has any idea on how to integrate such a function when these v1, v2, and v3 are not integers on regions such as 0 < x <1/3 -y and 0 < y < 1/3?

the dirichlet distribution is given by: gamma(v1+v2+v3)/(gamma(v1)gamma(v2)gamma(v3))*x^(v1-1)*y^(v2-1)*(1-x-y)^(v3-1)

thanks a lot.

frankie

Subject: integrate a dirichlet pdf when n is not an integer

From: Peter Perkins

Date: 23 Jun, 2010 14:29:58

Message: 2 of 5

On 6/23/2010 3:39 AM, frankie tsoi wrote:
> Hi,
>
> I have generated a bunch of data with the dirichlet distribution and
> then fit by mle. The problem is the resulting estimators of v1, v2, and
> v3 are not integers. May I ask if anyone has any idea on how to
> integrate such a function when these v1, v2, and v3 are not integers on
> regions such as 0 < x <1/3 -y and 0 < y < 1/3?
>
> the dirichlet distribution is given by:
> gamma(v1+v2+v3)/(gamma(v1)gamma(v2)gamma(v3))*x^(v1-1)*y^(v2-1)*(1-x-y)^(v3-1)

Nothing in that expression requires an integer, although you are
obviously missing some *'s. But you almost certainly want to rewrite
your PDF by exponentiating a sum of gammaln terms and terms like
(v1-1)*log(x). And you will need to use ".*".

Hope this helps.

Subject: integrate a dirichlet pdf when n is not an integer

From: Roger Stafford

Date: 23 Jun, 2010 22:30:37

Message: 3 of 5

"frankie tsoi" <frank_tsoi@hotmail.com> wrote in message <hvsdmq$jb1$1@fred.mathworks.com>...
> Hi,
>
> I have generated a bunch of data with the dirichlet distribution and then fit by mle. The problem is the resulting estimators of v1, v2, and v3 are not integers. May I ask if anyone has any idea on how to integrate such a function when these v1, v2, and v3 are not integers on regions such as 0 < x <1/3 -y and 0 < y < 1/3?
>
> the dirichlet distribution is given by: gamma(v1+v2+v3)/(gamma(v1)gamma(v2)gamma(v3))*x^(v1-1)*y^(v2-1)*(1-x-y)^(v3-1)
>
> thanks a lot.
>
> frankie
- - - - - - - - - -
  For the particular triangle you have described, 0<=x, 0<=y, x+y<=1/3, there is no need to do numerical integration. You can find the integral of the dirichlet pdf over it directly using matlab's betainc function.

  The density function is

 f(x,y) = gamma(a+b+c)/gamma(a)/gamma(b)/gamma(c) * ...
          x^(a-1)*y^(b-1)*(1-x-y)^(c-1)

For convenience, temporarily drop the gamma factors and define a change of variables: t = x + y, u = y. The Jacobian is 1 and the integrand becomes

 (t-u)^(a-1)*u^(b-1)*(1-t)^(c-1)

so the integral over the given triangle can be expressed as the iterated integrals:

 int( (1-t)^(c-1) * int( (t-u)^(a-1)*u^(b-1),'u',0,t),'t',0,1/3)

Then substitute v = u/t getting du = t*dv and

 int( (1-t)^(c-1)*int( (t-v*t)^(a-1)*(v*t)^(b-1)*t,'v',0,1),'t',0,1/3) =
 int( (1-t)^(c-1)*t^(a+b-1) * int( (1-v)^(a-1)*v^(b-1),0,1),'t',0,1/3) =
 gamma(a)*gamma(b)/gamma(a+b)* int( (1-t)^(c-1)*t^(a+b-1),'t',0,1/3)
 
If we bring back the original gamma factors, this gives us

 gamma(a+b+c)/gamma(a+b)/gamma(c)*int( t^(a+b-1)*(1-t)^(c-1),'t',0,1/3) =
 betainc(1/3,a+b,c)

as the final integral value over the given triangle. Therefore all you need to do is make the single call:

 betainc(1/3,a+b,c)

and you have your integral over the triangle.

  Of course the value 1/3 can actually be any value between 0 and 1. Note that at 1, you will get 1 as the integral value, as you would require for a valid pdf.

Roger Stafford

Subject: integrate a dirichlet pdf when n is not an integer

From: frankie tsoi

Date: 24 Jun, 2010 14:22:09

Message: 4 of 5

Peter Perkins <Peter.Perkins@MathRemoveThisWorks.com> wrote in message <hvt5p6$neo$1@fred.mathworks.com>...
> On 6/23/2010 3:39 AM, frankie tsoi wrote:
> > Hi,
> >
> > I have generated a bunch of data with the dirichlet distribution and
> > then fit by mle. The problem is the resulting estimators of v1, v2, and
> > v3 are not integers. May I ask if anyone has any idea on how to
> > integrate such a function when these v1, v2, and v3 are not integers on
> > regions such as 0 < x <1/3 -y and 0 < y < 1/3?
> >
> > the dirichlet distribution is given by:
> > gamma(v1+v2+v3)/(gamma(v1)gamma(v2)gamma(v3))*x^(v1-1)*y^(v2-1)*(1-x-y)^(v3-1)
>
> Nothing in that expression requires an integer, although you are
> obviously missing some *'s. But you almost certainly want to rewrite
> your PDF by exponentiating a sum of gammaln terms and terms like
> (v1-1)*log(x). And you will need to use ".*".
>
> Hope this helps.

Thanks a lot, Peter.

Subject: integrate a dirichlet pdf when n is not an integer

From: frankie tsoi

Date: 24 Jun, 2010 14:30:26

Message: 5 of 5

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <hvu1ud$71b$1@fred.mathworks.com>...
> "frankie tsoi" <frank_tsoi@hotmail.com> wrote in message <hvsdmq$jb1$1@fred.mathworks.com>...
> > Hi,
> >
> > I have generated a bunch of data with the dirichlet distribution and then fit by mle. The problem is the resulting estimators of v1, v2, and v3 are not integers. May I ask if anyone has any idea on how to integrate such a function when these v1, v2, and v3 are not integers on regions such as 0 < x <1/3 -y and 0 < y < 1/3?
> >
> > the dirichlet distribution is given by: gamma(v1+v2+v3)/(gamma(v1)gamma(v2)gamma(v3))*x^(v1-1)*y^(v2-1)*(1-x-y)^(v3-1)
> >
> > thanks a lot.
> >
> > frankie
> - - - - - - - - - -
> For the particular triangle you have described, 0<=x, 0<=y, x+y<=1/3, there is no need to do numerical integration. You can find the integral of the dirichlet pdf over it directly using matlab's betainc function.
>
> The density function is
>
> f(x,y) = gamma(a+b+c)/gamma(a)/gamma(b)/gamma(c) * ...
> x^(a-1)*y^(b-1)*(1-x-y)^(c-1)
>
> For convenience, temporarily drop the gamma factors and define a change of variables: t = x + y, u = y. The Jacobian is 1 and the integrand becomes
>
> (t-u)^(a-1)*u^(b-1)*(1-t)^(c-1)
>
> so the integral over the given triangle can be expressed as the iterated integrals:
>
> int( (1-t)^(c-1) * int( (t-u)^(a-1)*u^(b-1),'u',0,t),'t',0,1/3)
>
> Then substitute v = u/t getting du = t*dv and
>
> int( (1-t)^(c-1)*int( (t-v*t)^(a-1)*(v*t)^(b-1)*t,'v',0,1),'t',0,1/3) =
> int( (1-t)^(c-1)*t^(a+b-1) * int( (1-v)^(a-1)*v^(b-1),0,1),'t',0,1/3) =
> gamma(a)*gamma(b)/gamma(a+b)* int( (1-t)^(c-1)*t^(a+b-1),'t',0,1/3)
>
> If we bring back the original gamma factors, this gives us
>
> gamma(a+b+c)/gamma(a+b)/gamma(c)*int( t^(a+b-1)*(1-t)^(c-1),'t',0,1/3) =
> betainc(1/3,a+b,c)
>
> as the final integral value over the given triangle. Therefore all you need to do is make the single call:
>
> betainc(1/3,a+b,c)
>
> and you have your integral over the triangle.
>
> Of course the value 1/3 can actually be any value between 0 and 1. Note that at 1, you will get 1 as the integral value, as you would require for a valid pdf.
>
> Roger Stafford

Thanks a lot, Roger. I ended up using the NIntegrate function with mathematica summing up tiny retangular regions. I will try this betainc function and see if the result is more accurate. Thanks again.

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us