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.

I know that your issue has already been solved by the other replies, but just for completeness, here is a possible solution involving the taylor command. As you already found out, the issue is that what you've written doesn't distinguish between the general variable, x, and the specific value for it, 3, once Maple enters the procedure f that you've written. One way to solve that is by using a variable instead of the value when you call f:

(**) expr := f(a);
                                            2    3
                          expr :=  1 + a + a  + a

(**) eval(%, a = 3);
                                      40

This gives you access to the formal expression, expr, first, which you can then evaluate at a point. If all you're interested in is the evaluation, you could use a local variable in f; then, you'll need to write it explicitly as a procedure instead of using arrow notation:

newF := proc(a)
local x;
  return eval(convert(taylor(1/(1-x), x=0, 4), polynom), x = a);
end proc;

This way, you know that x will just be a variable (with no value assigned to it) when taylor gets to do its thing. Now newF(3) will return 40 directly. (And actually, the convert/polynom step is unnecessary now, because evaluating a series at a point drops the big-Oh term.) Note that newF(a) still returns the formal expression.

Hope this helps,

Erik Postma
Maplesoft.

I know that your issue has already been solved by the other replies, but just for completeness, here is a possible solution involving the taylor command. As you already found out, the issue is that what you've written doesn't distinguish between the general variable, x, and the specific value for it, 3, once Maple enters the procedure f that you've written. One way to solve that is by using a variable instead of the value when you call f:

(**) expr := f(a);
                                            2    3
                          expr :=  1 + a + a  + a

(**) eval(%, a = 3);
                                      40

This gives you access to the formal expression, expr, first, which you can then evaluate at a point. If all you're interested in is the evaluation, you could use a local variable in f; then, you'll need to write it explicitly as a procedure instead of using arrow notation:

newF := proc(a)
local x;
  return eval(convert(taylor(1/(1-x), x=0, 4), polynom), x = a);
end proc;

This way, you know that x will just be a variable (with no value assigned to it) when taylor gets to do its thing. Now newF(3) will return 40 directly. (And actually, the convert/polynom step is unnecessary now, because evaluating a series at a point drops the big-Oh term.) Note that newF(a) still returns the formal expression.

Hope this helps,

Erik Postma
Maplesoft.

Hi,

It looks like you're better off making dTvalues a list instead of a vector. Replace your definition by

dTvalues := [seq(seq([d, T], d=-30..30, 10), T=1..20)];

and then you should use nops instead of rtable_dims (since dTvalues is not an rtable anymore), so that would look like

i=1..nops(dTvalues)

where you had

i=rtable_dims(dTvalues)

before.

The other worksheet looks like quite a different problem and I'm not sure I have the time now to look at it. Since it's a different problem, you could try posting a new post for that and seeing if someone else can help you with that.

Hope this helps,

Erik Postma
Maplesoft.

Hi,

It looks like you're better off making dTvalues a list instead of a vector. Replace your definition by

dTvalues := [seq(seq([d, T], d=-30..30, 10), T=1..20)];

and then you should use nops instead of rtable_dims (since dTvalues is not an rtable anymore), so that would look like

i=1..nops(dTvalues)

where you had

i=rtable_dims(dTvalues)

before.

The other worksheet looks like quite a different problem and I'm not sure I have the time now to look at it. Since it's a different problem, you could try posting a new post for that and seeing if someone else can help you with that.

Hope this helps,

Erik Postma
Maplesoft.

Or, as a slight improvement, in this case you can also use ?LinearFit , since the expression is linear in the parameters (p[1] and p[2]). That will generally lead to a more accurate result and use less time. In fact, if you use ?Statistics[Fit] , then Maple will decide whether LinearFit is applicable or not; you could call it like this, for example:

independentDataMatrix := <<IC50> | <REDOX>>;
H := Statistics[Fit](p1 * x1 + p2 * x2 + p3, independentDataMatrix, LIPO, [x1, x2]);

You can add the output=solutionmodule option, like HerClau did above, if you like. Otherwise it will just give you the fitted function.

