# [Maxima] Integrate Dirichlet distribution

Barton Willis willisb at unk.edu
Sat Jun 5 06:16:57 CDT 2010

```-----maxima-bounces at math.utexas.edu wrote: -----

>dirichlet_dist = x_1**(a-1) * x_2**(b-1) * (1-x_2-x_1)**(c-1)
>e = integral(dirichlet_dist, x_1, 0, 1)
>f = integral(e, x_2, 0, 1)

The optional package hyperint can evaluate this definite integral. Also,
abs_integrate provides a mechanism for extending Maxima's integrate function. (Pay no
attention to the warnings Maxima prints when it loads hyperint.)

(%i2) extra_integration_methods : cons(hyper_int, extra_integration_methods)\$

(%i3) f : integrate(x_1**(a-1) * x_2**(b-1) * (1-x_2-x_1)**(c-1),x_1,0,1);
"Is  "a"  positive, negative, or zero?"pos;
(%o3) -(hypergeometric([a,1-c],[a+1],-1/(x_2-1))*x_2^(b-1)*%e^(c*log(-x_2)))/((x_2/(x_2-1))^c*(a*x_2-a))

Check a special case:

(%i4) subst([c=5,a=1,x_2=1/2,b=1/3], x_1**(a-1) * x_2**(b-1) * (1-x_2-x_1)**(c-1));
(%o4) 2^(2/3)*(1/2-x_1)^4

(%i5) integrate(%,x_1,0,1);
(%o5) 1/(5*2^(10/3))

(%i6) float(%);
(%o6) 0.019842513149602

Checks OK:

(%i7) subst([c=5,a=1,x_2=1/2,b=1/3],f);
(%o7) 1/(5*2^(10/3))

Floating point evaluation can be a problem :( This gives a spurious imaginary part,
and requires a call to rectform to fully evaluate to a float:

(%i8) g(a,b,c,x_2) := ''f\$
(%i9) g(1,1/3,5,0.5);

(%o9) -0.63496042078728*%e^(5*(3.141592653589793*%i-0.69314718055995))

(%i10) rectform(%);
(%o10) 0.019842513149602-1.2149633839403642*10^-17*%i

A workaround is to apply radcan to f.

OK:

(%i12) g(1,1/3,5,0.5);
(%o12) 0.019842513149603

And this is a bug in the hypergeometric code, I think:

(%i13) g(1.0,1/3.0,5.0,0.5);
CLIENT: Lost socket connection ...
Restart Maxima with 'Maxima->Restart Maxima'.

Sorry about that (it's my fault)--if I'm able, I'll fix it.

--Barton

```