FindRoot vs Solve, NSolve and Reduce

Question

First some non-essential context for fun. My real question is far below. Please don't touch the dial.

I'm playing with the new probabilistic functions of Mathematica 8. Goal is to do a simple power analysis. The power of an experiment is 1 minus the probability of a type II error (i.e., anouncing 'no effect', whereas there is an effect in reality).

As an example I chose an experiment to determine whether a coin is fair. Suppose the probability to throw tails is given by b (a fair coin has b=0.5), then the power to determine that the coin is biased for an experiment with n coin flips is given by

1 - Probability[-in <= x - n/2 <= in, x \[Distributed] BinomialDistribution[n, b]]

with in the size of the deviation from the expected mean for a fair coin that I an willing to call not suspicious (in is chosen so that for a fair coin flipped n times the number of tails will be about 95% of the time within mean +/- in ; this, BTW, determines the size of the type I error, the probability to incorrectly claim the existence of an effect).

Mathematica nicely draws a plot of the calculated power:

n = 40;
in = 6;
Plot[1-Probability[-in<=x-n/2<=in,x \[Distributed] BinomialDistribution[n, b]], {b, 0, 1},
 Epilog -> Line[{{0, 0.85}, {1, 0.85}}], Frame -> True,
 FrameLabel -> {"P(tail)", "Power", "", ""},
 BaseStyle -> {FontFamily -> "Arial", FontSize -> 16, 
   FontWeight -> Bold}, ImageSize -> 500]

I drew a line at a power of 85%, which is generally considered to be a reasonable amount of power. Now, all I want is the points where the power curve intersects with this line. This tells me the minimum bias the coin must have so that I have a reasonable expectation to find it in an experiment with 40 flips.

So, I tried:

In[47]:= Solve[ Probability[-in <= x - n/2 <= in, 
    x \[Distributed] BinomialDistribution[n, b]] == 0.15 && 
  0 <= b <= 1, b]

Out[47]= {{b -> 0.75}}

This fails miserably, because for b = 0.75 the power is:

In[54]:= 1 - Probability[-in <= x - n/2 <= in, x \[Distributed] BinomialDistribution[n, 0.75]]

Out[54]= 0.896768

NSolve finds the same result. Reducedoes the following:

In[55]:= res =  Reduce[Probability[-in <= x - n/2 <= in, 
     x \[Distributed] BinomialDistribution[n, b]] == 0.15 && 
   0 <= b <= 1, b, Reals]

Out[55]= b == 0.265122 || b == 0.73635 || b == 0.801548 || 
 b == 0.825269 || b == 0.844398 || b == 0.894066 || b == 0.932018 || 
 b == 0.957616 || b == 0.987099

In[56]:= 1 -Probability[-in <= x - n/2 <= in, 
              x \[Distributed] BinomialDistribution[n, b]] /. {ToRules[res]}

Out[56]= {0.85, 0.855032, 0.981807, 0.994014, 0.99799, 0.999965, 1., 1., 1.}

So, Reduce manages to find the two solutions, but it finds quite a few others that are dead wrong.

FindRoot works best here:

In[57]:= FindRoot[{Probability[-in <= x - n/2 <= in, 
             x \[Distributed] BinomialDistribution[n, b]] - 0.15`}, {b, 0.2, 0, 0.5}]
         FindRoot[{Probability[-in <= x - n/2 <= in, 
             x \[Distributed] BinomialDistribution[n, b]] - 0.15`}, {b, 0.8, 0.5, 1}]

Out[57]= {b -> 0.265122}

Out[58]= {b -> 0.734878}

OK, long introduction. My question is: why do Solve, NSolve, and Reduce fail so miserably (and silently!) here? IMHO, it can't be numerical accuracy since the power values found for the various solutions seem to be correct (they lie perfectly on the power curve) and are considerably removed from the real solution.

For the mma8-deprived Mr.Wizard: The expression for the power is a heavy one:

In[42]:= Probability[-in <= x - n/2 <= in, 
 x \[Distributed] BinomialDistribution[n, b]]

