FindRoot vs Solve,NSolve和Reduce

Question

首先是一些有趣的非必要背景.我真正的问题远远低于.请勿触摸表盘.

我正在玩Mathematica 8的新概率函数.目标是进行简单的功率分析.实验的力量是1减去II型错误的概率(即,反弹'没有效果',而实际上有效果).

作为一个例子,我选择了一个实验来确定硬币是否公平.假设抛出尾巴的概率由b给出(一个公平的硬币具有b = 0.5),那么确定硬币偏向于用n个硬币翻转的实验的能力由下式给出:

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

与在与预期均值的偏差一个公平的硬币,我一个愿意叫不怀疑(的大小在选择,这样对于一个公平的硬币翻转ň次尾的数量将是平均内的时间约95% +/- in ;这个,BTW,确定I型错误的大小,错误地声称存在效果的概率).

Mathematica很好地绘制了计算能力的图:

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]

我以85%的功率绘制了一条线,这通常被认为是合理的功率量.现在,我想要的只是功率曲线与此线相交的点.这告诉我硬币必须具有的最小偏差,以便我有合理的期望在40次翻转的实验中找到它.

所以,我尝试过:

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}}

这失败了,因为对于b = 0.75,功率是:

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

Out[54]= 0.896768

NSolve找到相同的结果.Reduce执行以下操作:

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

因此,Reduce设法找到这两个解决方案,但它找到了其他一些错误的解决方案.

FindRoot 在这里效果最好:

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}

好的,长期介绍.我的问题是:为什么Solve,NSolve和Reduce失败如此悲惨(并且默默地!)在这里？恕我直言,它不能是数字精度,因为各种解决方案的功率值似乎是正确的(它们完全位于功率曲线上),并且从实际解决方案中大大消除.

对于被剥夺了mma8的Mr.Wizard:权力的表达很重要:

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

我也没有想到Solve来处理这个问题,但我寄予厚望NSolve和Reduce.请注意,ñ = 30,在 = 5 Solve,NSolve,Reduce和FindRoot所有发现同样的,正确的解决方案(当然,多项式阶较低那里).

Answer 1

我认为问题只是找到高阶多项式的根的数值不稳定性:

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}

(nb MaxExtraConditions->0实际上是默认选项,因此它可能被排除在上面.)

二者Solve并Reduce简单地产生Root对象和给定的不精确系数时,它们会自动进行数值计算.如果你查看(缩短的)输出,Out[7]那么你将看到Root完整的40阶多项式:

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]

在此图表中,您可以确认零点位于(约){{b - > 0.265122},{b - > 0.734878}}.但是,要在凸起的右侧获得平坦部件需要大量的数字取消.这是没有显式WorkingPrecision选项的情况:

该曲线图清楚为什么Reduce(或Solve用MaxConditions->1,见In[5]上文)发现(从左至右)适当地将第一溶液和几乎正确的第二溶液中,然后污物的整个负载.

Answer 2

处理此问题时,不同的数字方法会有所不同.

(1)找到所有多项式根的那些具有最困难的工作,因为它们可能需要处理缩小的多项式.FindRoot已经摆​​脱困境.

(2)多项式是具有实质多重性的多项式的扰动.我希望数字方法有麻烦.

(3)根的大小都在1-2个数量级之内.因此,这与通常在单位圆周围有根的"坏"多项式相差不大.

(4)最困难的是处理Solve [数字eqn和ineq].这必须将不等式求解方法(即圆柱分解)与机器算法结合起来.期待一点怜悯.好吧,这是单变量的,所以它相当于Sturm序列或笛卡尔的符号规则.仍然没有数字表现良好.

以下是使用各种方法设置的一些实验.

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}

看起来NSolve正在使用带有Aberth方法的NRoots,而Solve可能只是在调用NSolve.

不同的解决方案集似乎遍布地图.实际上许多声称是真实的(但不是)的数字可能并不那么糟糕.我将比较一个这样的集合的大小与通过数字化精确根对象形成的集合(通常是安全的过程).

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