Monty hall e diagramas de influência

Redes bayesianas

As redes Bayesianas são o resultado da combinação de conceitos probabilísticos e conceitos da teoria dos grafos. Segundo Pearl, tal união tem como consequências três benefícios: i) prover formas convenientes para expressar suposições do modelo; ii) facilitar a representação de funções de probabilidade conjuntas; e iii) facilitar o cálculo eficiente de inferências a partir de observações.

Da teoria de probabilidades precisamos apenas de alguns resultados básicos sobre probabilidade condicional. Primeiramente, pela definição de probabilidade condicional, sabemos que

\[ p(x_1, x_2) = p(x_1)p(x_2|x_1). \]

Aplicando essa regra iterativamente para \(n\) variáveis, temos

\[ p(x_1, \dots, x_p) = \prod_j p(x_j|x_1,\dots, x_{j-1}). \]

Agora, imagine que, no seu problema, a variável aleatória \(X_j\) não dependa probabilisticamente de todas as variáveis \(X_1,\dots, X_{j-1}\), e sim apenas de um subconjunto \(\Pi_j\) dessas variáveis. Fazendo isso, a equação pode ser escrita como

\[ p(x_1, \dots, x_p) = \prod_j p(x_j|\pi_j). \]

Chamamos \(\Pi_j\) de pais de \(X_j\). Esse conjunto pode ser pensado como as variáveis que são suficientes para determinar as probabilidades de \(X_j\).

A parte mais legal das redes Bayesianas é que elas podem ser representadas a partir de DAGs (grafos direcionados acíclicos). No grafo, se \(X_1\) aponta para \(X_2\), então \(X_1\) é pai de \(X_2\). Por exemplo, esse grafo aqui

representa a distribuição de probabilidades \(p(x_1, \dots, x_5)\) com

\[ p(x_1, \dots, x_5) = p(x_1)p(x_2|x_1)p(x_3|x_1)p(x_4|x_3,x_2)p(x_5|x_4). \]

Diagrama de influências

Um diagrama e influências é uma rede Bayesiana com nós de decisão e utilidade (ganhos). Ou seja, é uma junção de três conceitos:

\[ \underbrace{\text{prob. condicional} + \text{grafos}}_{\text{rede Bayesiana}} + \text{teoria da decisão} = \text{diagrama de influências} \]

Na teoria da decisão, usualmente estamos interessados em maximizar a utilidade esperada. No diagrama, considerando a estrutura de probabilidades dada pela rede Bayesiana e as informações disponíveis, queremos escolher a decisão que faz com que, em média, nosso retorno seja mais alto.

Com diagramas de influências, é possível organizar sistemas complexos com múltiplas decisões, considerando diferentes conjuntos de informações disponíveis. É uma ferramenta realmente muito poderosa.

Voltando ao Monty Hall

Agora que sabemos um pouquinho de diagramas de influência, podemos desenhar o do Monty Hall:

O jogador tem duas decisões a tomar:

• \(D_1\) (escolha_inicial): A escolha da porta inicial (1, 2, 3).
• \(D_2\) (trocar): Trocar a porta ou não (s, n).

Também temos duas fontes de incerteza:

• \(X_1\) (ferrari): Em qual porta está a Ferrari (1, 2, 3).
• \(X_2\) (porta_retirada): Qual porta foi retirada (1, 2, 3). Essa variável não é sempre aleatória: se eu escolho a porta 1 e a Ferrari está em 2, o apresentador é obrigado a retirar a porta 3. Se o apresentador tiver a opção de escolher (que acontece no caso da escolha inicial ser a Ferrari), o apresentador escolhe uma porta para retirar aleatoriamente.

Finalmente, temos um nó de utilidade:

• \(U_1\) (result): Ganhei a Ferrari (ganhei, perdi).

Em R, podemos construir a rede Bayesiana do problema utilizando o pacote bnlearn:

# nós do grafo
nodes <- c("escolha_inicial", "ferrari", "porta_retirada", "trocar", "result")

# matriz de adjacências
edges <- matrix(
  c("escolha_inicial", "porta_retirada",
    "ferrari", "porta_retirada",
    "porta_retirada", "trocar",
    "trocar", "result",
    "ferrari", "result",
    "escolha_inicial", "result"),
  ncol = 2, 
  byrow = TRUE)

edges

##      [,1]              [,2]            
## [1,] "escolha_inicial" "porta_retirada"
## [2,] "ferrari"         "porta_retirada"
## [3,] "porta_retirada"  "trocar"        
## [4,] "trocar"          "result"        
## [5,] "ferrari"         "result"        
## [6,] "escolha_inicial" "result"

# criando o grafo a partir de um grafo vazio
g <- bnlearn::empty.graph(nodes)
bnlearn::arcs(g) <- edges

O output desse conjunto de operações é um objeto do tipo bn com várias propriedades pré calculadas pelo pacote bnlearn:

