Como programar uma simulação de Monte Carlo do paradoxo da caixa de Bertrand?

12

O seguinte problema foi publicado na página do Mensa International no Facebook:

insira a descrição da imagem aqui

O post em si recebeu mais de 1000 comentários, mas não vou entrar em detalhes sobre o debate, pois sei que esse é o paradoxo da caixa de Bertrand e a resposta é . O que me interessa aqui é como responder a esse problema usando uma abordagem de Monte Carlo? Como é o algoritmo para resolver esse problema?23

Aqui está a minha tentativa:

  1. Gere números aleatórios distribuídos uniformemente entre e .N01
  2. Deixe que o evento da caixa contenha 2 bolas de ouro (caixa 1) selecionadas com menos da metade.
  3. Contar os números que menos do que e chamar o resultado como .0.5S
  4. Como é certo obter uma bola de ouro se a caixa 1 for selecionada e apenas 50% de chance de obter uma bola de ouro se a caixa 2 for selecionada, portanto, a probabilidade de obter uma sequência GG é
    P(B2=G|B1=G)=SS+0.5(NS)

Implementando o algoritmo acima em R:

N <- 10000
S <- sum(runif(N)<0.5)
S/(S+0.5*(N-S))

A saída do programa acima é de cerca de que quase corresponde à resposta correta, mas não tenho certeza de que seja o caminho correto. Existe uma maneira adequada de resolver esse problema programaticamente?0.67

Anastasiya-Romanova 秀
fonte

Respostas:

14

Como @ Henry , eu realmente não acho que sua solução seja um Monte Carlo. Certamente, você tira amostras da distribuição, mas isso não tem muito a ver com imitar o processo de geração de dados. Se você quiser usar Monte Carlo para convencer alguém de que a solução teórica está correta, será necessário usar uma solução que imite o processo de geração de dados. Eu imagino que seja algo como abaixo:

boxes <- list(
  c(0, 0),
  c(0, 1),
  c(1, 1)
)

count_successes = 0
count_valid_samples = 0

for (i in 1:5000) {
  sampled_box <- unlist(sample(boxes, 1)) # sample box
  sampled_balls <- sample(sampled_box)    # shuffle balls in the box

  if (sampled_balls[1] == 1) {            # if first ball is golden
    if (sampled_balls[2] == 1) {          # if second ball is golden
      count_successes = count_successes + 1
    }
    count_valid_samples = count_valid_samples + 1
  }
}
count_successes / count_valid_samples

ou usando o código "vetorizado":

mean(replicate(5000, {       # repeat 5000 times, next calculate empirical probability
  x <- boxes[[sample(3, 1)]] # pick a box
  if (x[sample(2, 1)] == 1)  # pick a ball, check if it is golden
    return(sum(x) == 2)      # check if you have two golden balls in the box
  else
    return(NA)               # ignore if sampled ball is silver
  }), na.rm = TRUE)          # not count if silver

Observe que, uma vez que você condiciona o fato de que a primeira bola já foi sacada e é dourada, o código acima poderia simplesmente usar duas caixas boxes <- list(c(0, 1), c(1, 1))e depois fazer uma amostra delas x <- boxes[[sample(2, 1)]], para que o código fosse mais rápido, pois não faria o 1/3 corridas vazias que descontamos. No entanto, como o problema é simples e o código é executado rapidamente, podemos simular explicitamente todo o processo de geração de dados "para ter certeza" de que o resultado está correto. Além disso, essa etapa não é necessária, pois produziria os mesmos resultados nos dois casos.

Tim
fonte
Será que x <- boxes[[sample(3, 1)]]significa que você tomar uma caixa de 3 caixas? Se sim, por que é necessário, pois sabemos que você já escolheu uma bola de ouro?
Anastasiya-Romanova 秀
7
@ Anastasiya-Romanova 秀 você poderia amostrar a partir de duas caixas boxes <- list(c(0, 1), c(1, 1))e x <- boxes[[sample(2, 1)]], em seguida , mas como esse é quase o mesmo tempo de computação, por que não usar a etapa extra que se assemelha exatamente ao processo de amostragem? Não altera nada sobre o resultado, mas torna a simulação explícita.
Tim
Tudo bem Tim, obrigado pela sua resposta. Dê-me tempo para entender sua resposta primeiro, já que sou relativamente novo em R. Por enquanto, +1 para você e para @Henry.
Anastasiya-Romanova #
1
@ Anastasiya-Romanova 秀 sim, exatamente. O código mostra uma caixa, depois uma bola da caixa, se é dourada (= 1) e depois verifica se a outra bola da caixa também é dourada (1 + 1 = 2); se sim, conta e no final, divide a soma pela contagem total (ou seja, usa mean).
Tim
1
@ Anastasiya-Romanova 秀return(NA)retorna o valor ausente e mean(, na.rm = TRUE)é usado, em que o na.rm = TRUEargumento diz à função para ignorar os valores ausentes. Em outras linguagens de programação, isso pode ser feito de maneira diferente, por exemplo, usando continueou passpalavras-chave.
Tim
4

Eu sinto que seu S/(S+0.5*(N-S))cálculo não é realmente simulação

Tente algo como isto

N <- 10^6
ballsinboxes <- c("G","G", "G","S", "S","S")
selectedbox <- sample(c(1,2,3), N, replace=TRUE)
selectedball <- sample(c(1,2), N, replace=TRUE)
selectedcolour <- ballsinboxes[(selectedbox-1)*2 + selectedball ]
othercolour <- ballsinboxes[(selectedbox-1)*2 + 3 - selectedball ]
sum(selectedcolour == "G" & othercolour == "G") / sum(selectedcolour == "G")
Henry
fonte
-2

Por que não simplesmente listar os casos?

Aqui: G é para "ouro", S é para "prata", capital é para a extração inicial:

Gg

gG

Gs

... todos os outros casos invocam uma extração inicial de prata (S) e não satisfazem a condicional (extração inicial G).

Assim, P (g | G) = 2/3.

ghuezt
fonte
7
A pergunta é sobre a solução de Monte Carlo.
Tim
Bem, listar TODAS as possibilidades é um caso extremo de Monte Carlo.
ghuezt
4
Não, não é, Monte Carlo é sobre simulação / randomização.
Tim
Tim, faça suas contas direito. Com infinitos empates, você obtém exatamente uma lista de todos os casos com probabilidades iguais. Eu triste "caso extremo" e significava limite.
ghuezt
1
Claro, mas listar todos os casos não é Monte Carlo. Quadrado é um retângulo, mas retângulo não é um quadrado.
Tim