would you mind to post that example? evalf(Int(f,t=0..1)) uses the NAG library i think ...

randomize():
B:=[stats[random,normald[0,1]](N+1)]:
W:=[stats[random,normald[0,1]](N+1)]:
Z:= zip((x,y) -> evalf( rho*x+sqrt(1-rho^2)*y ),B,W): # or evalhf?

Then Z and B will do, where rho is the desired correlation. In Maple 10 it is a bit different, but much (!) faster.

Once i translated Ooura's double exponential integration into Maple (find it uploaded, but i have neither cleaned up that port nor tested it carefully). That method is quite robust. So you can use Maple's new ability to generate a DLL and use that as a function directly in Maple. This might give you some of the speed improvements you are looking for. The additional advantage is: you can test values against Maple if you have some doubts. For distributions usually a Gauss-Legendre method also works well (using fixed numbers of sub-divisions), 16 point suffice. What i can not see is your chi^2, i thought it contains an exp in the integrand? Just a last thought: when i try to fit some pdf against data (to get parameters) i pass to the logarithm (and refine later if neccessary). If you want to fit against a cdf then i can imagine troubles (i am not sure whether this is a stable approach, especially if you are not absolutely convinced that your cdf is the correct one). PS: as i do not recognize how to attach files use that link http://www.axelvogt.de/axalom/maple/Ooura_DoubleExp_Integration.mws

the only thing i can imagine is that Maple methods can be worked out ... having not traced through the procs it is clear that they involve external callings for NAG routines what exactly are you looking for?