¿El generador de números aleatorios de Mathematica se desvía de la probabilidad binomial?

9

Entonces, digamos que lanzas una moneda 10 veces y llamas a ese 1 "evento". Si ejecuta, 1,000,000 de estos "eventos", ¿cuál es la proporción de eventos que tienen cabezas entre 0.4 y 0.6? La probabilidad binomial sugeriría que esto sea aproximadamente 0.65, pero mi código de Mathematica me dice aproximadamente 0.24

Aquí está mi sintaxis:

In[2]:= X:= RandomInteger[];
In[3]:= experiment[n_]:= Apply[Plus, Table[X, {n}]]/n;
In[4]:= trialheadcount[n_]:= .4 < Apply[Plus, Table[X, {n}]]/n < .6
In[5]:= sample=Table[trialheadcount[10], {1000000}]
In[6]:= Count[sample2,True];
Out[6]:= 245682

¿Dónde está el contratiempo?

Tim McKnight
fuente
3
tal vez esto sería más adecuado a Mathematica StackExchange mathematica.stackexchange.com
Jeromy Anglim
1
@JeromyAnglim En este caso, sospecho que el problema probablemente se deba al razonamiento y no estrictamente a la codificación.
Glen_b -Reinstala a Monica
@Glen_b Supongo que lo principal es que hay una buena respuesta en algún lugar de Internet, que parece haber proporcionado. :-)
Jeromy Anglim

Respuestas:

19

El contratiempo es el uso de estricto menor que.

Con diez lanzamientos, la única forma de obtener un resultado de proporción de cabezas estrictamente entre 0.4 y 0.6 es si obtiene exactamente 5 caras. Eso tiene una probabilidad de aproximadamente 0.246 ( ), que es sobre cuáles son sus simulaciones (correctamente ) dar.(105)(12)100.246

Si incluye 0.4 y 0.6 en sus límites (es decir, 4, 5 o 6 caras en 10 lanzamientos), el resultado tiene una probabilidad de aproximadamente 0.656, más de lo que esperaba.

Su primer pensamiento no debería ser un problema con el generador de números aleatorios. Ese tipo de problema habría sido obvio en un paquete muy usado como Mathematica mucho antes.

Glen_b -Reinstate a Monica
fuente
Irónicamente, @TimMcKnight demostró la probabilidad binomial para nosotros.
Simon Kuang
8

Algunos comentarios sobre el código que escribiste:

  • Lo definió experiment[n_]pero nunca lo usó, sino que repitió su definición en trialheadcount[n_].
  • experiment[n_]podría programarse de manera mucho más eficiente (sin usar el comando incorporado BinomialDistribution) como Total[RandomInteger[{0,1},n]/ny esto también haría Xinnecesario.
  • Contar el número de casos donde experiment[n_]está estrictamente entre 0.4 y 0.6 se logra más eficientemente escribiendo Length[Select[Table[experiment[10],{10^6}], 0.4 < # < 0.6 &]].

Pero, para la pregunta en sí misma, como señala Glen_b, la distribución binomial es discreta. De cada 10 lanzamientos de monedas con caras observadas, la probabilidad de que la proporción muestral de caras esté estrictamente entre 0.4 y 0.6 es en realidad el caso ; es decir, Mientras que, si calculara la probabilidad de que la proporción de la muestra esté entre 0.4 y 0.6 inclusive , sería Por lo tanto, solo necesita modificar su código para usarxp^=x/10x=5

Pr[X=5]=(105)(0.5)5(10.5)50.246094.
Pr[4X6]=x=46(10x)(0.5)x(10.5)10x=67210240.65625.
0.4 <= # <= 0.6en lugar. Pero claro, también podríamos escribir
Length[Select[RandomVariate[BinomialDistribution[10,1/2],{10^6}], 4 <= # <= 6 &]]

Este comando es aproximadamente 9.6 veces más rápido que su código original. Me imagino que alguien aún más competente que yo en Mathematica podría acelerarlo aún más.

heropup
fuente
2
Puede acelerar su código por otro factor de 10 usando Total@Map[Counts@RandomVariate[BinomialDistribution[10, 1/2], 10^6], {4, 5, 6}]. Sospecho que Counts[], al ser una función incorporada, está altamente optimizado en comparación con Select[], que tiene que trabajar con predicados arbitrarios.
David Zhang
1

Haciendo experimentos de probabilidad en Mathematica

Mathematica ofrece un marco muy cómodo para trabajar con probabilidades y distribuciones y, si bien se ha abordado el problema principal de los límites apropiados, me gustaría utilizar esta pregunta para aclarar esto y tal vez sea útil como referencia.

Simplemente hagamos que los experimentos sean repetibles y definamos algunas opciones de trama que se ajusten a nuestro gusto:

SeedRandom["Repeatable_151115"];
$PlotTheme = "Detailed";
SetOptions[Plot, Filling -> Axis];
SetOptions[DiscretePlot, ExtentSize -> Scaled[0.5], PlotMarkers -> "Point"];

Trabajando con distribuciones paramétricas

Ahora podemos definir la distribución asintótica para un evento que es la proporción de caras en tiros de una moneda (regular):πn

distProportionTenCoinThrows = With[
    {
        n = 10, (* number of coin throws *)
        p = 1/2 (* fair coin probability of head*)
    },
    (* derive the distribution for the proportion of heads *)
    TransformedDistribution[
        x/n,
        x \[Distributed] BinomialDistribution[ n, p ]
    ];

With[
    {
        pr = PlotRange -> {{0, 1}, {0, 0.25}}
    },
    theoreticalPlot = DiscretePlot[
        Evaluate @ PDF[ distProportionTenCoinThrows, p ],
        {p, 0, 1, 0.1},
        pr
    ];
    (* show plot with colored range *)
    Show @ {
        theoreticalPlot,
        DiscretePlot[
            Evaluate @ PDF[ distProportionTenCoinThrows, p ],
            {p, 0.4, 0.6, 0.1},
            pr,
            FillingStyle -> Red,
            PlotLegends -> None
        ]
    }
]

Lo que nos da la gráfica de la distribución discreta de proporciones: Distribución teórica

Podemos usar la distribución inmediatamente para calcular las probabilidades de Pr[0.4π0.6|πB(10,12)] y :Pr[0.4<π<0.6|πB(10,12)]

{
    Probability[ 0.4 <= p <= 0.6, p \[Distributed] distProportionTenCoinThrows ],
    Probability[ 0.4 < p < 0.6, p \[Distributed] distProportionTenCoinThrows ]
} // N

{0.65625, 0.246094}

Haciendo experimentos de Monte Carlo

Podemos usar la distribución de un evento para muestrearlo repetidamente (Monte Carlo).

distProportionsOneMillionCoinThrows = With[
    {
        sampleSize = 1000000
    },
    EmpiricalDistribution[
        RandomVariate[
            distProportionTenCoinThrows,
            sampleSize
        ]
    ]
];

empiricalPlot = 
    DiscretePlot[
        Evaluate@PDF[ distProportionsOneMillionCoinThrows, p ],
        {p, 0, 1, 0.1}, 
        PlotRange -> {{0, 1}, {0, 0.25}} , 
        ExtentSize -> None, 
        PlotLegends -> None, 
        PlotStyle -> Red
    ]
]

EmpirialDistributionPlot

Comparar esto con la distribución teórica / asintótica muestra que todo encaja:

Show @ {
   theoreticalPlot,
   empiricalPlot
}

Comparación de distribuciones

gwr
fuente
Puede encontrar un puesto similar con más información de fondo con respecto a Mathematica en Mathematica SE .
gwr