g

  Random/Generated Bayesian network

  model:
   [escolha_inicial][ferrari][porta_retirada|escolha_inicial:ferrari][trocar|porta_retirada]
   [result|escolha_inicial:ferrari:trocar]
  nodes:                                 5 
  arcs:                                  6 
    undirected arcs:                     0 
    directed arcs:                       6 
  average markov blanket size:           3.60 
  average neighbourhood size:            2.40 
  average branching factor:              1.20 

  generation algorithm:                  Empty 

Com as especificação do problema dada, se gerarmos aleatoriamente todos os cenários, chegamos à essa combinação de casos equiprováveis (ver Extra 2)

Agora, vamos escrever todas as combinações possíveis de cenários e guardar num data.frame chamado dados:

Finalmente, ajustamos nossa rede Bayesiana, usando a função bnlearn::bn.fit().

fit <- bnlearn::bn.fit(g, dados)

A função bnlearn::cpquery() (conditional probability query) serve para realizar uma consulta de probabilidades dada a rede ajustada. No nosso caso, a partir de uma escolha inicial qualquer \(d_1\), queremos saber o ganho ao trocar é maior que o ganho ao não trocar.

\[ \mathbb E(U_1\; |\; D_2 = \text{s}, D_1 = d_1) > \mathbb E(U_1\; |\; D_2 = \text{n}, D_1 = d_1). \]

Fazendo contas, isso equivale matematicamente a consultar se

\[ \mathbb P(U_1=\text{ganhei}\; |\; D_2 = \text{s}) > \mathbb P(U_1=\text{ganhei}\; |\; D_2 = \text{n}) \]

Agora, podemos consultar \(\mathbb P(U_1=\text{ganhei}\; |\; D_2 = \text{s})\) com nosso modelo!

set.seed(13)                    # reprodutibilidade
bnlearn::cpquery(
  fitted = fit, 
  event = (result == "ganhei"), # o que queremos saber?
  evidence = (trocar == "s"),   # qual informação adicionar?
  n = 5e6)                      # n grande para aumentar a precisão

[1] 0.6666704

E não é que dá 2/3 mesmo? Da mesma forma, temos

bnlearn::cpquery(fit, (result == "ganhei"), (trocar == "n"), n = 5e6)

[1] 0.3333187

Resolvido!

Wrap-up

• Vale à pena trocar a porta!
• Redes Bayesianas juntam grafos e probabilidades condicionais
• Diagramas de influência juntam redes Bayesianas e teoria da decisão
• Essas ferramentas podem ser utilizadas tanto para resolver Monty Hall quanto para ajudar em sistemas complexos.

É isso pessoal. Happy coding ;)

diagrama <- "
graph LR;
A{D1Escolha inicial}-->B(X2Porta retirada);
C(X1Ferrari)-->B;
B-->D{D2Trocar porta};
D-->E[U1Ganhar]
C-->E
A-->E
"
# tweak para centralizar e grifar as variáveis
diagrama <- stringr::str_replace_all(
  diagrama, 
  pattern = "([XDU][0-9])", 
  replacement = "<center><b>\\1</b></center><br/>")

DiagrammeR::DiagrammeR(diagrama)

set.seed(13)

simular_monty_hall <- function(z = 0) {
  v <- 1:3                                  # opcoes
  escolha_inicial <- sample(v, 1)           # escolha inicial aleatoria
  ferrari <- sample(v, 1)                   # ferrari aleatoria
  
  # qual porta retirar?
  if (escolha_inicial == ferrari) {
    porta_retirada <- sample(setdiff(v, ferrari), 1)
  } else {
    porta_retirada <- setdiff(v, c(escolha_inicial, ferrari))
  }
  
  # trocar porta?
  trocar <- sample(c("s", "n"), 1)
  
  # calculando resultado
  if (trocar == "s") {
    escolha_final <- setdiff(v, c(escolha_inicial, porta_retirada))
  } else {
    escolha_final <- escolha_inicial
  }
  result <- ifelse(escolha_final == ferrari, "ganhei", "perdi")
  
  # guardando no BD
  tibble::tibble(escolha_inicial, ferrari, porta_retirada, trocar, result)
}

dados_simulados <- purrr::map_dfr(seq_len(1e4), simular_monty_hall) %>% 
  dplyr::mutate_all(as.factor)

dplyr::glimpse(dados_simulados)

Observations: 10,000
Variables: 5
$ escolha_inicial <fct> 3, 1, 2, 1, 1, 1, 3, 1, 2, 3, 3, 1, 3, 1, 2, 2, 2,...
$ ferrari         <fct> 1, 1, 2, 1, 1, 2, 3, 3, 1, 2, 3, 3, 2, 1, 1, 3, 1,...
$ porta_retirada  <fct> 2, 3, 1, 3, 2, 3, 2, 2, 3, 1, 1, 2, 1, 2, 3, 1, 3,...
$ trocar          <fct> n, s, s, n, s, n, n, n, n, s, s, s, s, n, n, s, n,...
$ result          <fct> perdi, perdi, perdi, ganhei, perdi, perdi, ganhei,...

dados_simulados %>% 
  dplyr::distinct() %>% 
  dplyr::arrange(escolha_inicial, ferrari)