Łamigłówka: Jak wygenerować 7 liczb całkowitych z jednakowym prawdopodobieństwem przy użyciu stronniczej monety, która ma pr (głowa) = p?

Oto pytanie, które znalazłem na Glassdoor : Jak wygenerować 7 liczb całkowitych z jednakowym prawdopodobieństwem, używając monety, która ma $\mathbb{Pr}(\text{Head}) = p\in(0,1)$ ?

Zasadniczo masz monetę, która może, ale nie musi być uczciwa, i jest to jedyny proces generowania liczb losowych, jaki masz, więc wymyśl generator liczb losowych, który generuje liczby całkowite od 1 do 7, w których prawdopodobieństwo uzyskania każdej z tych liczb całkowitych wynosi 1/7.

Wydajność danych generuje ważne procesy.

Rzuć monetą dwa razy. Jeśli wyląduje HHlub TTzignoruj ​​go i dwukrotnie odwróć.

Teraz moneta ma jednakowe prawdopodobieństwo pojawienia się HTlub TH. Jeśli się pojawi HT, zadzwoń H1. Jeśli się pojawi TH, zadzwoń T1.

Zdobywaj kolejne H1lub T1dopóki nie będziesz mieć trzech z rzędu. Te trzy wyniki dają liczbę na podstawie poniższej tabeli:

H1 H1 H1 -> 1
H1 H1 T1 -> 2
H1 T1 H1 -> 3
H1 T1 T1 -> 4
T1 H1 H1 -> 5
T1 H1 T1 -> 6
T1 T1 H1 -> 7
T1 T1 T1 -> [Throw out all results so far and repeat]

Twierdzę, że to działałoby idealnie dobrze, chociaż miałbyś wiele zmarnowanych rzutów!

Załóżmy, że .$p \in (0,1)$

Etap 1: . Rzuć monetą 5 razy.

Jeśli wynik jest

, zwróć 1 i zatrzymaj.$(H, H, H, T, T)$$1$

, zwróć 2 i zatrzymaj się.$(H, H, T, T, H)$$2$

, zwróć 3 i zatrzymaj.$(H, T, T, H, H)$$3$

, zwróć 4 i zatrzymaj.$(T, T, H, H, H)$$4$

, zwróć 5 i zatrzymaj.$(T, H, H, H, T)$$5$

, wróć 6 i zatrzymaj się.$(H, H, T, H, T)$$6$

, zwróć 7 i zatrzymaj się.$(H, T, H, T, H)$$7$

Etap 2: . Jeśli wynik nie jest żadnym z powyższych, powtórz krok 1.

Pamiętaj, że niezależnie od wartości , każdy z siedmiu wymienionych powyżej wyników ma prawdopodobieństwo q = p 3 ( 1 - p ) 2 , a oczekiwana liczba rzutów monetą wynosi $p \in (0,1)$$q = p^3(1-p)^2$ . Rzutnik nie musi znać wartościp(z wyjątkiem tego, żep≠0ip≠1); gwarantowane jest, że siedem liczb całkowitych zostanie jednakowo zwróconych przez eksperyment po jego zakończeniu (i gwarantuje się, że zakończy się z prawdopodobieństwem1).$\displaystyle \frac{5}{7q}$$p$$p\neq 0$$p\neq 1$$1$

Uogólnienie przypadku opisanego przez Dilip Sarwate

Niektóre metody opisane w innych odpowiedziach wykorzystują schemat, w którym rzucasz sekwencją $n$ monet w „turze” iw zależności od wyniku wybierasz liczbę od 1 do 7 lub odrzucasz turę i rzucasz ponownie.

Sztuką jest znalezienie w rozszerzeniu możliwości wielu 7 wyników o tym samym prawdopodobieństwie $p^k(1-p)^{n-k}$ i dopasowanie ich do siebie.

Ponieważ łączna liczba wyników nie jest wielokrotnością liczby 7, mamy kilka wyników, których nie możemy przypisać do liczby, i istnieje pewne prawdopodobieństwo, że musimy odrzucić wyniki i zacząć od nowa.

---

Przypadek użycia 7 rzutów monetą na turę

Intuicyjnie możemy powiedzieć, że rzucie kostką siedem razy byłoby bardzo interesujące. Ponieważ musimy wyrzucić tylko 2 z $2^7$ możliwości. Mianowicie 7 razy głowy i 0 razy głowy.

