How to calculate the dirichlet cdf

Dear Fellows,

Do you know how to write the cdf of dirichlet distribution? I could write the pdf but do not know how to formulate the cdf. I would like to calculate the area between x-0.5 to x+0.5, which is F(x+0.5)-F(x-0.5) where F refers to the cdf. Could you offer me some advice about this?

Thanks very much.

Xueqi

Answer 1

[Star Strider]

Wikipedia doesn’t give an analytic expression for the cdf. Since you already have a function written for the pdf, I suggest simply integrating it using the integral or integral2 function.

[xueqi]

Hey. Thanks. I tried to use integral2.

if true
  % a=[2,5,4];
X=[0.4,0.3,0.3];
Y=dirpdf(X,a)
f1=@(x)dirpdf(x,a);
f2=@(x1,x2) dirpdf([x1,x2,1-x1-x2],a);
f2(0.3,0.2);
f3=@(x1,x2) integral2(f2,x1-0.005,x1+0.005,x2-0.005,x2+0.005);
f3(0.3,0.2)

end

But I get the error message

if true
  % 
Error using dirpdf (line 46)
Columns in X must be equal or one less than columns in ALPHA
Error in @(x1,x2)dirpdf([x1,x2,1-x1-x2],a)
Error in integral2Calc>integral2t/tensor (line 229)
      Z = FUN(X,Y);  NFE = NFE + 1;
Error in integral2Calc>integral2t (line 56)
[Qsub,esub] = tensor(thetaL,thetaR,phiB,phiT);
Error in integral2Calc (line 10)
  [q,errbnd] = integral2t(fun,xmin,xmax,ymin,ymax,optionstruct);
Error in integral2 (line 107)
  Q = integral2Calc(fun,xmin,xmax,yminfun,ymaxfun,opstruct);
Error in @(x1,x2)integral2(f2,x1-0.005,x1+0.005,x2-0.005,x2+0.005)
Error in testdircdf (line 8)
f3(0.3,0.2)
end

The function dirpdf and f2 are working fine as I can get f2(0.3,0.2)=1.512. So there is something wrong when I use integral2. Could you find out where is wrong?

[Star Strider]

It is difficult for me to follow your code. I have no experience with dirpdf.

If you are using vector-valued functions, or if your call to dirpdf returns a vector, it might be necessary for you to use two nested integral calls, one to evaluate the value for x1 and the other for x2. The integral function will take vectors, while integral2 will not.

If that is not the problem, you will likely need to go through your code and your calls to dirpdf with the arguments you present to it in your call to integral or integral2 to understand what it returns.

If all else fails, generate a vector of values for x, evaluate dirpdf in as many nested loops as necessary, and then use trapz to integrate it.