O quebra-cabeça ‘moeda do prisioneiro lançando’: simulação organizada em R

cupom com desconto - o melhor site de cupom de desconto cupomcomdesconto.com.br


[Esteartigofoipublicadopelaprimeiravezem[Thisarticlewasfirstpublishedon Variação explicada, e gentilmente contribuiu para os R-blogueiros]. (Você pode relatar um problema sobre o conteúdo desta página aqui)


Deseja compartilhar seu conteúdo com R-blogueiros? clique aqui se você tiver um blog ou aqui se não tiver.

Anteriormente nesta série

  • O quebra-cabeça “cartão de embarque perdido”
  • O quebra-cabeça do “jogo mortal”
  • Puzle O cavaleiro em um tabuleiro de xadrez infinito
  • O quebra-cabeça “maior lucro ou perda de ações”
  • O quebra-cabeça “paradoxo do aniversário”
  • O quebra-cabeça “Spelling Bee honeycomb”
  • O quebra-cabeça “jogar moedas” de Feller
  • O quebra-cabeça “comentários de spam”

Adoro a coluna Riddler do 538 e gostei de resolver o quebra-cabeça de 1º de maio. Vou citar:

Você está trancado na masmorra de um castelo distante com três companheiros de prisão (ou seja, há quatro presos no total), cada um em uma cela separada, sem meios de comunicação. Mas acontece que todos vocês são lógicos (é claro)….

Cada prisioneiro receberá uma moeda justa, que pode ser lançada razoavelmente uma vez ou devolvida aos guardas sem ser lançada. Se todas as moedas lançadas aparecerem, você será libertado! Mas se qualquer uma das moedas lançadas aparecer, ou se ninguém escolher jogar uma moeda, todos estarão condenados a passar o resto de suas vidas na masmorra do castelo.

As únicas ferramentas que você e seus companheiros de prisão têm para ajudá-lo são geradores de números aleatórios, o que dará a cada preso um número aleatório, escolhido de maneira uniforme e independente entre zero e um.

Quais são suas chances de ser libertado?

Eu resolvo isso com uma simulação organizada em R, em particular usando uma das minhas funções favoritas, a função tidyr crossing(). Em um apêndice, mostrarei como obter uma solução de formulário fechado para N = 4.

Também publiquei um screencast de 30 minutos sobre como me aproximei da simulação e visualização.

Simulando quatro prisioneiros

Antes de pularmos para a nossa simulação, podemos começar com um pouco de lógica. Os quatro prisioneiros não conseguem se comunicar e estão em situações simétricas. Isso significa que todos teremos que seguir a mesma estratégia, confiando no fato de que todos os outros lógicos escolherão a mesma.

Se todos nós decidíssemos não jogar uma moeda, nunca seríamos livres. Se todos decidíssemos jogar a moeda, nossa chance de liberdade seria frac {1} {2 ^ 4} = frac {1} {16}, a chance de quatro caras. Portanto, temos um botão para controlar nossa estratégia: a probabilidade de decidirmos jogar nossa moeda em vez de devolvê-la aos guardas.

Como essa probabilidade será a mesma nos quatro prisioneiros, podemos simular isso para muitas estratégias possíveis entre 1% e 100% usando rbinom().

library(tidyverse)
library(scales)
theme_set(theme_light())
set.seed(2020-05-04)

sim <- crossing(trial = 1:100000,
                probability = seq(.01, 1, .01)) %>%
  mutate(num_flips = rbinom(n(), 4, probability),
         num_tails = rbinom(n(), num_flips, .5),
         set_free = num_flips !=  & num_tails == )

sim
## # A tibble: 10,000,000 x 5
##    trial probability num_flips num_tails set_free
##                         
##  1     1      0.01           0         0 FALSE   
##  2     1      0.02           0         0 FALSE   
##  3     1      0.03           0         0 FALSE   
##  4     1      0.04           0         0 FALSE   
##  5     1      0.05           0         0 FALSE   
##  6     1      0.06           1         0 TRUE    
##  7     1      0.0700         0         0 FALSE   
##  8     1      0.08           0         0 FALSE   
##  9     1      0.09           1         0 TRUE    
## 10     1      0.10           0         0 FALSE   
## # … with 9,999,990 more rows

O exemplo acima realiza 10 milhões de simulações (100.000 tentativas para cada probabilidade), mas como é vetorizado é muito rápido: cerca de 1,5 segundos na minha máquina. Observe que os prisioneiros são libertados somente se jogarem pelo menos uma moeda e não receberem coroa (num_flips != 0 & num_tails == 0) Nas dez primeiras simulações acima, a estratégia pareceu funcionar duas vezes (a 6ª e a 9ª observações).

Leia Também  Entendendo a aritmética de incorporação de palavras: Por que não há uma resposta única para "Rei - homem + mulher =?"

Cada valor de probability é uma estratégia: a probabilidade de cada prisioneiro decidir usar para ver se eles jogam sua moeda. Podemos assim resumir e visualizar a chance de liberdade dentro de cada probabilidade.

summarized <- sim %>%
  group_by(probability) %>%
  summarize(pct_free = mean(set_free)) 

summarized %>%
  ggplot(aes(probability, pct_free)) +
  geom_line() +
  expand_limits(y = )

Centro

Essa curva faz algum sentido intuitivo. Se a probabilidade de jogar for muito baixa, há um risco alto de que ninguém jogue uma moeda, mas se a probabilidade for muito alta, ela se aproxima de frac {1} {16}. Também sabíamos que o pico não poderia estar acima de 50%, pois pelo menos uma moeda precisará ser lançada.

summarized %>%
  arrange(desc(pct_free))
## # A tibble: 100 x 2
##    probability pct_free
##              
##  1        0.35    0.286
##  2        0.37    0.286
##  3        0.36    0.285
##  4        0.3     0.284
##  5        0.34    0.284
##  6        0.32    0.283
##  7        0.33    0.283
##  8        0.31    0.282
##  9        0.38    0.282
## 10        0.39    0.282
## # … with 90 more rows

Parece que o ideal é de cerca de 35% (embora exista alguma incerteza), e que quando eles usarem essa estratégia, os prisioneiros terão 28% de chance de libertação.

Solução exata com otimização

Para passar de uma simulação para uma solução exata, vamos começar obtendo a fórmula exata para essa curva. Qual a probabilidade de os presos ficarem livres se a chance de cada um deles jogar sua moeda for p?

Existem quatro maneiras pelas quais acabamos vencendo: 1 prisioneiro pode virar 1 cabeça, 2 prisioneiros podem virar 2 cabeças, 3 prisioneiros podem virar 3 cabeças e 4 prisioneiros podem virar 4 cabeças. Esses são eventos disjuntos (é impossível que dois deles aconteçam juntos), para que possamos resumir as probabilidades. Deixei F ser o número de moedas lançadas e T seja o número de caudas viradas. A probabilidade de obter liberdade é

Leia Também  desafio # explicarCovid19 | R-bloggers

sum_ {k = 1} ^ 4 {P (F = k | p) P (T = 0 | F = k)}

Ambas as probabilidades seguem uma distribuição binomial: a probabilidade de algum número de sucessos em um conjunto de tentativas de identificação. E a probabilidade de não haver caudas é frac {1} {2 ^ k}. Isso pode ser escrito em R da seguinte maneira, usando o dbinom() função:

probability_exact <- function(p, n = 4) {
  sum(dbinom(1:n, n, p) / 2 ^ (1:n))
}

# Probability all heads if each player has 20% chance
probability_exact(.2)
## [1] 0.2465

Poderíamos adicionar esses valores exatos em nossa simulação anterior para verificar nossa matemática.

# map_dbl lets us calculate probability of freedom for each strategy 
summarized %>%
  mutate(exact = map_dbl(probability, probability_exact)) %>%
  ggplot(aes(probability, pct_free)) +
  geom_line() +
  geom_line(aes(y = exact), color = "red", lty = 2) +
  expand_limits(y = )

Centro

Isso corresponde à simulação, então parece que acertamos!

Estamos especialmente interessados ​​no pico: qual é a estratégia ideal e a probabilidade correspondente de se libertar? Podemos usar o built-in optimize , criada para otimização unidimensional dentro de um intervalo.

opt <- optimize(probability_exact, c(, 1), maximum = TRUE)
opt
## $maximum
## [1] 0.3420391
## 
## $objective
## [1] 0.2848424

A maior chance de fuga é de 28,5%, quando os prisioneiros usam o gerador de números aleatórios para ter uma chance de 34,2% de jogar a moeda.

Se você quiser ver algumas equações em vez de simulações, o Apêndice abaixo mostra como calcular a forma exata (um pouco confusa) e obtém algumas dicas sobre como é a aparência de um N. arbitrário.

Crédito extra: N arbitrário

Crédito extra: em vez de quatro presos, suponha que haja N presos. Agora, quais são suas chances de ser libertado?

O que é maravilhoso sobre o crossing() A função é que sempre podemos adicionar outra variável ao nosso cálculo. Vamos adicionar n, variando de 2 prisioneiros a 8 prisioneiros.

sim_n <- crossing(trial = 1:100000,
                  probability = seq(.02, 1, .02),
                  n = 2:8) %>%
  mutate(num_flips = rbinom(n(), n, probability),
         num_tails = rbinom(n(), num_flips, .5),
         set_free = num_flips !=  & num_tails == )

Uma vez que nossa probability_exact() A função recebe dois argumentos (p e n), também podemos calcular todas as probabilidades exatas com map2_dbl.

cupom com desconto - o melhor site de cupom de desconto cupomcomdesconto.com.br
probabilities_n <- sim_n %>%
  group_by(probability, n) %>%
  summarize(simulated = mean(set_free)) %>%
  ungroup() %>%
  mutate(exact = map2_dbl(probability, n, probability_exact))

probabilities_n %>%
  ggplot(aes(probability, exact, color = factor(n))) +
  geom_line() +
  geom_point(aes(y = simulated), size = .4) +
  scale_x_continuous(labels = percent) +
  scale_y_continuous(labels = percent) +
  labs(x = "p: Probability of flipping the coin",
       y = "Probability of freedom",
       color = "N: # of prisoners",
       title = "What's the chance of escaping with n prisoners?",
       subtitle = "Points show simulations of 100,000 prisoners each; lines are exact solution")

Centro

Isso permite verificar nossos resultados em busca de um $ N $ arbitrário. Parece que nossa simulação e cálculos exatos estão alinhados, o que é uma boa maneira de verificar nosso trabalho!

Leia Também  Salvando gráficos R entre sistemas operacionais

A probabilidade de liberdade tem um pico para qualquer valor de $ N $. Parece que o melhor valor de p quando há dois prisioneiros tem cerca de 2/3 de chance de virar, e isso diminui à medida que N aumenta. A chance de sucesso também diminui à medida que N aumenta, mas não muito, e parece que isso pode ser uma hipótese.

Vamos usar optimize para encontrar a melhor estratégia para cada valor de N, até (digamos) 60 prisioneiros.

# Function that takes n and runs the optimise step
optimize_for_n <- function(n) {
  optimize(function(p) probability_exact(p, n), c(, 1), maximum = TRUE)
}

optimal_n <- tibble(n = 2:60) %>%
  mutate(optimal = map(n, optimize_for_n)) %>%
  unnest_wider(optimal)

optimal_n %>%
  gather(metric, value, -n) %>%
  mutate(metric = ifelse(metric == "maximum", "Optimal probability to flip", "Probability of escape")) %>%
  ggplot(aes(n, value, color = metric)) +
  geom_line() +
  geom_hline(lty = 2, yintercept = .25) +
  scale_y_continuous(labels = percent) +
  expand_limits(y = ) +
  labs(x = "N: # of prisoners",
       y = "Probability",
       color = "",
       title = "How does the optimal strategy and outcome change with N?")

Centro

O p ótimo de fato diminui à medida que N aumenta e parece estar se aproximando de zero. A probabilidade de fuga (você prefere jogar este jogo com apenas mais um prisioneiro do que com muitos), mas observe que ele está se aproximando de uma assíntota, que parece ser de 25% (mostrada como uma linha tracejada).

Em vez de pensar no valor ideal de p, pode fazer sentido pensar em Np, o número esperado de lançamentos. Ou seja, quantos disparos você está buscando em todos os N prisioneiros?

optimal_n %>%
  arrange(desc(n)) %>%
  mutate(expected_coins_flipped = n * maximum) %>%
  ggplot(aes(n, expected_coins_flipped)) +
  geom_line() +
  labs(y = "Expected # of coins getting flipped")

Centro

Com um N pequeno, você deseja lançar cerca de 1,34 moedas e, à medida que N cresce, esse alvo parece se aproximar de 1,386.

Faz sentido intuitivo que você esteja visando um número de lançamentos um pouco acima de 1. Você realmente não quer acabar lançando zero (nesse caso, você perderá imediatamente), mas está equilibrando isso contra não querer virar demais; nesse caso, você terá um alto risco de virar pelo menos uma coroa.

À medida que N cresce, como p permanece pequeno, o número de moedas lançadas se aproximará de uma distribuição de Poisson, uma distribuição útil de contagens em que a variação é igual à média. Nesse ponto, não acaba importando o que é N, por isso faz sentido que a probabilidade de escapar assintote em vez de continuar a declinar.

# Summing the Poisson probabilities up to 1000
optimize(function(p) sum(dpois(1:1000, p) / 2 ^ (1:1000)),
         c(, 10),
         maximum = TRUE)
## $maximum
## [1] 1.386312
## 
## $objective
## [1] 0.25

Em outras palavras: mesmo se você estivesse jogando esse jogo com um bilhão de prisioneiros, ainda teria 25% de chance de escapar: cada prisioneiro teria apenas uma frac {1.386} {1.000.000.000} de jogar a moeda. Muito legal!

(Eu não tenho uma intuição de por que o número ideal esperado de inversões é de 1.386, ou por que a probabilidade de fuga se aproxima de frac {1} {4}, mas aposto que alguém que é bom com séries infinitas poderia resolver isso com base em a densidade da distribuição de Poisson).

Apêndice: Solução de formulário fechado para N = 4

Eu me concentro mais no código do que nas equações deste blog, mas pensei em tentar obter uma forma exata para a estratégia ideal. Afinal, eu teria que ser um dos prisioneiros e não ter acesso a R.

sum_ {k = 1} ^ 4 {P (F = k | p) P (T = 0 | F = k)}

= sum_ {k = 1} ^ 4 frac {{4 escolha k} p ^ k (1-p) ^ {4-k}} {2 ^ k}

= frac {4p (1-p) ^ 3} {2} + frac {6p ^ 2 (1-p) ^ 2} {4} + frac {4p ^ 3 (1-p)} {8} + frac {p ^ 4} {16}

= 2p (1-p) ^ 3 + frac {3} {2} p ^ 2 (1-p) ^ 2 + frac {1} {2} p ^ 3 (1-p) + frac {1 } {16} p ^ 4

= 2p (1-p) ^ 3 + frac {3} {2} p ^ 2 (1-p) ^ 2 + frac {1} {2} p ^ 3 (1-p) + frac {1 } {16} p ^ 4

= – frac {15} {16} p ^ 4 + frac {7} {2} p ^ 3- frac {9} {2} p ^ 2 + 2p

Podemos ver no gráfico acima que existe um máximo entre 0 e 1. Para encontrar esse máximo, pegaremos a derivada e a definiremos como zero:

0 = – frac {60} {16} p ^ 3 + frac {21} {2} p ^ 2-9p + 2

A vida é muito curta para que eu lembre como encontrar as raízes das equações cúbicas (embora eu espere que os prisioneiros possam resolver isso), mas Wolfram Alpha me diz que a única solução é:

– frac {2} {15} (- 7 + 2 ^ {1/3} +2 vezes 2 ^ {2/3})

Então … se você tinha essa solução de formulário fechado no seu cartão de bingo, parabéns. Mas, na verdade, é quase igual a 0,334, portanto, verifica com a nossa simulação! Nunca me canso de como a simulação pode ajudar a verificar nossa matemática.

Corri a solução (não mostrada) para n = 3 en = 2 e descobri que eles são, respectivamente, frac {1} {7} (6-2 sqrt {2}) e ​​ frac {2} {3}. Portanto, existem algumas dicas sobre estruturas em comum (2 ^ N-1 no denominador, uma sequência de N-1as raízes de 2), mas não é algo que planejo resolver!

var vglnk = {key: ‘949efb41171ac6ec1bf7f206d57e90b8’};

(função (d, t) {
var s = d.createElement
s.src = ‘//cdn.viglink.com/api/vglnk.js’;
var r = d.getElementsByTagName
} (documento, ‘script’));

Para Deixe um comentário para o autor, siga o link e comente no blog: Variação explicada.

R-bloggers.com oferece atualizações diárias por email sobre notícias e tutoriais do R sobre o aprendizado do R e muitos outros tópicos. Clique aqui se você deseja publicar ou encontrar um emprego em ciência da dados / R.


Deseja compartilhar seu conteúdo com R-blogueiros? clique aqui se você tiver um blog ou aqui se não tiver.



cupom com desconto - o melhor site de cupom de desconto cupomcomdesconto.com.br