FindRoot vs Solve, NSolve and Reduce

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

开发者_如何学Python

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).

---

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.

---

Different numeric methods will fare differently when handling this.

(1) The ones that find all polynomial roots have the most difficult job, in that they may need to deal with deflated polynomials. FindRoot is off the hook there.

(2) The polynomial is a perturbation of one with substantial multiplicity. I would expect numeric methods to have trouble.

(3) The roots are all within 1-2 orders of magnitude in size. SO this is not so far from generally "bad" polynomials with roots around the unit circle.

(4) Most difficult is handling Solve[numeric eqn and ineq]. This must combine inequality solving methods (i.e. cylindrical decomposition) with machine arithmetic. Expect little mercy. Okay, this is univariate, so it amounts to Sturm sequences or Descartes' Rule of Signs. Still not numerically well behaved.

Here are some experiments using various method settings.

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

r1 = NRoots[p[b] == .15, b, Method -> "JenkinsTraub"];
r2 = NRoots[p[b] == .15, b, Method -> "Aberth"];
r3 = NRoots[p[b] == .15, b, Method -> "CompanionMatrix"];
r4 = NSolve[p[b] == .15, b];
r5 = Solve[p[b] == 0.15, b];
r6 = Solve[p[b] == 0.15 && Element[b, Reals], b];
r7 = N[Solve[p[b] == 15/100 && Element[b, Reals], b]]; 
r8 = N[Solve[p[b] == 15/100, b]];

Sort[Cases[b /. {ToRules[r1]}, _Real]]
Sort[Cases[b /. {ToRules[r2]}, _Real]]
Sort[Cases[b /. {ToRules[r3]}, _Real]]
Sort[Cases[b /. r4, _Real]]
Sort[Cases[b /. r5, _Real]]
Sort[Cases[b /. r6, _Real]]
Sort[Cases[b /. r7, _Real]]
Sort[Cases[b /. r8, _Real]]

{-0.128504, 0.265122, 0.728, 1.1807, 1.20794, 1.22063}

{-0.128504, 0.265122, 0.736383, 0.801116, 0.825711, 0.845658, \
0.889992, 0.931526, 0.958879, 0.986398, 1.06506, 1.08208, 1.18361, \
1.19648, 1.24659, 1.25157}

{-0.128504, 0.265122, 0.733751, 0.834331, 0.834331, 0.879148, \
0.879148, 0.910323, 0.97317, 0.97317, 1.08099, 1.08099, 1.17529, \
1.17529, 1.23052, 1.23052}

{-0.128504, 0.265122, 0.736383, 0.801116, 0.825711, 0.845658, \
0.889992, 0.931526, 0.958879, 0.986398, 1.06506, 1.08208, 1.18361, \
1.19648, 1.24659, 1.25157}

{-0.128504, 0.265122, 0.736383, 0.801116, 0.825711, 0.845658, \
0.889992, 0.931526, 0.958879, 0.986398, 1.06506, 1.08208, 1.18361, \
1.19648, 1.24659, 1.25157}

{-0.128504, 0.75}

{-0.128504, 0.265122, 0.734878, 1.1285}

{-0.128504, 0.265122, 0.734878, 1.1285}

It looks like NSolve is using NRoots with Aberth's method, and Solve might just be calling NSolve.

The distinct solution sets seem to be all over the map. Actually many of the numeric ones that claim to be real (but aren't) might not be so bad. I'll compare magnitudes of one such set vs a set formed from numericizing exact root objects (a generally safe process).

mags4 = Sort[Abs[b /. r4]]

Out[77]= {0.128504, 0.129867, 0.129867, 0.13413, 0.13413, 0.141881, \
0.141881, 0.154398, 0.154398, 0.174443, 0.174443, 0.209069, 0.209069, \
0.265122, 0.543986, 0.543986, 0.575831, 0.575831, 0.685011, 0.685011, \
0.736383, 0.801116, 0.825711, 0.845658, 0.889992, 0.902725, 0.902725, \
0.931526, 0.958879, 0.986398, 1.06506, 1.08208, 1.18361, 1.19648, \
1.24659, 1.25157, 1.44617, 1.44617, 4.25448, 4.25448}

mags8 = Sort[Abs[b /. r8]]

Out[78]= {0.128504, 0.129867, 0.129867, 0.13413, 0.13413, 0.141881, \
0.141881, 0.154398, 0.154398, 0.174443, 0.174443, 0.209069, 0.209069, \
0.265122, 0.543985, 0.543985, 0.575831, 0.575831, 0.685011, 0.685011, \
0.734878, 0.854255, 0.854255, 0.902725, 0.902725, 0.94963, 0.94963, \
1.01802, 1.01802, 1.06769, 1.06769, 1.10183, 1.10183, 1.12188, \
1.12188, 1.1285, 1.44617, 1.44617, 4.25448, 4.25448}

Chop[mags4 - mags8, 10^(-6)]

Out[82]= {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \
0.00150522, -0.0531384, -0.0285437, -0.0570674, -0.0127339, \
-0.0469044, -0.0469044, -0.0864986, -0.0591449, -0.0812974, \
-0.00263812, -0.0197501, 0.0817724, 0.0745959, 0.124706, 0.123065, 0, \
0, 0, 0}

Daniel Lichtblau

---

Well, not a proper answer, but an interesting observation. Solve[ ] has the same behavior than Reduce[ ] when the magic (aka MaxExtraConditions) option is used:

n=40;
in=6;
Solve[Probability[-in<=x-n/2<=in,
      x\[Distributed]BinomialDistribution[n,b]]==0.15 &&
      0<=b<=1,b, MaxExtraConditions->1]

{{b -> 0.265122}, {b -> 0.736488}, {b -> 0.80151}, {b -> 0.825884}, 
 {b -> 0.84573}, {b -> 0.890444}, {b -> 0.931972}, {b -> 0.960252}, 
 {b -> 0.985554}}