Dla wszystkich pozostałych $2^7-2$ możliwości zawsze występuje wiele 7 skrzynek z tą samą liczbą głowic. Mianowicie 7 skrzyń z 1 głowami, 21 skrzyń z 2 głowami, 35 skrzynek z 3 głowami, 35 skrzynek z 4 głowami, 21 skrzynek z 5 głowami i 7 skrzynek z 6 głowami.

Więc jeśli policzysz liczbę (odrzucając 0 głów i 7 głów)

$$X = \sum_{k=1}^{7} (k-1) \cdot C_k$$

z $C_k$ Bernoulliego rozproszone zmiennych (wartość 0 lub 1), wówczas X modulo 7 jest jednolity zmiennej z siedmiu możliwych wyników.

---

Porównywanie różnej liczby rzutów monetą na turę

Pozostaje pytanie, jaka byłaby optymalna liczba rzutów na turę. Rzucanie większą ilością kości na turę kosztuje więcej, ale zmniejszasz prawdopodobieństwo ponownego rzutu.

Poniższy obraz pokazuje ręczne obliczenia dla pierwszych kilku liczb rzutów monetą na turę. (być może istnieje rozwiązanie analityczne, ale uważam, że można bezpiecznie powiedzieć, że system z 7 rzutami monetą stanowi najlepszą metodę pod względem wartości oczekiwanej dla wymaganej liczby rzutów monetą)

# plot an empty canvas
plot(-100,-100,
     xlab="flips per turn",
     ylab="E(total flips)",
     ylim=c(7,400),xlim=c(0,20),log="y")
