It turns out that most of the time is spent in overhead in the Statistics package. There is fast external C code for shuffling, but finding that code takes a while. If you use the debugger to step through a Statistics[Shuffle] call, you can see that in the end, it calls an external function called MapleShuffle. You can shave off another factor of 3 or so by calling that directly. The following example was directly adapted from your code; it may of course not work in past or future versions of Maple (it does work in Maple 14, at least).

restart;
with(Statistics):
opacity := kernelopts(opaquemodules=false):
Statistics:-ExternalSupport:-initializeHwLibrary():
shuffle := Statistics:-ExternalSupport:-DefineExternal('MapleShuffle');
kernelopts(opaquemodules=opacity):
interface(rtablesize=12):
L := Array([2, 2, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0]):
Y := Array(1..6,[seq(0,i=2..7)]):
X := Array(1..6,[seq(0,i=2..7)]):
n:=792000:
# Digits:=2: # why?
st:=time():
for i to n do
L := shuffle(L):
E1:=add(L[j],j=1..5):
X[E1-1] := X[E1-1]+1:
E2:=add(L[j],j=6..10):
Y[E1-1] := Y[E1-1] + E2:
end do:
time()-st;

Hope this helps,

Erik Postma
Maplesoft.

It turns out that most of the time is spent in overhead in the Statistics package. There is fast external C code for shuffling, but finding that code takes a while. If you use the debugger to step through a Statistics[Shuffle] call, you can see that in the end, it calls an external function called MapleShuffle. You can shave off another factor of 3 or so by calling that directly. The following example was directly adapted from your code; it may of course not work in past or future versions of Maple (it does work in Maple 14, at least).

restart;
with(Statistics):
opacity := kernelopts(opaquemodules=false):
Statistics:-ExternalSupport:-initializeHwLibrary():
shuffle := Statistics:-ExternalSupport:-DefineExternal('MapleShuffle');
kernelopts(opaquemodules=opacity):
interface(rtablesize=12):
L := Array([2, 2, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0]):
Y := Array(1..6,[seq(0,i=2..7)]):
X := Array(1..6,[seq(0,i=2..7)]):
n:=792000:
# Digits:=2: # why?
st:=time():
for i to n do
L := shuffle(L):
E1:=add(L[j],j=1..5):
X[E1-1] := X[E1-1]+1:
E2:=add(L[j],j=6..10):
Y[E1-1] := Y[E1-1] + E2:
end do:
time()-st;

Hope this helps,

Erik Postma
Maplesoft.

Thanks Alec. I always do that, and it always shows up perfectly for me in the preview, using either the HTML options. Then in the posted comment itself, boom, it's gone. Maybe there's a gremlin in the system that picks on my browser/OS combination or something.

Erik.

Thanks Alec. I always do that, and it always shows up perfectly for me in the preview, using either the HTML options. Then in the posted comment itself, boom, it's gone. Maybe there's a gremlin in the system that picks on my browser/OS combination or something.

Erik.

Here's one option, that works especially pleasantly with Acer's evalhf'able integrand (thanks Acer) :

n := 7;
plot(numericIntegral, 3 .. 7);

Now it may look like there's a zero at t = 5, but if you make Maple look closely, you can see that there's a reason why fsolve won't tell you that:

plot(numericIntegral, 0 .. 5);

It shows you that the integral is positive but extremely small, probably all due to roundoff error.

OK - I can't help myself, I have to give you an option to do it non-graphically :) If you're interested in just finding a spot where the integral is close to zero, you might want to use some custom code (instead of fsolve or RootFinding[NextZero]). For example, here's a very simple piece of code (using bisection) that finds the maximum value (up to an x tolerance, xtol) in a range where the supplied function is within a tolerance (ytol) of zero - or at least, a value that is a right end point of a range where the supplied function is near zero, there could be other such ranges further to the right:

maxApproximateRoot := proc(f :: procedure, r :: range(numeric), xtol :: numeric, ytol :: numeric, $)
local a, b, c, fa, fb, fc;
  a, b := op(evalf(r));
  fb := f(b);
  if abs(fb) <= ytol then
    return b;
  end if;

  fa := f(a);
  if fa * fb > 0 and abs(fa) > ytol then
    error "f has the same sign on both sides of the supplied interval and is not close to zero on either end - aborting";
  end if;

  # Invariant: (fa * fb < 0 or abs(fa) <= ytol) and abs(fb) > ytol
  while b - a > xtol do
    c := (a + b)/2;
    fc := f(c);

    if abs(fc) <= ytol or fc * fb < 0 then
      a := c;
      fa := fc;
    else
      b := c;
      fb := fc;
    end if;
  end do;

  return (a + b)/2;
