SR2 Capítulo 3 Médio | R-bloggers

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


[Esteartigofoipublicadopelaprimeiravezem[Thisarticlewasfirstpublishedon Brian Callander, 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.

Aqui está minha solução para os exercícios médios no capítulo 3 do Statistical Reethinking da McElreath, 2ª edição.

( DeclareMathOperator { dbinomial} {Binomial} DeclareMathOperator { dbernoulli} {Bernoulli} DeclareMathOperator { dpoisson} {Poisson} DeclareMathOperator { dnormal} {Normal} DeclareMathOperator { dnormal} {Normal} DeclareMathOperator { dnormalom { dcauchy} {Cauchy} DeclareMathOperator { dexponential} {Exp} DeclareMathOperator { duniform} {Uniform} DeclareMathOperator { dgamma} {Gamma} DeclareMathOperator { dinvpamma} {Invpathma} {InvpathmaO} Logit } DeclareMathOperator { logit} {Logit} DeclareMathOperator { ddirichlet} {Dirichlet} DeclareMathOperator { dbeta} {Beta} )

Supondo que a Terra tenha 70% de cobertura de água e observamos a água 8 vezes em 15 lançamentos do globo, vamos calcular algumas quantidades posteriores com duas opções de priorização: uniforme e passo.

p_true  0.7

W  8
N  15

granularity  1000 # points on the grid

Uniforme anterior

Calculamos a aproximação da grade posterior, como mostrado no livro.

m1_grid  tibble(p = seq(, 1, length.out = granularity)) %>% 
  mutate(prior = 1)

m1_posterior  m1_grid %>% 
  mutate(
    likelihood = dbinom(W, N, p),
    posterior = prior * likelihood
  )
Solução para exercitar 3M1
Solução para exercitar 3M1

Podemos obter empates do nosso posterior, amostrando os valores de cobertura da água muitas vezes com a substituição, cada valor sendo desenhado proporcionalmente à probabilidade posterior. Podemos então resumir esses empates para obter o intervalo desejado.

m2_samples  m1_posterior %>% 
  sample_n(10000, replace = T, weight = posterior)

m2_hpdi  HPDI(m2_samples$p, prob = 0.9)
m2_hpdi
     |0.9      0.9| 
0.3223223 0.7097097 

O histograma tem a seguinte aparência. É o mesmo que o gráfico anterior, mas calculado a partir das amostras.

Solução para o exercício 3M2
Solução para exercitar 3M2

Para obter a amostra preditiva posterior, fazemos nossos sorteios posteriores (p ), use-os para desenhar um número aleatório de lançamentos de água observados em 15. A fração de amostras preditivas posteriores com um determinado valor é a probabilidade preditiva posterior desse valor.

m3_prob  m2_samples %>% 
  mutate(W = rbinom(n(), 15, p)) %>% 
  group_by(W) %>% 
  tally() %>% 
  mutate(probability = n / sum(n))
Solução para o exercício 3M3
Solução para exercitar 3M3

Também podemos calcular as probabilidades preditivas posteriores com um número diferente de lançamentos. Aqui com 9 lançamentos.

m4_prob  m2_samples %>% 
  mutate(W = rbinom(n(), 9, p)) %>% 
  group_by(W) %>% 
  tally() %>% 
  mutate(probability = n / sum(n))
Solução para o exercício 3M4
Solução para exercitar 3M4

Etapa anterior

Agora repetimos os mesmos passos, mas com o passo anterior em vez do uniforme anterior. Vamos repetir sem comentar.

Leia Também  Como verificar sua necessidade de seguro?
cupom com desconto - o melhor site de cupom de desconto cupomcomdesconto.com.br
m5_grid  m1_grid %>% 
  mutate(prior = if_else(p  0.5, , 1))

m5_posterior  m5_grid %>% 
  mutate(
    likelihood = dbinom(W, N, p),
    posterior = prior * likelihood
  )
Solução para o exercício 3M5 parte 1
Solução para o exercício 3M5 parte 1
m5_samples  m5_posterior %>% 
  sample_n(10000, replace = T, weight = posterior)

m5_hpdi  HPDI(m5_samples$p, prob = 0.9)
m5_hpdi
     |0.9      0.9| 
0.5005005 0.7107107 
Solução para o exercício 3M5 parte 2
Solução para o exercício 3M5 parte 2
m5_prob  m5_samples %>% 
  mutate(W = rbinom(n(), 15, p)) %>% 
  group_by(W) %>% 
  tally() %>% 
  mutate(probability = n / sum(n))
Solução para o exercício 3M5 parte 3
Solução para o exercício 3M5 parte 3
m5_prob  m5_samples %>% 
  mutate(W = rbinom(n(), 9, p)) %>% 
  group_by(W) %>% 
  tally() %>% 
  mutate(probability = n / sum(n))
Solução para o exercício 3M5 parte 4
Solução para o exercício 3M5 parte 4

Vamos comparar a proporção de amostras dentro de 0,05 do valor verdadeiro de cada anterior.

p_close_uniform  m2_samples %>% 
  group_by(close = p %>% between(p_true - 0.05, p_true + 0.05)) %>% 
  tally() %>% 
  mutate(probability = n / sum(n)) %>% 
  filter(close) %>% 
  pull(probability)

p_close_step  m5_samples %>% 
  group_by(close = p %>% between(p_true - 0.05, p_true + 0.05)) %>% 
  tally() %>% 
  mutate(probability = n / sum(n)) %>% 
  filter(close) %>% 
  pull(probability)

A probabilidade de estar próximo do valor real sob o uniforme e as etapas anteriores é de 0,1316 e 0,2157, respectivamente. O passo anterior, portanto, tem mais massa em torno do valor verdadeiro.

Exercício 3M6

Os modelos bayesianos são generativos, o que significa que podemos simular novos conjuntos de dados de acordo com nossas probabilidades anteriores. Simularemos 10 conjuntos de dados para cada valor de N de interesse. Simulamos um conjunto de dados escolhendo aleatoriamente um p_true do nosso uniforme anterior, e então escolher aleatoriamente W da distribuição binomial correspondente.

m6_prior_predictive  crossing(
    N = 200 * (1:16), 
    iter = 1:10
  ) %>% 
  mutate(
    p_true = runif(n(), min=, max=1), 
    W = rbinom(n(), N, p_true)
  )

Para cada um desses conjuntos de dados simulados, fazemos uma grade aproximada da posterior, coletamos amostras posteriores e calculamos o HPDI.

m6_grid  tibble(p = seq(, 1, length.out = granularity)) %>% 
  mutate(prior = 1)

m6_posteriors  m6_prior_predictive %>% 
  crossing(m6_grid) %>% 
  group_by(N, p_true, iter) %>% 
  mutate(
    likelihood = dbinom(W, N, p),
    posterior = prior * likelihood
  )

m6_samples  m6_posteriors %>% 
  sample_n(1000, replace = TRUE, weight = posterior) 

m6_hpdi  m6_samples %>% 
  summarise(lo = HPDI(p, 0.99)[1], hi = HPDI(p, 0.99)[2]) %>% 
  mutate(width = abs(hi - lo))

Agora, para cada valor de N, verificamos quantos dos intervalos têm a largura desejada.

m6_n  m6_hpdi %>% 
  group_by(N) %>% 
  summarise(fraction = mean(width  0.05)) 
Solução para o exercício 3M6
Solução para o exercício 3M6

Portanto, esperamos que um tamanho de amostra em torno de 2600-3000 nos forneça uma estimativa posterior suficientemente precisa.

Leia Também  Recriar um aplicativo brilhante com o Flask



Se você chegou até aqui, por que não inscreva-se para atualizações do site? Escolha seu sabor: e-mail, Twitter, RSS ou facebook …



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