title("expectation value for total number of coin flips
(number of turns times flips per turn)")

# loop 1
# different values p from fair to very unfair 
# since this is symmetric only from 0 to 0.5 is necessary 

# loop 2
# different values for number of flips per turn
# we can only use a multiple of 7 to assign 
#   so the modulus will have to be discarded
#   from this we can calculate the probability that the turn succeeds
#   the expected number of flips is 
#       the flips per turn 
#             divided by 
#       the probability for the turn to succeed 

for (p in c(0.5,0.2,0.1,0.05)) {
  Ecoins <- rep(0,16)
  for (dr in (5:20)){
    Pdiscards = 0
    for (i in c(0:dr)) { 
      Pdiscards = Pdiscards + p^(i)*(1-p)^(dr-i) * (choose(dr,i) %% 7)
    }
    Ecoins[dr-4] = dr/(1-Pdiscards)
  }
  lines(5:20, Ecoins)
  points(5:20, Ecoins, pch=21, col="black", bg="white", cex=0.5)
  text(5, Ecoins[1], paste0("p = ",p), pos=2)
}

---

Korzystanie z reguły wczesnego zatrzymywania

Uwaga: poniższe obliczenia, dla wartości oczekiwanej liczby rzutów, dotyczą uczciwej monety $p=0.5$ , zrobienie tego byłoby bałaganem dla różnych $p$ , ale zasada pozostaje ta sama (chociaż różne księgowanie skrzynie są potrzebne)

Powinniśmy być w stanie wybrać przypadki (zamiast formuły dla $X$ ), tak abyśmy mogli zatrzymać wcześniej.

• Z 5 rzutami monet mamy sześć możliwych różnych nieuporządkowanych zestawów głów i reszek:

  1 + 5 + 10 + 10 + 5 + 1 zamówionych zestawów

  I możemy użyć grup z dziesięcioma przypadkami (to jest grupa z 2 głowami lub grupa z 2 ogonami), aby wybrać (z jednakowym prawdopodobieństwem) liczbę. Dzieje się tak w 14 z 2 ^ 5 = 32 przypadków. To pozostawia nam:

  1 + 5 + 3 + 3 + 5 + 1 zamówionych zestawów

• Z dodatkowym (6-tym) rzutem monetą mamy siedem możliwych różnych nieuporządkowanych zestawów głów i reszek:

  1 + 6 + 8 + 6 + 8 + 6 + 1 zamówionych zestawów

  I możemy użyć grup z ośmioma przypadkami (to jest grupa z 3 głowami lub grupa z 3 ogonami), aby wybrać (z jednakowym prawdopodobieństwem) liczbę. Dzieje się tak w 14 na 2 * (2 ^ 5-14) = 36 przypadków. To pozostawia nam:

  1 + 6 + 1 + 6 + 1 + 6 + 1 zamówionych zestawów

• Z kolejnym (7-tym) dodatkowym rzutem monetą mamy osiem możliwych różnych nieuporządkowanych zestawów głów i reszek:

  1 + 7 + 7 + 7 + 7 + 7 + 7 + 1 zamówionych zestawów

  I możemy użyć grup z siedmioma przypadkami (wszystkie z wyjątkiem wszystkich ogonów i wszystkich głów), aby wybrać (z jednakowym prawdopodobieństwem) liczbę. Dzieje się tak w 42 z 44 przypadków. To pozostawia nam:

  1 + 0 + 0 + 0 + 0 + 0 + 0 + 1 zamówionych zestawów

  (moglibyśmy to kontynuować, ale tylko w 49 kroku daje to nam przewagę)

Prawdopodobieństwo wyboru liczby

• przy 5 rzutach jest $\frac{14}{32} = \frac{7}{16}$
• przy 6 rzutach jest $\frac{9}{16}\frac{14}{36} = \frac{7}{32}$
• przy 7 rzutach jest $\frac{11}{32}\frac{42}{44} = \frac{231}{704}$
• nie w 7 rzutach jest $1-\frac{7}{16}-\frac{7}{32}-\frac{231}{704} = \frac{2}{2^7}$

To sprawia, że ​​wartość oczekiwana dla liczby przerzutów w jednej turze jest uzależniona od powodzenia, a p = 0,5:

$$5 \cdot \frac{7}{16}+ 6 \cdot \frac{7}{32} + 7 \cdot \frac{231}{704} = 5.796875$$

Wartość oczekiwana dla całkowitej liczby przewrotów (do momentu sukcesu), pod warunkiem, że p = 0,5, staje się:

$$\left(5 \cdot \frac{7}{16}+ 6 \cdot \frac{7}{32} + 7 \cdot \frac{231}{704}\right) \frac{2^7}{2^7-2} = \frac{53}{9} = 5.88889$$

---

Odpowiedź NcAdams wykorzystuje odmianę strategii zatrzymania (za każdym razem wymyślają dwa nowe rzuty monetą), ale nie optymalnie wybiera wszystkie rzuty.

Odpowiedź Clida może być podobna, chociaż może istnieć nierówna zasada wyboru, że każda z dwóch monet może zostać wybrana, ale niekoniecznie z jednakowym prawdopodobieństwem (rozbieżność, która jest naprawiana podczas późniejszych rzutów monetą)

---

Porównanie z innymi metodami

Inne metody wykorzystujące podobną zasadę to NcAdams i AdamO.

Zasada jest następująca : decyzja o liczbie od 1 do 7 jest podejmowana po określonej liczbie głów i ogonów. Po liczbie $x$ przewrotów dla każdej decyzji prowadzącej do liczby $i$ istnieje podobna, równie prawdopodobna decyzja, która prowadzi do liczby $j$ (ta sama liczba głów i ogonów, ale tylko w innej kolejności). Niektóre serie głów i ogonów mogą prowadzić do podjęcia decyzji o rozpoczęciu od nowa.

W przypadku tego rodzaju metod, ta umieszczona tutaj jest najbardziej wydajna, ponieważ podejmuje decyzje tak wcześnie, jak to możliwe (tak szybko, jak istnieje możliwość 7 równych sekwencji prawdopodobieństwa głów i ogonów, po $x$ przerzucie możemy użyć aby podjąć decyzję w sprawie liczby i nie musimy przewracać dalej, jeśli napotkamy jeden z tych przypadków).

Dowodzi tego poniższy obraz i symulacja:

#### mathematical part #####
set.seed(1)


#plotting this method
p <- seq(0.001,0.999,0.001)
tot <- (5*7*(p^2*(1-p)^3+p^3*(1-p)^2)+
       6*7*(p^2*(1-p)^4+p^4*(1-p)^2)+
       7*7*(p^1*(1-p)^6+p^2*(1-p)^5+p^3*(1-p)^4+p^4*(1-p)^3+p^5*(1-p)^2+p^6*(1-p)^1)+
        7*1*(0+p^7+(1-p)^7) )/
             (1-p^7-(1-p)^7)
plot(p,tot,type="l",log="y",
     xlab="p",
     ylab="expactation value number of flips"
     )

#plotting method by AdamO
tot <- (7*(p^20-20*p^19+189*p^18-1121*p^17+4674*p^16-14536*p^15+34900*p^14-66014*p^13+99426*p^12-119573*p^11+114257*p^10-85514*p^9+48750*p^8-20100*p^7+5400*p^6-720*p^5)+6*
          (-7*p^21+140*p^20-1323*p^19+7847*p^18-32718*p^17+101752*p^16-244307*p^15+462196*p^14-696612*p^13+839468*p^12-806260*p^11+610617*p^10-357343*p^9+156100*p^8-47950*p^7+9240*p^6-840*p^5)+5*
          (21*p^22-420*p^21+3969*p^20-23541*p^19+98154*p^18-305277*p^17+733257*p^16-1389066*p^15+2100987*p^14-2552529*p^13+2493624*p^12-1952475*p^11+1215900*p^10-594216*p^9+222600*p^8-61068*p^7+11088*p^6-1008*p^5)+4*(-
          35*p^23+700*p^22-6615*p^21+39235*p^20-163625*p^19+509425*p^18-1227345*p^17+2341955*p^16-3595725*p^15+4493195*p^14-4609675*p^13+3907820*p^12-2745610*p^11+1592640*p^10-750855*p^9+278250*p^8-76335*p^7+13860*p^6-
          1260*p^5)+3*(35*p^24-700*p^23+6615*p^22-39270*p^21+164325*p^20-515935*p^19+1264725*p^18-2490320*p^17+4027555*p^16-5447470*p^15+6245645*p^14-6113275*p^13+5102720*p^12-3597370*p^11+2105880*p^10-999180*p^9+371000
           *p^8-101780*p^7+18480*p^6-1680*p^5)+2*(-21*p^25+420*p^24-3990*p^23+24024*p^22-103362*p^21+340221*p^20-896679*p^19+1954827*p^18-3604755*p^17+5695179*p^16-7742301*p^15+9038379*p^14-9009357*p^13+7608720*p^12-
           5390385*p^11+3158820*p^10-1498770*p^9+556500*p^8-152670*p^7+27720*p^6-2520*p^5))/(7*p^27-147*p^26+1505*p^25-10073*p^24+49777*p^23-193781*p^22+616532*p^21-1636082*p^20+3660762*p^19-6946380*p^18+11213888*p^17-
           15426950*p^16+18087244*p^15-18037012*p^14+15224160*p^13-10781610*p^12+6317640*p^11-2997540*p^10+1113000*p^9-305340*p^8+55440*p^7-5040*p^6)
lines(p,tot,col=2,lty=2)

#plotting method by NcAdam
lines(p,3*8/7/(p*(1-p)),col=3,lty=2)

legend(0.2,500,
       c("this method calculation","AdamO","NcAdams","this method simulation"),
       lty=c(1,2,2,0),pch=c(NA,NA,NA,1),col=c(1,2,3,1))


##### simulation part ######

#creating decision table
mat<-matrix(as.numeric(intToBits(c(0:(2^5-1)))),2^5,byrow=1)[,c(1:12)]
colnames(mat) <- c("b1","b2","b3","b4","b5","b6","b7","sum5","sum6","sum7","decision","exit")

# first 5 rolls
mat[,8] <- sapply(c(1:2^5), FUN = function(x) {sum(mat[x,1:5])})

mat[which((mat[,8]==2)&(mat[,11]==0))[1:7],12] = rep(5,7) # we can stop for 7 cases with 2 heads
mat[which((mat[,8]==2)&(mat[,11]==0))[1:7],11] = c(1:7)   
mat[which((mat[,8]==3)&(mat[,11]==0))[1:7],12] = rep(5,7) # we can stop for 7 cases with 3 heads
mat[which((mat[,8]==3)&(mat[,11]==0))[1:7],11] = c(1:7)    

# extra 6th roll
mat <- rbind(mat,mat)
mat[c(33:64),6] <- rep(1,32)
mat[,9] <- sapply(c(1:2^6), FUN = function(x) {sum(mat[x,1:6])})

mat[which((mat[,9]==2)&(mat[,11]==0))[1:7],12] = rep(6,7) # we can stop for 7 cases with 2 heads
mat[which((mat[,9]==2)&(mat[,11]==0))[1:7],11] = c(1:7)   
mat[which((mat[,9]==4)&(mat[,11]==0))[1:7],12] = rep(6,7) # we can stop for 7 cases with 4 heads
mat[which((mat[,9]==4)&(mat[,11]==0))[1:7],11] = c(1:7)    

# extra 7th roll
mat <- rbind(mat,mat)
mat[c(65:128),7] <- rep(1,64)
mat[,10] <- sapply(c(1:2^7), FUN = function(x) {sum(mat[x,1:7])})

for (i in 1:6) {
  mat[which((mat[,10]==i)&(mat[,11]==0))[1:7],12] = rep(7,7) # we can stop for 7 cases with i heads
  mat[which((mat[,10]==i)&(mat[,11]==0))[1:7],11] = c(1:7)   
}


mat[1,12] = 7           # when we did not have succes we still need to count the 7 coin tosses
mat[2^7,12] = 7


draws = rep(0,100)
num = rep(0,100)
# plotting simulation
for (p in seq(0.05,0.95,0.05)) {
  n <- rep(0,1000)
  for (i in 1:1000) {
    coinflips <- rbinom(7,1,p)  # draw seven numbers
    I <- mat[,1:7]-matrix(rep(coinflips,2^7),2^7,byrow=1) == rep(0,7)                      # compare with the table
    Imatch = I[,1]*I[,2]*I[,3]*I[,4]*I[,5]*I[,6]*I[,7]        # compare with the table 
      draws[i] <- mat[which(Imatch==1),11]                 # result which number
      num[i]   <- mat[which(Imatch==1),12]                 # result how long it took
  }
  Nturn <- mean(num)                   #how many flips we made
  Sturn <- (1000-sum(draws==0))/1000   #how many numbers we got (relatively)
  points(p,Nturn/Sturn)
}

inny obraz skalowany przez $p*(1-p)$ dla lepszego porównania:

powiększyć metody porównywania opisane w tym poście i komentarzach

„warunkowe pominięcie siódmego kroku” to niewielka poprawa, którą można wprowadzić na zasadzie wczesnego zatrzymania. W takim przypadku nie wybierasz grup z jednakowymi prawdopodobieństwami po 6-tym rzucie. Masz 6 grup z jednakowymi prawdopodobieństwami i 1 grupy z nieco innym prawdopodobieństwem (dla tej ostatniej grupy musisz przerzucić jeszcze jeden dodatkowy czas, gdy masz 6 głów lub ogonów, a ponieważ odrzucisz 7 głów lub 7 ogonów, skończysz w końcu z takim samym prawdopodobieństwem)

---

Napisane przez StackExchangeStrike

Podziel pole na siedem regionów o równej powierzchni, każdy oznaczony liczbą całkowitą. Wrzuć monetę do pudełka w taki sposób, aby prawdopodobieństwo wylądowania było równe w każdym regionie.

To tylko połowa żartu - jest to w zasadzie ta sama procedura, co oszacowanie $\pi$ za pomocą fizycznej procedury Monte Carlo, na przykład upuszczenie ziaren ryżu na papier z narysowanym kółkiem.

Jest to jedna z jedynych odpowiedzi, która działa w przypadku lub p = 0 .$p = 1$$p=0$

EDYCJA: w oparciu o opinie innych.

Oto interesująca myśl:

ustaw listę {1,2,3,4,5,6,7}. Rzuć monetą kolejno dla każdego elementu na liście. Jeśli wyląduje głową do góry dla określonego elementu, usuń numer z listy. Jeśli wszystkie liczby z określonej iteracji listy zostaną usunięte, powtórz próbkowanie. Rób tak, aż pozostanie tylko jedna liczba.

drop.one <- function(x, p) {
  drop <- runif(length(x)) < p
  if (all(drop))
    return(x)
  return(x[!drop])
}

sample.recur <- function(x, p) {
  if (length(x) > 1)
    return(sample.recur(drop.one(x, p), p))
  return(x)
}

# x <- c(1:7,7:1)
x <- 1:7
p <- 0.01

out <- replicate(1e5, sample.recur(x, p))

round(prop.table(table(out)), 2)

daje mi w przybliżeniu równomierny rozkład

> round(prop.table(table(out)), 2)
out
   1    2    3    4    5    6    7 
0.14 0.14 0.15 0.14 0.14 0.14 0.14 

$N$

---

Ocena wartości oczekiwanej liczby rzutów monetą

$x$$y$

$$M = \begin{bmatrix} q^7 & 0 & 0 & 0 &0 & 0 & 1 & 1\\ 7p^1q^6 & q^6 & 0 & 0 & 0 & 0 & 0 & 0 \\ 21p^2q^5 & 6p^1q^5 & q^5 & 0 & 0 & 0 & 0 & 0\\ 35 p^3q^4 & 15 p^2q^4 & 5q^4 & q^4 & 0 & 0 & 0 & 0\\ 35 p^4q^3 & 20 p^3q^3 & 10 p^2q^3 & 4 p^1q^3 & q^3 & 0 & 0 & 0\\ 21p^5q^2 & 15 p^4q^2 & 10 p^3q^2 & 6 p^2q^2 & 3 p^1q^2 & q^2 & 0 & 0\\ 7p^6q^1 & 6 p^5q^1 & 5 p^4q^1 & 4 p^3q^1 & 3p^2q^1 & 2p^1q^1 & 0 & 0\\ p^7 & p^6 & p^5 & p^4 & p^3 & p^2 & 0 & 0 \end{bmatrix}$$

$(M-I)v=0$ ) pokazuje, ile czasu relatywnie spędza się w jakim stanie. W takim razie stan siódmy to, jak często będziesz w stanie narysować liczbę od 1 do 7. Inne stany mówią, ile kosztuje rzut monetą.

