FindRoot与Solve，NSolve和Reduce

如何解决FindRoot与Solve，NSolve和Reduce

| 首先是一些非必要的乐趣。我的真正问题远不止于此。请不要触摸拨盘。我正在使用Mathematica 8的新概率函数。目标是进行简单的功效分析。实验的功效为1减去II型错误的概率（即宣布“无效果”，而实际上是有效果的）。例如，我选择了一个实验来确定硬币是否公平。假设掷尾巴的概率由b（一个公平的硬币的b = 0.5）给出，那么确定该硬币对n次掷硬币的实验有偏见的能力由下式给出

1 - Probability[-in <= x - n/2 <= in,x \\[distributed] Binomialdistribution[n,b]]

我愿意称其为不怀疑的公平硬币与预期均值的偏差大小（选择该值是为了使公平硬币翻转n倍，则尾数将在均值范围内的时间约为95％ +/- in；这是BTW，它确定了I类错误的大小，即错误地声明存在效果的概率）。 Mathematica很好地绘制了计算出的功率图：

n = 40;
in = 6;
Plot[1-Probability[-in<=x-n/2<=in,b]],{b,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,b]] == 0.15 && 
  0 <= b <= 1,b]

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

这很惨，因为b = 0.75的功效是：

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

Out[54]= 0.896768

NSolve发现相同的结果。 Reduce执行以下操作：

In[55]:= res =  Reduce[Probability[-in <= x - n/2 <= in,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,b]] /. {ToRules[res]}

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

因此，ѭ5设法找到了两个解决方案，但发现了很多其他完全错误的解决方案。 FindRoot在这里效果最好：

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

Out[57]= {b -> 0.265122}

Out[58]= {b -> 0.734878}

好的，很长的介绍。我的问题是：为什么“解决”，“ NSolve”和“归约”在这里如此悲惨（且无声！）失败？恕我直言，这不是数值精度，因为为各种解决方案找到的功效值似乎是正确的（它们完美地位于功效曲线上），并且已从实际解决方案中大大移除。对于被mma8剥夺的Wizard先生：权力的表达很沉重：

In[42]:= Probability[-in <= x - n/2 <= in,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

我不会期望expected11来解决这个问题，但是我对ѭ4和ѭ5寄予了很高的期望。注意，对于n = 30，in = 5Solve，NSolve，Reduce和FindRoot都找到相同的正确解（当然，多项式阶次较低）。

解决方法

处理此数值方法时，不同的数值方法会有所不同。（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,Method -> \"Aberth\"];
r3 = NRoots[p[b] == .15,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,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.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.733751,0.834331,0.879148,\\
0.879148,0.910323,0.97317,1.08099,1.17529,\\
1.17529,1.23052,1.23052}

{-0.128504,0.75}

{-0.128504,0.734878,1.1285}

{-0.128504,1.1285}

看起来NSolve正在将NRoots与Aberth的方法一起使用，而Solve可能只是在调用NSolve。不同的解决方案集似乎遍布整个地图。实际上，许多声称是实数（但不是实数）的数字可能并不那么糟糕。我将比较一个这样的集合与通过对精确的根对象进行数字化处理而形成的集合（通常是安全的过程）的大小。

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

Out[77]= {0.128504,0.129867,0.13413,0.141881,\\
0.141881,0.154398,0.174443,0.209069,\\
0.265122,0.543986,0.575831,0.685011,\\
0.736383,0.889992,0.902725,\\
0.931526,1.19648,\\
1.24659,1.25157,1.44617,4.25448,4.25448}

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

Out[78]= {0.128504,0.543985,\\
0.734878,0.854255,0.94963,\\
1.01802,1.01802,1.06769,1.10183,1.12188,\\
1.12188,1.1285,4.25448}

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

Out[82]= {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}

丹尼尔·里奇布劳 ,我认为问题仅仅是找到高阶多项式的根的数字不稳定性：

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,MaxExtraConditions->0]
        1-p[b]/.%
Out[3]= {{b->0.75}}
Out[4]= {0.896768}

In[5]:= Solve[p[b]==0.15 && 0<=b<=1,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.}

In[7]:= Solve[p[b]==3/20 && 0<=b<=1,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,MaxExtraConditions->0]//N
        1-p[b]/.%
Out[9]= {{b->0.265122},{b->0.734878}}
Out[10]= {0.85,0.85}

（n.b.MaxExtraConditions->0实际上是默认选项，因此它可能不在上面。） Solve和ѭ5simply都只是生成Root对象当给出不精确的系数时，它们会自动进行数值评估。如果查看（缩短的）输出Out[7]，那么您将看到完整的40阶多项式的Root：

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}}。但是，要使平坦部件位于凸块的右侧，需要进行大量的数值抵消。这是不带显式ѭ28it选项的情况：该图清楚地说明了为什么Reduce（或Solve与MaxConditions->1，请参见上面的ѭ32properly）正确地（从左到右）找到第一个解决方案，并且找到第二个解决方案几乎正确，然后找到整个负载。 , 好吧，这不是一个正确的答案，而是一个有趣的发现。当使用magic（aka35ѭ）选项时，Solve[ ]与Reduce[ ]具有相同的行为：

n=40;
in=6;
Solve[Probability[-in<=x-n/2<=in,b]]==0.15 &&
      0<=b<=1,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}}

FindRoot与Solve，NSolve和Reduce

如何解决FindRoot与Solve，NSolve和Reduce

解决方法

相关推荐