Primero un contexto no esencial para la diversión. Mi verdadera pregunta está muy abajo. Por favor, no toque el dial.FindRoot vs Solve, NSolve and Reduce
Estoy jugando con las nuevas funciones probabilísticas de Mathematica 8. El objetivo es hacer un análisis de potencia simple. El poder de un experimento es 1 menos la probabilidad de un error tipo II (es decir, anunciar 'no efecto', mientras que hay un efecto en la realidad).
Como ejemplo elegí un experimento para determinar si una moneda es justa. Supongamos que la probabilidad de tirar colas está dada por b (una moneda tiene b = 0,5), entonces la potencia para determinar que la moneda es sesgada para un experimento con n moneda voltea está dada por
1 - Probability[-in <= x - n/2 <= in, x \[Distributed] BinomialDistribution[n, b]]
con en el tamaño de la desviación de la media esperada para una moneda que un dispuestos a llamar a sospechoso (en se elige de forma que para una moneda volcó n veces el número de colas será de unos 95% del tiempo dentro de la media +/- en; esto, por cierto, de termina el tamaño del error tipo I, la probabilidad de reclamar incorrectamente la existencia de un efecto).
Mathematica bien dibuja un gráfico de la potencia calculada:
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 dibujó una línea a una potencia de 85%, lo que generalmente se considera que es una cantidad razonable de poder. Ahora, todo lo que quiero es los puntos donde la curva de potencia se cruza con esta línea. Esto me dice el sesgo mínimo que la moneda debe tener para que tenga una expectativa razonable de encontrarlo en un experimento con 40 volteos.
tanto, he intentado:
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}}
Esta falla estrepitosamente, ya que para b = 0,75 la potencia es:
In[54]:= 1 - Probability[-in <= x - n/2 <= in, x \[Distributed] BinomialDistribution[n, 0.75]]
Out[54]= 0.896768
NSolve encuentra el mismo resultado. Reduce realiza lo siguiente:
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.}
Así, Reduce las arregla para encontrar las dos soluciones, sino que encuentra un buen número de otros que son totalmente equivocado.
FindRoot funciona mejor aquí:
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}
bien, larga introducción. Mi pregunta es: ¿por qué Solve, NSolve y Reduce fallan tan miserablemente (y silenciosamente) aquí? En mi humilde opinión, no puede ser precisión numérica ya que los valores de potencia encontrados para las diversas soluciones parecen ser correctos (se encuentran perfectamente en la curva de potencia) y se eliminan considerablemente de la solución real.
Para el Sr. mma8 privado.Asistente: La expresión de la potencia es uno pesado:
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
y no habría esperado Solve para manejar esto, pero tenía grandes esperanzas para NSolve y Reduce. Tenga en cuenta que para n = 30, en = 5 Solve, NSolve, Reduce y FindRoot todos buscan el mismo, las soluciones correctas (por supuesto, el orden polinomio es menor allí).
me dejó un comentario en virtud de la respuesta de Daniel. Por lo que vale, Simplificar [desordenado && 0 <= b <= 1] es una ruta rápida hacia las soluciones. – telefunkenvf14