end proc;

You could use it as follows:

maxApproximateRoot(numericIntegral, 3 .. 7, 1e-5, 1e-5);
# returns 5.019527435 (for n = 7)

Indeed, running

plot(numericIntegral, 5.01 .. 5.03);

gives a very interesting, non-continuous plot... numeric issues again!

Hope this helps,

Erik Postma
Maplesoft.

Here's one option, that works especially pleasantly with Acer's evalhf'able integrand (thanks Acer) :

n := 7;
plot(numericIntegral, 3 .. 7);

Now it may look like there's a zero at t = 5, but if you make Maple look closely, you can see that there's a reason why fsolve won't tell you that:

plot(numericIntegral, 0 .. 5);

It shows you that the integral is positive but extremely small, probably all due to roundoff error.

OK - I can't help myself, I have to give you an option to do it non-graphically :) If you're interested in just finding a spot where the integral is close to zero, you might want to use some custom code (instead of fsolve or RootFinding[NextZero]). For example, here's a very simple piece of code (using bisection) that finds the maximum value (up to an x tolerance, xtol) in a range where the supplied function is within a tolerance (ytol) of zero - or at least, a value that is a right end point of a range where the supplied function is near zero, there could be other such ranges further to the right:

maxApproximateRoot := proc(f :: procedure, r :: range(numeric), xtol :: numeric, ytol :: numeric, $)
local a, b, c, fa, fb, fc;
  a, b := op(evalf(r));
  fb := f(b);
  if abs(fb) <= ytol then
    return b;
  end if;

  fa := f(a);
  if fa * fb > 0 and abs(fa) > ytol then
    error "f has the same sign on both sides of the supplied interval and is not close to zero on either end - aborting";
  end if;

  # Invariant: (fa * fb < 0 or abs(fa) <= ytol) and abs(fb) > ytol
  while b - a > xtol do
    c := (a + b)/2;
    fc := f(c);

    if abs(fc) <= ytol or fc * fb < 0 then
      a := c;
      fa := fc;
    else
      b := c;
      fb := fc;
    end if;
  end do;

  return (a + b)/2;
end proc;

You could use it as follows:

maxApproximateRoot(numericIntegral, 3 .. 7, 1e-5, 1e-5);
# returns 5.019527435 (for n = 7)

Indeed, running

plot(numericIntegral, 5.01 .. 5.03);

gives a very interesting, non-continuous plot... numeric issues again!

Hope this helps,

Erik Postma
Maplesoft.

Hi!

There seem to be fairly severe numerical difficulties that arise in evaluating the integral. If I experiment for a bit, I can see that there are values where the integrand is almost exactly zero for most of the domain t - 5 .. t + 5, and clearly nonzero for another part. It may be that this causes the difficulties, or it could be something else. One option to make the situation a bit more robust is to change the numericIntegral procedure to set the epsilon and digits options explicitly; for example,

numericIntegral := proc(t)
local d, result;
  for d in [0, 2, 5, 10] do
    result := evalf[Digits](Int(integrand(t), t-5 .. t+5, digits = Digits + d));
    if type(result, float) then
      return result;
    end if;
  end do;

  return result;
end proc;

Another thing that would help would be, if you could say beforehand what the range will be where the integrand is zero, and restrict the integral to the nonzero part of the range. For example, with n = 9 and t = -4, it looks like the integrand is nonzero only between 0 and 1, as you can see by running n := 9; plot(integrand(-4), -9 .. 1); But this information may be tricky to come by.

I also noticed that the ceil(10 - x/5) - 2 can be rewritten to 8 + ceil(-x/5); or did you mean ceil((10 - x)/5) - 2?

Finally, I've had some success where I failed otherwise, and some failures where I had successes otherwise, replacing the Sum in integrand by add, and as a third option, just solving the symbolic sum, as follows:

symsum := x*(1 - eval(sum(binomial(n, k) * ((x-t+5)/10)^k * (1-((x-t+5)/10))^(n-k), k=0..K), K = ceil((10 - x)/5) - 2));
symsumfun := unapply(symsum, t, x);
integrand := t -> (x -> symsumfun(t, x));

