O gerador de números aleatórios do Mathematica se desvia da probabilidade binomial?

9

Então, digamos que você jogue uma moeda 10 vezes e chame esse 1 de "evento". Se você executar 1.000.000 desses "eventos", qual a proporção de eventos com cabeças entre 0,4 e 0,6? A probabilidade binomial sugere que isso seja cerca de 0,65, mas meu código do Mathematica está me dizendo cerca de 0,24

Aqui está a minha sintaxe:

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

Onde está o acidente?

Tim McKnight
fonte
3
talvez isso seria mais adequado para mathematica Stackexchange mathematica.stackexchange.com
Jeromy Anglim
11
@ JeremyAnglim Neste caso, suspeito que o problema provavelmente esteja no raciocínio, e não estritamente na codificação.
Glen_b -Reinstala Monica
@ Glen_b Eu acho que o principal é que existe uma boa resposta em algum lugar da internet, que você parece ter fornecido. :-)
Jeromy Anglim

Respostas:

19

O contratempo é o uso de menos de estrito.

Com dez lançamentos, a única maneira de obter um resultado de proporção de cabeças estritamente entre 0,4 e 0,6 é se você obtiver exatamente 5 cabeças. Isso tem uma probabilidade de cerca de 0,246 ( ), que é sobre o que suas simulações (corretamente ) dar.(105)(12)100.246

Se você incluir 0,4 e 0,6 em seus limites (ou seja, 4, 5 ou 6 cabeças em 10 lançamentos), o resultado terá uma probabilidade de cerca de 0,656, conforme o esperado.

Seu primeiro pensamento não deve ser um problema com o gerador de números aleatórios. Esse tipo de problema teria sido óbvio em um pacote muito usado como o Mathematica, muito antes de agora.

Glen_b -Reinstate Monica
fonte
Ironicamente, o @TimMcKnight demonstrou probabilidade binomial para nós.
Simon Kuang 02/02
8

Alguns comentários sobre o código que você escreveu:

  • Você definiu, experiment[n_]mas nunca o usou, em vez disso, repete sua definição em trialheadcount[n_].
  • experiment[n_]poderia ser programado com muito mais eficiência (sem usar o comando interno BinomialDistribution), pois Total[RandomInteger[{0,1},n]/nisso também tornaria Xdesnecessário.
  • Contar o número de casos em que experiment[n_]é estritamente entre 0,4 e 0,6 é realizado com mais eficiência por escrito Length[Select[Table[experiment[10],{10^6}], 0.4 < # < 0.6 &]].

Mas, para a pergunta propriamente dita, como Glen_b aponta, a distribuição binomial é discreta. De 10 lançamentos de moedas com cabeças observadas, a probabilidade de que a proporção da amostra de cabeças seja estritamente entre 0,4 e 0,6 é realmente apenas o caso ; ou seja, Enquanto que, se você calcular a probabilidade de que a proporção da amostra esteja entre 0,4 e 0,6 inclusive , isso seria Portanto, você só precisa modificar seu código para usarxp^=x/10x=5Pr[4X6]=6 x=4 ( 10

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.6em vez de. Mas é claro, também poderíamos escrever
Length[Select[RandomVariate[BinomialDistribution[10,1/2],{10^6}], 4 <= # <= 6 &]]

Este comando é aproximadamente 9,6 vezes mais rápido que o seu código original. Imagino que alguém ainda mais competente do que eu sou no Mathematica possa acelerar ainda mais.

heropup
fonte
2
Você pode acelerar o seu código por outro fator de 10 usando Total@Map[Counts@RandomVariate[BinomialDistribution[10, 1/2], 10^6], {4, 5, 6}]. Eu suspeito que Counts[], sendo uma função interna, seja altamente otimizada em comparação com Select[], que tem que trabalhar com predicados arbitrários.
David Zhang
1

Fazendo experimentos de probabilidade no Mathematica

O Mathematica oferece uma estrutura muito confortável para trabalhar com probabilidades e distribuições e - embora a questão principal dos limites apropriados tenha sido abordada - eu gostaria de usar esta pergunta para tornar isso mais claro e talvez útil como referência.

Vamos simplesmente tornar os experimentos repetíveis e definir algumas opções de plotagem que se ajustem ao nosso gosto:

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

Trabalhando com Distribuições Paramétricas

Agora podemos definir a distribuição assintótica de um evento que é a proporção de cabeças em arremessos de uma moeda (justa):nπ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
        ]
    }
]

O que nos dá o gráfico da distribuição discreta de proporções: TheoreticalDistributionPlot

Podemos usar a distribuição imediatamente para calcular probabilidades para e :Prr[Pr[0.4π0.6|πB(10,12)]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}

Fazendo experiências em Monte Carlo

Podemos usar a distribuição para um evento para amostrá-lo 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

Comparando isso com a distribuição teórica / assintótica, mostra que tudo se encaixa:

Show @ {
   theoreticalPlot,
   empiricalPlot
}

ComparingDistributions

gwr
fonte
Você pode encontrar uma postagem semelhante com mais informações sobre o Mathematica no Mathematica SE .
gwr