Lancer un dé à 6 faces jusqu'à ce que le total

11

Voici la question:

Vous lancez un dé équitable à 6 faces de manière itérative jusqu'à ce que la somme des lancers de dés soit supérieure ou égale à M. Quelle est la moyenne et l'écart-type de la somme moins M lorsque M = 300?

Dois-je écrire un code pour répondre à ce genre de questions?

Veuillez me donner quelques indices à ce sujet. Merci!

eli
la source
1
Veuillez ajouter la [self-study]balise et lire son wiki . Ensuite, dites-nous ce que vous comprenez jusqu'à présent, ce que vous avez essayé et où vous êtes coincé. Nous vous fournirons des conseils pour vous aider à vous décoller.
gung - Rétablir Monica
2
Je soupçonne que pourrait être lu comme "très grand " car je crois que ou donnerait presque exactement le même résultat. Ce que je voudrais faire est de trouver la répartition de la somme moins . M=300MM=301M=999M
Henry

Réponses:

13

Vous pouvez certainement utiliser du code, mais je ne simulerais pas.

Je vais ignorer la partie "moins M" (vous pouvez le faire assez facilement à la fin).

Vous pouvez calculer les probabilités de manière récursive très facilement, mais la réponse réelle (avec un très haut degré de précision) peut être calculée à partir d'un simple raisonnement.

Que les rouleaux soient . Soit .S t = t i = 1 X iX1,X2,...St=i=1tXi

Que soit l'indice le plus petit où .S τMτSτM

P(Sτ=M)=P(got to M6 at τ1 and rolled a 6)+P(got to M5 at τ1 and rolled a 5)++P(got to M1 at τ1 and rolled a 1)=16j=16P(Sτ1=Mj)

entrez la description de l'image ici

De même

P(Sτ=M+1)=16j=15P(Sτ1=Mj)

P(Sτ=M+2)=16j=14P(Sτ1=Mj)

P(Sτ=M+3)=16j=13P(Sτ1=Mj)

P(Sτ=M+4)=16j=12P(Sτ1=Mj)

P(Sτ=M+5)=16P(Sτ1=M1)

Des équations similaires à la première ci-dessus pourraient alors (au moins en principe) être exécutées jusqu'à ce que vous atteigniez l'une des conditions initiales pour obtenir une relation algébrique entre les conditions initiales et les probabilités que nous voulons (ce qui serait fastidieux et pas particulièrement instructif) , ou vous pouvez construire les équations avancées correspondantes et les exécuter à partir des conditions initiales, ce qui est facile à faire numériquement (et c'est ainsi que j'ai vérifié ma réponse). Cependant, nous pouvons éviter tout cela.

Les probabilités des points sont des moyennes pondérées des probabilités précédentes; ceux-ci lisseront (géométriquement rapidement) toute variation de probabilité par rapport à la distribution initiale (toute probabilité au point zéro dans le cas de notre problème). le

Pour une approximation (très précise), nous pouvons dire que à M - 1 devraient être presque également probables au temps τ - 1 (très proche de lui), et ainsi de ce qui précède, nous pouvons écrire que les probabilités être très proche d'être dans des ratios simples, et comme ils doivent être normalisés, nous pouvons simplement écrire les probabilités.M6M1τ1

Autrement dit, nous pouvons voir que si les probabilités de commencer de à M - 1 étaient exactement égales, il y a 6 façons tout aussi probables d'arriver à M , 5 d'arriver à M + 1 , et ainsi de suite jusqu'à 1 façon d'arriver à M + 5 .M6M1MM+1M+5

Autrement dit, les probabilités sont dans le rapport 6: 5: 4: 3: 2: 1, et la somme de 1, donc elles sont triviales à noter.