$E(n) = \frac{24}{7}p(1-p)$

$p>2/3$ . Ale także wydajność jest niesymetryczna. Symetryczna i lepsza ogólna wydajność mogłaby zostać osiągnięta, gdyby utworzona została probabilistyczna reguła przełączania, która zmienia regułę decyzji z ogonów na głowy, gdy głowy okazują się nieprawdopodobne.

Rozwiązanie znalezione w wxMaxima

M: matrix(
 [(1-p)^7,        0,          0,0,0,0,1,1], 
 [7* p*(1-p)^6,   (1-p)^6,        0,0,0,0,0,0], 
 [21*p^2*(1-p)^5, 6*p*(1-p)^5,    (1-p)^5,0,0,0,0,0], 
 [35*p^3*(1-p)^4, 15*p^2*(1-p)^4, 5*p*(1-p)^4,(1-p)^4,0,0,0,0], 
 [35*p^4*(1-p)^3, 20*p^3*(1-p)^3, 10*p^2*(1-p)^3,4*p*(1-p)^3,(1-p)^3,0,0,0], 
 [21*p^5*(1-p)^2, 15*p^4*(1-p)^2, 10*p^3*(1-p)^2,6*p^2*(1-p)^2,3*p*(1-p)^2,(1-p)^2,0,0], 
 [7* p^6*(1-p)^1, 6*p^5*(1-p),    5*p^4*(1-p),4*p^3*(1-p),3*p^2*(1-p),2*(1-p)*p,0,0], 
 [p^7,        p^6,        p^5,p^4,p^3,p^2,0,0]
);
z: nullspace(M-diagmatrix(8,1));
x : apply (addcol, args (z));
t : [7,6,5,4,3,2,0,0];
plot2d(t.x/x[7],[p,0,1],logy);