Out[42]= 23206929840 (1 - b)^26 b^14 + 40225345056 (1 - b)^25 b^15 + 
 62852101650 (1 - b)^24 b^16 + 88732378800 (1 - b)^23 b^17 + 
 113380261800 (1 - b)^22 b^18 + 131282408400 (1 - b)^21 b^19 + 
 137846528820 (1 - b)^20 b^20 + 131282408400 (1 - b)^19 b^21 + 
 113380261800 (1 - b)^18 b^22 + 88732378800 (1 - b)^17 b^23 + 
 62852101650 (1 - b)^16 b^24 + 40225345056 (1 - b)^15 b^25 + 
 23206929840 (1 - b)^14 b^26

and I wouldn't have expected Solve to handle this, but I had high hopes for NSolve and Reduce. Note that for n=30, in=5 Solve, NSolve, Reduce and FindRoot all find the same, correct solutions (of course, the polynomial order is lower there).

Answer 1

I think the problem is just the numeric instablitity of finding roots to high order polynomials:

In[1]:= n=40; in=6;
        p[b_]:= Probability[-in<=x-n/2<=in,
                            x\[Distributed]BinomialDistribution[n,b]]

In[3]:= Solve[p[b]==0.15 && 0<=b<=1, b, MaxExtraConditions->0]
        1-p[b]/.%
Out[3]= {{b->0.75}}
Out[4]= {0.896768}

In[5]:= Solve[p[b]==0.15 && 0<=b<=1, b, MaxExtraConditions->1]
        1-p[b]/.%
Out[5]= {{b->0.265122},{b->0.736383},{b->0.801116},{b->0.825711},{b->0.845658},{b->0.889992},{b->0.931526},{b->0.958879},{b->0.986398}}
Out[6]= {0.85,0.855143,0.981474,0.994151,0.998143,0.999946,1.,1.,1.}

In[7]:= Solve[p[b]==3/20 && 0<=b<=1, b, MaxExtraConditions->0]//Short
        1-p[b]/.%//N
Out[7]//Short= {{b->Root[-1+<<39>>+108299005920 #1^40&,2]},{b->Root[<<1>>&,3]}}
Out[8]= {0.85,0.85}

In[9]:= Solve[p[b]==0.15`100 && 0<=b<=1, b, MaxExtraConditions->0]//N
        1-p[b]/.%
Out[9]= {{b->0.265122},{b->0.734878}}
Out[10]= {0.85,0.85}

(n.b. MaxExtraConditions->0 is actually the default option, so it could have been left out of the above.)

Both Solve and Reduce are simply generating Root objects and when given inexact coefficients, they are automatically numerically evaluated. If you look at the (shortened) output Out[7] then you'll see the Root of the full 40th order polynomial:

In[12]:= Expand@(20/3 p[b] - 1)
Out[12]= -1 + 154712865600 b^14 - 3754365538560 b^15 + 43996471155000 b^16 - 
         331267547520000 b^17 + 1798966820560000 b^18 - 
         7498851167808000 b^19 + 24933680132961600 b^20 - 
         67846748661120000 b^21 + 153811663157880000 b^22 - 
         294248399084640000 b^23 + 479379683508726000 b^24 - 
         669388358063093760 b^25 + 804553314979680000 b^26 - 
         834351666126339200 b^27 + 747086226686186400 b^28 - 
         577064755104364800 b^29 + 383524395817442880 b^30 - 
         218363285636496000 b^31 + 105832631433929400 b^32 - 
         43287834659596800 b^33 + 14776188957129600 b^34 - 
         4150451102878080 b^35 + 942502182076000 b^36 - 
         168946449235200 b^37 + 22970789150400 b^38 -
         2165980118400 b^39 + 108299005920 b^40
In[13]:= Plot[%, {b, -1/10, 11/10}, WorkingPrecision -> 100]

From this graph you can confirm that the zeros are at (approx) {{b -> 0.265122}, {b -> 0.734878}}. But, to get the flat parts on the right hand side of the bump requires lots of numerical cancellations. Here's what it looks like without the explicit WorkingPrecision option:

This graph makes it clear why Reduce (or Solve with MaxConditions->1, see In[5] above) finds (from left to right) the first solution properly and the second solution almost correctly, followed by a whole load of crud.