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.

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.

Hm. Two "whoopses", really. First was that my previous comment seems to show up with the html codes escaped - I'm not sure why, and I can't seem to get it fixed.

Second is that there may be an issue in the code I posted above. I tried to apply the fully-numerical approach to the formula that the original poster posted, as follows:

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;

This leads to different results than the procedure I posted above, though the same qualitative result arises. I get g(1) = 0.087, g(10) = 0.149, g(100) = 0.162, g(1000)=0.164. I'll file an SCR about and look at it.

Erik Postma
Maplesoft.

Hm. Two "whoopses", really. First was that my previous comment seems to show up with the html codes escaped - I'm not sure why, and I can't seem to get it fixed.

Second is that there may be an issue in the code I posted above. I tried to apply the fully-numerical approach to the formula that the original poster posted, as follows:

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;

This leads to different results than the procedure I posted above, though the same qualitative result arises. I get g(1) = 0.087, g(10) = 0.149, g(100) = 0.162, g(1000)=0.164. I'll file an SCR about and look at it.

Erik Postma
Maplesoft.

That's a great solution, Scott. Here's an extension that also includes the minimization aspect.

What we're looking for is the minimal positive root of that right hand side. We can find that using the RootFinding:-NextZero command, which will (numerically) find the first root of a function to the right of a given point. So we need to translate the solution that we have into a function. Scott already showed how to translate it into an expression; to go from an expression to a function, we'll use the unapply function. That leads to the following, after defining X1 as before:

theFunction := unapply(rhs(X1), t);
RootFinding:-NextZero(theFunction, 0);

The answer that Maple finds is 0.1002991917. Looking at the plot, I think that's possibly the only zero, so in this case there's not really any difference between the two solutions. But if there are multiple zeros, this might work better.

Hope this helps,

Erik Postma
Maplesoft.

That's a great solution, Scott. Here's an extension that also includes the minimization aspect.

What we're looking for is the minimal positive root of that right hand side. We can find that using the RootFinding:-NextZero command, which will (numerically) find the first root of a function to the right of a given point. So we need to translate the solution that we have into a function. Scott already showed how to translate it into an expression; to go from an expression to a function, we'll use the unapply function. That leads to the following, after defining X1 as before:

theFunction := unapply(rhs(X1), t);
RootFinding:-NextZero(theFunction, 0);

The answer that Maple finds is 0.1002991917. Looking at the plot, I think that's possibly the only zero, so in this case there's not really any difference between the two solutions. But if there are multiple zeros, this might work better.

Hope this helps,

Erik Postma
Maplesoft.