Obliczenia w R.

# plotting empty canvas
plot(-100,-100,
     xlab="p",
     ylab="E(total flips)",
     ylim=c(10,1000),xlim=c(0,1),log="y")

# plotting simulation
for (p in seq(0.1,0.9,0.05)) {

  n <- rep(0,10000)
  for (i in 1:10000) {
    success  = 0
    tests = c(1,1,1,1,1,1,1)     # start with seven numbers in the set
    count = 0
    while(success==0) {
      for (j in 1:7)  {
        if (tests[j]==1) {
          count = count + 1
          if  (rbinom(1,1,p) == 1) {
            tests[j] <- 0        # elliminate number when we draw heads
          }
        }
      }
      if (sum(tests)==1) {
        n[i] = count
        success = 1              # end     when 1 is left over
      }
      if (sum(tests)==0) {
        tests = c(1,1,1,1,1,1,1) # restart when 0 are left over
      }
    }
  }
  points(p,mean(n))
}

# plotting formula
p <- seq(0.001,0.999,0.001)

tot <- (7*(p^20-20*p^19+189*p^18-1121*p^17+4674*p^16-14536*p^15+34900*p^14-66014*p^13+99426*p^12-119573*p^11+114257*p^10-85514*p^9+48750*p^8-20100*p^7+5400*p^6-720*p^5)+6*
    (-7*p^21+140*p^20-1323*p^19+7847*p^18-32718*p^17+101752*p^16-244307*p^15+462196*p^14-696612*p^13+839468*p^12-806260*p^11+610617*p^10-357343*p^9+156100*p^8-47950*p^7+9240*p^6-840*p^5)+5*
    (21*p^22-420*p^21+3969*p^20-23541*p^19+98154*p^18-305277*p^17+733257*p^16-1389066*p^15+2100987*p^14-2552529*p^13+2493624*p^12-1952475*p^11+1215900*p^10-594216*p^9+222600*p^8-61068*p^7+11088*p^6-1008*p^5)+4*(-
    35*p^23+700*p^22-6615*p^21+39235*p^20-163625*p^19+509425*p^18-1227345*p^17+2341955*p^16-3595725*p^15+4493195*p^14-4609675*p^13+3907820*p^12-2745610*p^11+1592640*p^10-750855*p^9+278250*p^8-76335*p^7+13860*p^6-
    1260*p^5)+3*(35*p^24-700*p^23+6615*p^22-39270*p^21+164325*p^20-515935*p^19+1264725*p^18-2490320*p^17+4027555*p^16-5447470*p^15+6245645*p^14-6113275*p^13+5102720*p^12-3597370*p^11+2105880*p^10-999180*p^9+371000
   *p^8-101780*p^7+18480*p^6-1680*p^5)+2*(-21*p^25+420*p^24-3990*p^23+24024*p^22-103362*p^21+340221*p^20-896679*p^19+1954827*p^18-3604755*p^17+5695179*p^16-7742301*p^15+9038379*p^14-9009357*p^13+7608720*p^12-
 5390385*p^11+3158820*p^10-1498770*p^9+556500*p^8-152670*p^7+27720*p^6-2520*p^5))/(7*p^27-147*p^26+1505*p^25-10073*p^24+49777*p^23-193781*p^22+616532*p^21-1636082*p^20+3660762*p^19-6946380*p^18+11213888*p^17-
  15426950*p^16+18087244*p^15-18037012*p^14+15224160*p^13-10781610*p^12+6317640*p^11-2997540*p^10+1113000*p^9-305340*p^8+55440*p^7-5040*p^6)