Le calcul exact (jusqu'aux erreurs d'arrondi numériques cumulées) en exécutant les récursions de probabilité vers l'avant à partir de zéro (je l'ai fait dans R) donne des différences de l'ordre de .Machine$double.eps( sur ma machine) par rapport à l'approximation ci-dessus (c'est-à-dire simple le raisonnement suivant les lignes ci-dessus donne des réponses effectivement exactes , car elles sont aussi proches des réponses calculées à partir de la récursivité que nous nous attendrions à ce que les réponses exactes soient2.22e-16

Voici mon code pour cela (la plupart c'est juste l'initialisation des variables, le travail est tout en une seule ligne). Le code démarre après le premier lancer (pour éviter de mettre une cellule 0, ce qui est une petite nuisance à gérer dans R); à chaque étape, il prend la cellule la plus basse qui pourrait être occupée et avance d'un jet de dé (répartissant la probabilité de cette cellule sur les 6 cellules suivantes):

 p = array(data = 0, dim = 305)
 d6 = rep(1/6,6)
 i6 = 1:6
 p[i6] = d6
 for (i in 1:299) p[i+i6] = p[i+i6] + p[i]*d6

(nous pourrions utiliser rollapply(à partir de zoo) pour le faire plus efficacement - ou un certain nombre d'autres fonctions de ce type - mais il sera plus facile à traduire si je le maintiens explicite)

Notez qu'il d6s'agit d'une fonction de probabilité discrète sur 1 à 6, donc le code à l'intérieur de la boucle de la dernière ligne construit des moyennes pondérées courantes de valeurs antérieures. C'est cette relation qui adoucit les probabilités (jusqu'aux dernières valeurs qui nous intéressent).

Voici donc les 50 premières valeurs impaires (les 25 premières valeurs marquées de cercles). À chaque , la valeur sur l'axe des y représente la probabilité qui s'est accumulée dans la cellule la plus en arrière avant de la faire avancer dans les 6 cellules suivantes.t

entrez la description de l'image ici

Comme vous le voyez, il se lisse (à , l'inverse de la moyenne du nombre de pas que vous lance chaque dé) assez rapidement et reste constant.1/μ

MM

entrez la description de l'image ici

M1M6

MM1M6τ1M

Rt=StMR0

A partir de la distribution de probabilité, la moyenne et la variance des probabilités sont alors simples.

M

53253M=300

Glen_b -Reinstate Monica
la source
+1 Je n'ai pas bien compris cette réponse avant d'avoir développé la mienne, qui semble désormais superflue. Certains lecteurs verront peut-être la valeur des résultats de l'illustration et de la simulation, alors je garderai ma réponse ouverte.
whuber
1
@whuber Ma réponse est beaucoup moins concrète que je ne l'aurais souhaité parce que je fonctionnais sous l'hypothèse que c'était des devoirs (j'ai donc évité de faire trop de dérivation ou de donner du code - c'était plutôt destiné à servir de plan). J'ai eu du mal à écrire une réponse claire sur ce problème (c'est celui où le concret aide plus que d'habitude). Puisque vous avez donné une réponse qui contient les chiffres et le code réels (quelle réponse je pense que cela devrait rester), je sens que je peux faire certaines choses qui, je l'espère, rendront ma réponse plus facile à comprendre (soyez plus explicite, donnez mon propre code) .
Glen_b -Reinstate Monica
J'ai écrit une bien meilleure explication de ce genre de problème il y a quelques années. Si je peux me souvenir / comprendre comment cela s'est passé, j'essaierai d'en inclure une partie ici.
Glen_b -Reinstate Monica
@Glen_b a un peu compris les équations. Je suis novice. comment commencer à penser comme ça? Y a-t-il des livres que vous pourriez recommander à cet effet? Votre réponse serait très utile.
Suspect habituel du
Suspect habituel - J'ai écrit les équations en imaginant un plateau de jeu comme une longue piste et en allant "de quelle manière pourrais-je accéder à cet espace d'une manière qui s'adapte aux conditions du problème, et avec quelles chances?"; Je l'ai fait pour un espace étiqueté "M", puis pour l'espace après, et ainsi de suite. J'ai écrit le calcul similaire à l'avenir pour le code en imaginant être près de la cellule de départ et en disant "si j'étais ici, serais-je le prochain, avec quelles chances?". Les équations ne sont que des réponses à ces questions.
Glen_b -Reinstate Monica
8

Ω0nEnn

En={ωΩ|nω}.

XM(ω)ωMXMMXM

XM(ω)M{0,1,2,3,4,5}XMM=kωp(i)=1/6ii=1,2,3,4,5,6

Pr(XMM=k)=j=k6Pr(EM+kj)p(j)=16j=k6Pr(EM+kj).

M

Pr(Ei)2/7.
(1+2+3+4+5+6)/6=7/2ω

EiEi11Ei22Ei66

Pr(Ei)=j=16Pr(Eij)p(j)=16j=16Pr(Eij).

Les valeurs initiales de cette séquence sont

Pr(E0)=1;Pr(Ei)=0,i=1,2,3,.

Figure: tracé de E_i

Pr(Ei)i2/7

Pr(Ei)ith

x6p(1)x5p(2)x4p(3)x3p(6)=x6(x5+x4+x3+x2+x+1)/6.

exp(0.314368)exp(36.05)i36.05/0.314368=1152/7

M=300115EM+kj=2/7

Pr(XMM=(0,1,2,3,4,5))=(27)(16)(6,5,4,3,2,1).

Le calcul de la moyenne et de la variance de cette distribution est simple et facile.


RM+5=305X300300χ20.1367

M <- 300
n.iter <- 1e5
set.seed(17)
n <- ceiling((2/7) * (M + 3*sqrt(M)))
dice <- matrix(ceiling(6*runif(n*n.iter)), n, n.iter)
omega <- apply(dice, 2, cumsum)
omega <- omega[, apply(omega, 2, max) >= M+5]
omega[omega < M] <- NA
x <- apply(omega, 2, min, na.rm=TRUE)
count <- tabulate(x)[0:5+M]
(cbind(count, expected=round((2/7) * (6:1)/6 * length(x), 1)))
chisq.test(count, p=(2/7) * (6:1)/6)
whuber
la source