You need to make a decision which quantities you consider dependent and which independent; mathematically, there's no easy way around that.

Hope this helps,

Erik Postma
Maplesoft.

Or, as a slight improvement, in this case you can also use ?LinearFit , since the expression is linear in the parameters (p[1] and p[2]). That will generally lead to a more accurate result and use less time. In fact, if you use ?Statistics[Fit] , then Maple will decide whether LinearFit is applicable or not; you could call it like this, for example:

independentDataMatrix := <<IC50> | <REDOX>>;
H := Statistics[Fit](p1 * x1 + p2 * x2 + p3, independentDataMatrix, LIPO, [x1, x2]);

You can add the output=solutionmodule option, like HerClau did above, if you like. Otherwise it will just give you the fitted function.

You need to make a decision which quantities you consider dependent and which independent; mathematically, there's no easy way around that.

Hope this helps,

Erik Postma
Maplesoft.

Hi emoh,

I think this is a case where we could possibly do better with documentation, explaining what the intent is, but I think that it does perform as intended. You wrote:

About VolumeOfRevolution it says:

"The VolumeOfRevolution(f(x), x=a..b) command returns the volume of revolution of the expression f(x) from a to b."

Which is not true.

If you plot
VolumeOfRevolution(2*x^2,x=0..2,axis=vertical,'output'='plot), the program instead returns the volume of revolution of 2*x^2 from y=0..8.

However, if you leave out the axis option (axis = vertical) or replace it with the default, axis = horizontal, then I think you get what you would like to get. (If not, feel free to reply again and correct me.) The axis option is explained further down on that page, ?Student,Calculus1,VolumeOfRevolution , as follows:

axis = horizontal or vertical Whether the expression is rotated horizontally or vertically. By default, the rotation is horizontal.

The explanation here is not very good - what is meant is that the axis of rotation is horizontal or vertical, respectively, not the direction in which the plot is swept around. (As an aside, it's also of course not the expression that is rotated, but the curve.) So if you add the axis=vertical option, then the curve which runs from x=0 to x=2 (and from y=0 to y=8) is rotated around the y-axis to obtain the "outer surface" of the volume described. Indeed, if you plot(2*x^2, x=0..2) and mentally rotate the resulting curve around either the horizontal axis (y=0) or the vertical axis (x=0), you get what I get if I call VolumeOfRevolution(2*x^2, x=0..2, axis=horizontal, output=plot) or the same with axis = vertical, respectively.

For the second point you make, I would say the same thing is happening: what you describe indicates to me that you want to rotate around the horizontal axis, not the vertical one. For the tutor, this may be what's going wrong again. Then you add the distancefromaxis=8 option; that again would make most sense to me in case axis = horizontal, because then the used axis of rotation would be the horizontal line y=8, which would touch the curve you are defining. Otherwise you would rotate around x=8, sweeping the curve in a big circle to (x, y) = (16, 0) and around.

Finally, for your example with computation, I get the answer Pi*128/3 when computing by hand for rotation around x=6: rotating around x=6 means that we need to view everything as a function of y, which means that there are two parts: from x=0 to x=6-y (for y=4..6) and from x=6-y to x=6 (for y=0..4). Since we're rotating around a shifted axis, we need to correct for that and obtain:

Pi * (int((6 - 0)^2 - (6 - (6 - y))^2, y=4..6) + int((6 - (6 - y))^2 - (6 - 6)^2, y = 0..4))

which simplifies to

Pi * (int(36 - y^2, y=4..6) + int(y^2, y=0..4))

and both by hand and using Maple that gives Pi * 128/3.

Your result by hand and the accompanying formula I could only obtain with different parameters: rotating x = 6-y, taken between y=0..4 (or equivalently, x=2..6) around x=6.

I hope this may clear things up a little bit,

Erik Postma
Maplesoft.

Hi emoh,

I think this is a case where we could possibly do better with documentation, explaining what the intent is, but I think that it does perform as intended. You wrote:

About VolumeOfRevolution it says:

"The VolumeOfRevolution(f(x), x=a..b) command returns the volume of revolution of the expression f(x) from a to b."

Which is not true.

If you plot
VolumeOfRevolution(2*x^2,x=0..2,axis=vertical,'output'='plot), the program instead returns the volume of revolution of 2*x^2 from y=0..8.

However, if you leave out the axis option (axis = vertical) or replace it with the default, axis = horizontal, then I think you get what you would like to get. (If not, feel free to reply again and correct me.) The axis option is explained further down on that page, ?Student,Calculus1,VolumeOfRevolution , as follows:

axis = horizontal or vertical Whether the expression is rotated horizontally or vertically. By default, the rotation is horizontal.

The explanation here is not very good - what is meant is that the axis of rotation is horizontal or vertical, respectively, not the direction in which the plot is swept around. (As an aside, it's also of course not the expression that is rotated, but the curve.) So if you add the axis=vertical option, then the curve which runs from x=0 to x=2 (and from y=0 to y=8) is rotated around the y-axis to obtain the "outer surface" of the volume described. Indeed, if you plot(2*x^2, x=0..2) and mentally rotate the resulting curve around either the horizontal axis (y=0) or the vertical axis (x=0), you get what I get if I call VolumeOfRevolution(2*x^2, x=0..2, axis=horizontal, output=plot) or the same with axis = vertical, respectively.

For the second point you make, I would say the same thing is happening: what you describe indicates to me that you want to rotate around the horizontal axis, not the vertical one. For the tutor, this may be what's going wrong again. Then you add the distancefromaxis=8 option; that again would make most sense to me in case axis = horizontal, because then the used axis of rotation would be the horizontal line y=8, which would touch the curve you are defining. Otherwise you would rotate around x=8, sweeping the curve in a big circle to (x, y) = (16, 0) and around.

Finally, for your example with computation, I get the answer Pi*128/3 when computing by hand for rotation around x=6: rotating around x=6 means that we need to view everything as a function of y, which means that there are two parts: from x=0 to x=6-y (for y=4..6) and from x=6-y to x=6 (for y=0..4). Since we're rotating around a shifted axis, we need to correct for that and obtain:

Pi * (int((6 - 0)^2 - (6 - (6 - y))^2, y=4..6) + int((6 - (6 - y))^2 - (6 - 6)^2, y = 0..4))

which simplifies to

Pi * (int(36 - y^2, y=4..6) + int(y^2, y=0..4))

and both by hand and using Maple that gives Pi * 128/3.

Your result by hand and the accompanying formula I could only obtain with different parameters: rotating x = 6-y, taken between y=0..4 (or equivalently, x=2..6) around x=6.

I hope this may clear things up a little bit,

Erik Postma
Maplesoft.

Ah, now I think I understand - I wasn't really focusing on the volume at all, just on the surface. And the surface that was generated was, I think, the correct one, but it is not really made clear what bounds the actual volume.

Of course one can do almost exactly what you did but "say it differently", somewhat more in line with the original poster's question, as follows:

with(Student[Calculus1]):
Q1 := VolumeOfRevolution(0, sqrt(x), x=0..4, axis=vertical, output=plot): Q1;
Q2 := SurfaceOfRevolution(4, y=0..2, axis=horizontal, output=plot, surfaceoptions=[color=blue]): Q2;
plots[display](Q1, plottools[rotate](Q2, 0, Pi/2, 0), title="All together now...");

Best,

Erik Postma
Maplesoft.

Ah, now I think I understand - I wasn't really focusing on the volume at all, just on the surface. And the surface that was generated was, I think, the correct one, but it is not really made clear what bounds the actual volume.

Of course one can do almost exactly what you did but "say it differently", somewhat more in line with the original poster's question, as follows:

with(Student[Calculus1]):
Q1 := VolumeOfRevolution(0, sqrt(x), x=0..4, axis=vertical, output=plot): Q1;
Q2 := SurfaceOfRevolution(4, y=0..2, axis=horizontal, output=plot, surfaceoptions=[color=blue]): Q2;
plots[display](Q1, plottools[rotate](Q2, 0, Pi/2, 0), title="All together now...");

Best,

Erik Postma
Maplesoft.