lines(p,tot)

#plotting comparison with alternative method
lines(p,3*8/7/(p*(1-p)),lty=2)

legend(0.2,500,
       c("simulation","calculation","comparison"),
       lty=c(0,1,2),pch=c(1,NA,NA))

Pytanie jest nieco dwuznaczne, czy pyta „wygeneruj losową liczbę całkowitą równą lub mniejszą niż 7 z jednakowym prawdopodobieństwem”, czy też pytasz „wygeneruj 7 losowych liczb całkowitych z jednakowym prawdopodobieństwem?” - ale jaka jest przestrzeń liczb całkowitych?!?

Zakładam, że jest to pierwsza, ale ta sama logika, którą stosuję, może zostać rozszerzona na tę drugą sprawę, gdy problem zostanie rozwiązany.

Za pomocą stronniczej monety możesz wyprodukować uczciwą monetę, postępując zgodnie z następującą procedurą: https://en.wikipedia.org/wiki/Fair_coin#Fair_results_from_a_part_coin

Liczbę 7 lub mniejszą można zapisać dwójkowo jako trzy cyfry {0,1}. Wystarczy więc trzykrotnie wykonać powyższą procedurę i przekonwertować wyprodukowaną liczbę binarną na dziesiętną.

Rozwiązanie, które nigdy nie marnuje flipów, co bardzo pomaga w przypadku bardzo tendencyjnych monet.

