-----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.)
(%i1) (load(abs_integrate), load(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.
(%i11) g(a,b,c,x_2) := ''(radcan(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