Mathematikas Zufallszahlengenerator, der von der Binomialwahrscheinlichkeit abweicht?

9

Nehmen wir also an, Sie werfen 10 Mal eine Münze und nennen diese 1 "Ereignis". Wenn Sie 1.000.000 dieser "Ereignisse" ausführen, wie hoch ist der Anteil der Ereignisse mit Köpfen zwischen 0,4 und 0,6? Die Binomialwahrscheinlichkeit würde bedeuten, dass dies ungefähr 0,65 ist, aber mein Mathematica-Code sagt mir ungefähr 0,24

Hier ist meine Syntax:

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

Wo ist das Missgeschick?

Tim McKnight
quelle
3
Vielleicht wäre dies besser geeignet für mathematica stackexchange mathematica.stackexchange.com
Jeromy Anglim
1
@ JeromyAnglim In diesem Fall vermute ich, dass das Problem wahrscheinlich eher in der Begründung als in der Kodierung liegt.
Glen_b -State Monica
@Glen_b Ich denke, die Hauptsache ist, dass es irgendwo im Internet eine gute Antwort gibt, die Sie anscheinend gegeben haben. :-)
Jeromy Anglim

Antworten:

19

Das Missgeschick ist die Verwendung von strengen weniger als.

Mit zehn Würfen ist die einzige Möglichkeit, ein Ergebnis mit einem Kopfanteil von genau zwischen 0,4 und 0,6 zu erzielen, wenn Sie genau 5 Köpfe erhalten. Das hat eine Wahrscheinlichkeit von ungefähr 0,246 ( ), was ungefähr Ihren Simulationen entspricht (richtig ) geben.(105)(12)100.246

Wenn Sie 0,4 und 0,6 in Ihre Limits aufnehmen (dh 4, 5 oder 6 Köpfe in 10 Würfen), hat das Ergebnis eine Wahrscheinlichkeit von ungefähr 0,656, so wie Sie es erwartet haben.

Ihr erster Gedanke sollte kein Problem mit dem Zufallszahlengenerator sein. Diese Art von Problem wäre in einem stark genutzten Paket wie Mathematica schon lange vorher offensichtlich gewesen.

Glen_b - Monica neu starten
quelle
Ironischerweise hat @TimMcKnight für uns eine Binomialwahrscheinlichkeit gezeigt.
Simon Kuang
8

Einige Kommentare zu dem Code, den Sie geschrieben haben:

  • Sie haben es definiert, experiment[n_]aber nie verwendet, sondern seine Definition in wiederholt trialheadcount[n_].
  • experiment[n_]könnte viel effizienter programmiert werden (ohne den eingebauten Befehl zu verwenden BinomialDistribution) als Total[RandomInteger[{0,1},n]/nund dies würde auch Xunnötig machen .
  • Das Zählen der Anzahl von Fällen, in denen experiment[n_]streng zwischen 0,4 und 0,6 liegt, wird durch Schreiben effizienter erreicht Length[Select[Table[experiment[10],{10^6}], 0.4 < # < 0.6 &]].

Für die eigentliche Frage selbst ist die Binomialverteilung jedoch diskret, wie Glen_b hervorhebt. Von 10 Münzwürfen mit beobachtet Köpfen, die Wahrscheinlichkeit , dass die Probe Anteil der Köpfe ist streng zwischen 0,4 und 0,6 ist eigentlich nur der Fall ; dh Wenn Sie dagegen die Wahrscheinlichkeit berechnen würden, dass der Stichprobenanteil zwischen 0,4 und 0,6 einschließlich liegt, wäre dies Daher müssen Sie nur Ihren Code ändern, um ihn zu verwendenxp^=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.6stattdessen. Aber wir könnten natürlich auch schreiben
Length[Select[RandomVariate[BinomialDistribution[10,1/2],{10^6}], 4 <= # <= 6 &]]

Dieser Befehl ist ungefähr 9,6-mal schneller als Ihr ursprünglicher Code. Ich stelle mir vor, jemand, der noch kompetenter ist als ich bei Mathematica, könnte es noch weiter beschleunigen.

Heropup
quelle
2
Sie können Ihren Code mithilfe von um einen weiteren Faktor 10 beschleunigen Total@Map[Counts@RandomVariate[BinomialDistribution[10, 1/2], 10^6], {4, 5, 6}]. Ich vermute Counts[], dass eine eingebaute Funktion im Vergleich zu Select[], die mit beliebigen Prädikaten arbeiten muss , stark optimiert ist.
David Zhang
1

Wahrscheinlichkeitsexperimente in Mathematica durchführen

Mathematica bietet einen sehr komfortablen Rahmen für die Arbeit mit Wahrscheinlichkeiten und Verteilungen, und obwohl das Hauptproblem der angemessenen Grenzwerte angesprochen wurde, möchte ich diese Frage verwenden, um dies klarer und möglicherweise als Referenz zu verwenden.

Lassen Sie uns die Experimente einfach wiederholbar machen und einige Handlungsoptionen definieren, die unserem Geschmack entsprechen:

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

Arbeiten mit parametrischen Verteilungen

Wir können nun die asymptotische Verteilung für definieren ein Ereignis , das ist der Anteil der Köpfe in führt eine (fair) coin:π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
        ]
    }
]

Was uns die Darstellung der diskreten Verteilung der Proportionen gibt: TheoreticalDistributionPlot

Wir können die Verteilung sofort verwenden, um Wahrscheinlichkeiten für und zu berechnen. :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}

Monte-Carlo-Experimente durchführen

Wir können die Verteilung für ein Ereignis verwenden, um wiederholt eine Stichprobe daraus zu erstellen (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

Ein Vergleich mit der theoretischen / asymptotischen Verteilung zeigt, dass alles ziemlich gut passt:

Show @ {
   theoreticalPlot,
   empiricalPlot
}

Verteilung vergleichen

gwr
quelle
Sie finden einen ähnlichen Beitrag mit weiteren Hintergrundinformationen zu Mathematica auf Mathematica SE .
GWR