Hope that some of these work for you,

Erik Postma
Maplesoft.

Hi!

There seem to be fairly severe numerical difficulties that arise in evaluating the integral. If I experiment for a bit, I can see that there are values where the integrand is almost exactly zero for most of the domain t - 5 .. t + 5, and clearly nonzero for another part. It may be that this causes the difficulties, or it could be something else. One option to make the situation a bit more robust is to change the numericIntegral procedure to set the epsilon and digits options explicitly; for example,

numericIntegral := proc(t)
local d, result;
  for d in [0, 2, 5, 10] do
    result := evalf[Digits](Int(integrand(t), t-5 .. t+5, digits = Digits + d));
    if type(result, float) then
      return result;
    end if;
  end do;

  return result;
end proc;

Another thing that would help would be, if you could say beforehand what the range will be where the integrand is zero, and restrict the integral to the nonzero part of the range. For example, with n = 9 and t = -4, it looks like the integrand is nonzero only between 0 and 1, as you can see by running n := 9; plot(integrand(-4), -9 .. 1); But this information may be tricky to come by.

I also noticed that the ceil(10 - x/5) - 2 can be rewritten to 8 + ceil(-x/5); or did you mean ceil((10 - x)/5) - 2?

Finally, I've had some success where I failed otherwise, and some failures where I had successes otherwise, replacing the Sum in integrand by add, and as a third option, just solving the symbolic sum, as follows:

symsum := x*(1 - eval(sum(binomial(n, k) * ((x-t+5)/10)^k * (1-((x-t+5)/10))^(n-k), k=0..K), K = ceil((10 - x)/5) - 2));
symsumfun := unapply(symsum, t, x);
integrand := t -> (x -> symsumfun(t, x));

Hope that some of these work for you,

Erik Postma
Maplesoft.

Note that if you had StringTools loaded, you would get StringTools[Select] if you use uppercase S - see ?StringTools,Select - whereas the lowercase should probably always give you the "regular" select - see ?select.

Hope this helps,

Erik Postma
Maplesoft.

Note that if you had StringTools loaded, you would get StringTools[Select] if you use uppercase S - see ?StringTools,Select - whereas the lowercase should probably always give you the "regular" select - see ?select.

Hope this helps,

Erik Postma
Maplesoft.

Yes, that helps a lot. I assume you want to know this value of t for a specific value of n? I'll show here for n = 45 first.

Looking at what you write, maybe the first thing to realize is that Maple uses [ and ] to indicate lists, not algebraic expressions. Also, I found that for this case it seems to work better to define the integrand as a procedure. Because I wanted to be able to vary t as well, I decided to use a sort of two-step approach - see below.

integrand := proc(t)
  return proc(x)
  local k;
    return x * (1 - Sum(binomial(n, k) * ((x - t + 5)/10)^k * (1 - ((x - t + 5)/10))^(n-k),
                        k = 0 .. ceil(10 - x/5) - 2));
  end proc;
end proc;
n := 45;

This allows you to write integrand(10)(7) to obtain the integrand with t = 10 and x = 7 (or evalf(integrand(10)(7)) for a numerical value); you get a procedure for the integrand (as a function of x) by writing e.g. integrand(10) for t = 10. Note furthermore that I used Sum instead of sum this time, since I don't expect anything useful to come out of an attempt to process the summation symbolically. For the same reason, we will use Int instead of int for the integral, then run evalf to obtain a numerical value, as follows:

numericIntegral := t -> evalf(Int(integrand(t), t-5 .. t+5));

Recall that integrand(t) is itself a procedure, not an expression, so we supply the range without giving a variable name. Now we can evaluate numericIntegral(10) to obtain the numerical value of the integral for t = 10. Finally, the last piece of the puzzle is to run the command fsolve to make Maple search for a value where the numericIntegral is zero. This takes about a minute or so, but it will find the value t = -1.064932177. If you call trace(numericIntegral) beforehand, you can see that most of the time is spent refining an approximate solution to the full number of digits you need; you can lower the setting of Digits to indicate that you don't need too many digits of precision, but it will still do some conservative checking of the solution to make sure that those digits are really, really correct. Note that setting Digits too low may not work very well. So, to tie it all together:

# either:
fsolve(numericIntegral);
# or (with more information):
trace(numericIntegral);
fsolve(numericIntegral);
# or (faster):
Digits := 6;
fsolve(numericIntegral);

