Esta é a pergunta de um iniciante em um exercício na “Computação Bayesiana com R” de Jim Albert. Observe que, embora isso possa ser um dever de casa, no meu caso não é, pois estou aprendendo métodos bayesianos em R porque acho que posso usá-lo em minhas análises futuras.
De qualquer forma, embora seja uma pergunta específica, provavelmente envolve uma compreensão básica dos métodos bayesianos.
Portanto, no exercício 2.2, Jim Albert nos pede para analisar o experimento de um centavo. Veja aqui. Devemos usar um histograma antes, ou seja, dividir o espaço de pvalores possíveis em 10 intervalos de comprimento .1e atribuir uma probabilidade prévia a eles.
Como sei que a verdadeira probabilidade será .5e acho altamente improvável que o universo tenha mudado as leis da probabilidade ou que o centavo seja robusto, meus priores são:
prior <- c(1,5,20,100,5000,5000,100,20,5,1)
prior <- prior/sum(prior)
ao longo dos pontos médios do intervalo
midpt <- seq(0.05, 0.95, by=0.1)Por enquanto, tudo bem. Em seguida, giramos o centavo 20 vezes e registramos o número de sucessos (cabeças) e falhas (cauda). Feito facilmente:
y <- rbinom(n=20,p=.5,size=1)
s <- sum(y==1)
f <- sum(y==0)
Na minha experiência, s == 7e f == 13. Em seguida, vem a parte que eu não entendo:
Simule a partir da distribuição posterior calculando (1) a densidade posterior de p em uma grade de valores em (0,1) e (2) colhendo uma amostra simulada com substituição da grade. (A função
histprioresamplesão úteis neste cálculo). Como as probabilidades de intervalo foram alteradas com base nos seus dados?
É assim que se faz:
p <- seq(0,1, length=500)
post <- histprior(p,midpt,prior) * dbeta(p,s+1,f+1)
post <- post/sum(post)
ps <- sample(p, replace=TRUE, prob = post)
Mas por que fazemos isso ?
Podemos facilmente obter a densidade posterior multiplicando o anterior pela probabilidade apropriada, conforme feito na linha dois do bloco acima. Este é um gráfico da distribuição posterior:
Como a distribuição posterior é ordenada, podemos obter resultados para os intervalos introduzidos no histograma anteriormente, resumindo os elementos da densidade posterior:
post.vector <- vector()
post.vector[1] <- sum(post[p < 0.1])
post.vector[2] <- sum(post[p > 0.1 & p <= 0.2])
post.vector[3] <- sum(post[p > 0.2 & p <= 0.3])
post.vector[4] <- sum(post[p > 0.3 & p <= 0.4])
post.vector[5] <- sum(post[p > 0.4 & p <= 0.5])
post.vector[6] <- sum(post[p > 0.5 & p <= 0.6])
post.vector[7] <- sum(post[p > 0.6 & p <= 0.7])
post.vector[8] <- sum(post[p > 0.7 & p <= 0.8])
post.vector[9] <- sum(post[p > 0.8 & p <= 0.9])
post.vector[10] <- sum(post[p > 0.9 & p <= 1])
(Os especialistas em R podem encontrar uma maneira melhor de criar esse vetor. Acho que pode ter algo a ver com isso sweep?)
round(cbind(midpt,prior,post.vector),3)
midpt prior post.vector
[1,] 0.05 0.000 0.000
[2,] 0.15 0.000 0.000
[3,] 0.25 0.002 0.003
[4,] 0.35 0.010 0.022
[5,] 0.45 0.488 0.737
[6,] 0.55 0.488 0.238
[7,] 0.65 0.010 0.001
[8,] 0.75 0.002 0.000
[9,] 0.85 0.000 0.000
[10,] 0.95 0.000 0.000
Além disso, temos 500 empates da distribuição posterior que não nos dizem nada de diferente. Aqui está um gráfico da densidade dos desenhos simulados:
Agora usamos os dados simulados para obter probabilidades para nossos intervalos, contando qual a proporção de simulações dentro do intervalo:
sim.vector <- vector()
sim.vector[1] <- length(ps[ps < 0.1])/length(ps)
sim.vector[2] <- length(ps[ps > 0.1 & ps <= 0.2])/length(ps)
sim.vector[3] <- length(ps[ps > 0.2 & ps <= 0.3])/length(ps)
sim.vector[4] <- length(ps[ps > 0.3 & ps <= 0.4])/length(ps)
sim.vector[5] <- length(ps[ps > 0.4 & ps <= 0.5])/length(ps)
sim.vector[6] <- length(ps[ps > 0.5 & ps <= 0.6])/length(ps)
sim.vector[7] <- length(ps[ps > 0.6 & ps <= 0.7])/length(ps)
sim.vector[8] <- length(ps[ps > 0.7 & ps <= 0.8])/length(ps)
sim.vector[9] <- length(ps[ps > 0.8 & ps <= 0.9])/length(ps)
sim.vector[10] <- length(ps[ps > 0.9 & ps <= 1])/length(ps)
(Mais uma vez: existe uma maneira mais eficiente de fazer isso?)
Resuma os resultados:
round(cbind(midpt,prior,post.vector,sim.vector),3)
midpt prior post.vector sim.vector
[1,] 0.05 0.000 0.000 0.000
[2,] 0.15 0.000 0.000 0.000
[3,] 0.25 0.002 0.003 0.000
[4,] 0.35 0.010 0.022 0.026
[5,] 0.45 0.488 0.737 0.738
[6,] 0.55 0.488 0.238 0.236
[7,] 0.65 0.010 0.001 0.000
[8,] 0.75 0.002 0.000 0.000
[9,] 0.85 0.000 0.000 0.000
[10,] 0.95 0.000 0.000 0.000
Não é de surpreender que a simulação não produza outros resultados além da posterior, na qual se baseou. Assim, por que desenhamos essas simulações em primeiro lugar ?
fonte
ps <- sample(p, replace=TRUE, prob = post)! Estou certo ao supor que isso mudará para técnicas de simulação mais avançadas?Respostas:
Para responder à sua subquestão: Como fazer o seguinte de maneira mais elegante?
A maneira mais fácil de fazer isso usando a base R é:
Observe que as quebras variam de 0 a 1. Isso gera:
E nós temos:
fonte
Meu entendimento é que, uma vez que a densidade posterior obtida a partir do produto de densidade e probabilidade anteriores é apenas uma APROXIMA da densidade posterior, portanto, não podemos fazer nenhuma inferência EXATA diretamente.
Consequentemente, precisamos coletar uma amostra aleatória e realizar inferência a partir da amostra, assim como o método de simulação para a família beta posterior.
fonte