Este es un principiante pregunta en un ejercicio de Jim Albert "Bayesiano de Cálculo con R". Tenga en cuenta que mientras que esto puede ser una tarea, en mi caso no es, como yo soy el aprendizaje Bayesiano métodos en R porque creo que podría utilizar en mi futuro análisis.
De todos modos, si bien esta es una pregunta específica, probablemente implica la comprensión básica de los métodos Bayesianos.
Así, en el ejercicio 2.2, Jim Albert nos pide analizar el experimento de un centavo de tiro. Vea aquí. Vamos a utilizar un histograma antes, es decir, dividir el espacio de lo posible p valores en 10 intervalos de longitud .1 y asignar una probabilidad anterior a estos.
Ya sé que la verdadera probabilidad será .5, y creo que es muy poco probable que el universo ha cambiado las leyes de la probabilidad o de la moneda es robusto, mis dudas son:
prior <- c(1,5,20,100,5000,5000,100,20,5,1)
prior <- prior/sum(prior)
a lo largo del intervalo de los puntos medios
midpt <- seq(0.05, 0.95, by=0.1)
Hasta ahora tan bueno. A continuación, nos girar la moneda 20 veces y registrar el número de éxitos (los jefes) y fracasos (la cola). Fácil de hacer:
y <- rbinom(n=20,p=.5,size=1)
s <- sum(y==1)
f <- sum(y==0)
En mi expeniment, s == 7 y f == 13. Luego viene la parte que no entiendo:
Simular a partir de la distribución posterior por (1) el cómputo de la posterior de la densidad de p en una cuadrícula de valores en (0,1) y (2) tomar una muestras simuladas con la sustitución de la red. (La función
histpriorysampleson útiles en este cálculo). ¿Cómo han el intervalo de probabilidades cambiado en la base de datos?
Esto es cómo se hace:
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)
Pero ¿por qué lo hacemos?
Podemos obtener fácilmente por la parte posterior de la densidad multiplicando el anterior con la correspondiente probabilidad, como se ha hecho en la línea dos del bloque anterior. Este es un gráfico de la distribución posterior:
Como la distribución posterior es ordenado, podemos obtener resultados para los intervalos introducido en el histograma antes con un resumen sobre los elementos de la parte posterior de la densidad:
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])
(R expertos podría encontrar una mejor manera de crear ese vector. Supongo que podría tener algo que ver con 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
Además, hemos 500 dibuja a partir de la distribución posterior que nos dicen nada diferente. Aquí es un gráfico de la densidad de la simulación de los sorteos:
Ahora vamos a utilizar los datos simulados para obtener las probabilidades para nuestro intervalos contando qué proporción de las simulaciones están dentro del 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)
(De nuevo: ¿hay una manera más eficiente de hacer esto?)
Resumen de 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
Se trata no es de sorprender que la simultion produce ningún otro resultado que el posterior, en el que se basa. Por lo tanto, ¿por qué hemos de extraer las simulaciones en el primer lugar?