Now you can do this for a different value of n by reassigning, say,

n := 100;

and re-evaluating the fsolve command.

Hope this helps,

Erik Postma
Maplesoft.

Yes, that helps a lot. I assume you want to know this value of t for a specific value of n? I'll show here for n = 45 first.

Looking at what you write, maybe the first thing to realize is that Maple uses [ and ] to indicate lists, not algebraic expressions. Also, I found that for this case it seems to work better to define the integrand as a procedure. Because I wanted to be able to vary t as well, I decided to use a sort of two-step approach - see below.

integrand := proc(t)
  return proc(x)
  local k;
    return x * (1 - Sum(binomial(n, k) * ((x - t + 5)/10)^k * (1 - ((x - t + 5)/10))^(n-k),
                        k = 0 .. ceil(10 - x/5) - 2));
  end proc;
end proc;
n := 45;

This allows you to write integrand(10)(7) to obtain the integrand with t = 10 and x = 7 (or evalf(integrand(10)(7)) for a numerical value); you get a procedure for the integrand (as a function of x) by writing e.g. integrand(10) for t = 10. Note furthermore that I used Sum instead of sum this time, since I don't expect anything useful to come out of an attempt to process the summation symbolically. For the same reason, we will use Int instead of int for the integral, then run evalf to obtain a numerical value, as follows:

numericIntegral := t -> evalf(Int(integrand(t), t-5 .. t+5));

Recall that integrand(t) is itself a procedure, not an expression, so we supply the range without giving a variable name. Now we can evaluate numericIntegral(10) to obtain the numerical value of the integral for t = 10. Finally, the last piece of the puzzle is to run the command fsolve to make Maple search for a value where the numericIntegral is zero. This takes about a minute or so, but it will find the value t = -1.064932177. If you call trace(numericIntegral) beforehand, you can see that most of the time is spent refining an approximate solution to the full number of digits you need; you can lower the setting of Digits to indicate that you don't need too many digits of precision, but it will still do some conservative checking of the solution to make sure that those digits are really, really correct. Note that setting Digits too low may not work very well. So, to tie it all together:

# either:
fsolve(numericIntegral);
# or (with more information):
trace(numericIntegral);
fsolve(numericIntegral);
# or (faster):
Digits := 6;
fsolve(numericIntegral);

Now you can do this for a different value of n by reassigning, say,

n := 100;

and re-evaluating the fsolve command.

Hope this helps,

Erik Postma
Maplesoft.

Yes. Sorry.

I generally try not to make procedures depend on what's loaded with "with", because then the user can decide what to load for themselves, and the best way to do that without typing way too much is using "uses". I should avoid writing replies when it's late and I'm tired though... it causes me to forget to copy and paste the code into a fresh session to see if I haven't screwed up :)

I fixed the post above to address the issues that Scott and acer uncovered. Sorry again for the clumsiness.

Erik Postma
Maplesoft.

Yes. Sorry.

I generally try not to make procedures depend on what's loaded with "with", because then the user can decide what to load for themselves, and the best way to do that without typing way too much is using "uses". I should avoid writing replies when it's late and I'm tired though... it causes me to forget to copy and paste the code into a fresh session to see if I haven't screwed up :)

I fixed the post above to address the issues that Scott and acer uncovered. Sorry again for the clumsiness.

Erik Postma
Maplesoft.

I'll try and at least make the code show up correctly here:

f := proc(n)                                             
uses Statistics = ST;                                    
local quantile, X, Y;                                    
    X := ST:-RandomVariable(FRatio(3, 3*n));                
    Y := ST:-RandomVariable(NonCentralFRatio(3, 3*n, 0.82));
    quantile := ST:-Quantile(X, 0.95, numeric);             
    return Probability(Y > quantile, numeric);              
end proc;    

g := proc(n)
uses Statistics = ST;
local X, Y, i;
  X := ST:-RandomVariable(FRatio(3, 3*n));
  return evalf(Sum(exp(-0.82)*0.82^i/i! *
      ST:-Probability(
        ST:-RandomVariable(FRatio(3+2*i, 3*n)) >
        3*ST:-Quantile(FRatio(3, 3*n), 0.95, numeric) / (3 + 2*i),
        numeric),
      i=0..infinity));
end proc;

Hope this works...

Erik.