Wadą tego algorytmu (przynajmniej tak jak napisano) jest to, że używa arytmetyki o dowolnej precyzji. Praktycznie prawdopodobnie chcesz tego używać, aż do przepełnienia liczb całkowitych, a dopiero potem wyrzuć go i zacznij od nowa.

Musisz także wiedzieć, na czym polega uprzedzenie ... czego możesz nie powiedzieć, jeśli jest zależne od temperatury, jak większość zjawisk fizycznych.

---

Zakładając, że szansa na głowę wynosi, powiedzmy, 30%.

• Zacznij od zakresu [1, 8).
• Rzuć monetą. Jeśli masz głowy, użyj lewego 30%, więc twój nowy zakres to [1, 3.1). W przeciwnym razie użyj właściwych 70%, więc nowy zakres to [3.1, 8).
• Powtarzaj, aż cały zakres będzie miał tę samą liczbę całkowitą.

---

Pełny kod:

#!/usr/bin/env python3
from fractions import Fraction
from collections import Counter
from random import randrange


BIAS = Fraction(3, 10)
STAT_COUNT = 100000


calls = 0
def biased_rand():
    global calls
    calls += 1
    return randrange(BIAS.denominator) < BIAS.numerator


def can_generate_multiple(start, stop):
    if stop.denominator == 1:
        # half-open range
        stop = stop.numerator - 1
    else:
        stop = int(stop)
    start = int(start)
    return start != stop


def unbiased_rand(start, stop):
    if start < 0:
        # negative numbers round wrong
        return start + unbiased_rand(0, stop - start)
    assert isinstance(start, int) and start >= 0
    assert isinstance(stop, int) and stop >= start
    start = Fraction(start)
    stop = Fraction(stop)
    while can_generate_multiple(start, stop):
        if biased_rand():
            old_diff = stop - start
            diff = old_diff * BIAS
            stop = start + diff
        else:
            old_diff = stop - start
            diff = old_diff * (1 - BIAS)
            start = stop - diff
    return int(start)


def stats(f, *args, **kwargs):
    c = Counter()
    for _ in range(STAT_COUNT):
        c[f(*args, **kwargs)] += 1

    print('stats for %s:' % f.__qualname__)
    for k, v in sorted(c.items()):
        percent = v * 100 / STAT_COUNT
        print('  %s: %f%%' % (k, percent))


def main():
    #stats(biased_rand)
    stats(unbiased_rand, 1, 7+1)
    print('used %f calls at bias %s' % (calls/STAT_COUNT, BIAS))


if __name__ == '__main__':
    main()

Jak wspomniano we wcześniejszych komentarzach, zagadka ta odnosi się do artykułu Johna von Neumanna z 1951 r. „Różne techniki stosowane w związku z przypadkowymi cyframi” opublikowanego w czasopiśmie badawczym National Bureau of Standards:

$p$$f(p)$$f$ $f(p)=\min\{1,2p\}$$N$ prób.

$p \neq 1$$p \neq 0$

Najpierw zamieniamy (prawdopodobnie) nieuczciwą monetę na monetę uczciwą, korzystając z procesu z NcAdams :

Rzuć monetą dwa razy. Jeśli wyląduje HHlub TTzignoruj ​​go i dwukrotnie odwróć.

Teraz moneta ma jednakowe prawdopodobieństwo pojawienia się HTlub TH. Jeśli się pojawi HT, zadzwoń H1. Jeśli się pojawi TH, zadzwoń T1.

$0$$1$H1$= 1$T1 $= 0$0.H1 H1 T1$0.110$

$1/7$

$1/7 = 0.001 001 001 ...$

$2/7 = 0.010 010 010 ...$

$3/7 = 0.011 011 011 ...$

$4/7 = 0.100 100 100 ...$

$5/7 = 0.101 101 101 ...$

$6/7 = 0.110 110 110 ...$

$n$$n/7$$(n-1)/7$$1$$7$

Zainspirowany odpowiedzią AdamO, oto rozwiązanie Python, które pozwala uniknąć stronniczości:

def roll(p, n):
    remaining = range(1,n+1)
    flips = 0
    while len(remaining) > 1:
        round_winners = [c for c in remaining if random.choices(['H','T'], [p, 1.0-p]) == ['H']]
        flips += len(remaining)
        if len(round_winners) > 0:
            remaining = round_winners
        p = 1.0 - p
    return remaining[0], flips

Istnieją tutaj dwie główne zmiany: Główna polega na tym, że jeśli wszystkie liczby zostaną odrzucone w rundzie, powtórz rundę. Również zmieniam wybór, czy głowy czy ogony oznaczają odrzucanie za każdym razem. Zmniejsza to liczbę potrzebnych przerzutów w przypadkach, gdy p jest bliskie 0 lub 1 o ~ 70%, gdy p = 0,999

Wygląda na to, że za każdym razem możemy zmieniać mapowanie wyniku każdego przerzucenia . Tak więc, używając dla wygody pierwszych siedmiu liczb całkowitych dodatnich, wydajemy następujące zamówienia:

$H \to 1$
$H \to 2$

$H \to 7$
$H \to 1$

itp

$T$

$AP$$T$

$$P_{AP}(\text{no integers generated}) = (1-p)^7$$

$N_b$

$$\text{Count}_{AP}(\text{useless flips}) \to 7\cdot N_b(1-p)^7$$

$B(p,n=5)$$p^3(1-p)^2$

$$P_{DS}(\text{no integers generated}) = 1-7\cdot p^3(1-p)^2$$

Liczba bezużytecznych przewrotów będzie tutaj miała tendencję do

$$\text{Count}_{DS}(\text{useless flips}) \to 5\cdot N_b\cdot [1-7\cdot p^3(1-p)^2]$$

$AP$

$$\text{Count}_{AP}(\text{useless flips}) < \text{Count}_{DS}(\text{useless flips})$$

$$\implies 7\cdot N_b(1-p)^7 < 5\cdot N_b\cdot [1-7\cdot p^3(1-p)^2]$$

$$\implies 7\cdot (1-p)^7 < 5\cdot [1-7\cdot p^3(1-p)^2]$$

$p>0.0467$$AP$ RNG generuje mniej bezużytecznych przewrotów.

$p$$AP$$DS$$p\approx 0.5967$

$$\frac{\text{Count}_{AP}(\text{useless flips})}{\text{Count}_{DS}(\text{useless flips})}$$

$0.67$$p=0.1$$0.3$$p=0.2$$0.127$$p=0.4$