Likelihood Poisson and prior Gamma
Marcos Lacasa
13/11/2020
Librerias
#Librerias a utilizar
library(ggplot2)## Warning: package 'ggplot2' was built under R version 4.0.4library(dplyr)## Warning: package 'dplyr' was built under R version 4.0.4##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union1 Uso de la función de verosimilitud Poisson con previa Gamma
Análisis de eventos que son independientes entre s, con probabilidad baja. Ejemplos: Recepción de WhatsApp en una hora, nº de llamadas en un da, nmero de accidentes en un determinado lugar en un perodo de tiempo …
2 Definición de Funciones
2.1 Verosimilitud - Linkelihood.
Función de probabilidad de Poisson: \(p(Y=y\mid \lambda)=\frac{e^{-\lambda} \lambda^y}{y!}\)
\(\text{E}(Y)=\lambda\)
\(\text{Var}(Y)=\lambda\)
Función de Verosimilitud \[L(\lambda,y)=\Pi\frac{e^{-\lambda} \lambda^{y_{i}}}{y_{i}!}\]
3 Ejemplo
Análisis del nmero de WhatsApp que reciben en una hora una población determinada.
OBSERVACIONES
Muestra:
#Cargamos un vector con la muestra de 25 observaciones de distintas personas con los datos de recepción de WhatsApps en 1 hora.
obs <- (c(134, 121, 150, 141, 143, 133, 145, 135, 141, 141, 137, 154, 133, 135, 154, 143, 126, 149, 143, 149, 130, 129, 142, 126, 130))
hist(obs,
breaks = 12,
main = "Histograma de las observaciones")
3.1 Función dpois
#dpois(x, # Valores del eje X (x = 0, 1, 2, ...)
# lambda, # Nmero medio de eventos que ocurren en el intervalo
# log = FALSE) # Si TRUE, las probabilidades se devuelven como logSi nos fijamos en la serie de las observaciones, vemos que el nmero minimo es 121, y el máximo 150. Tiene sentido crear un vector x, que será un vector para el cual queremos trazar nuestra curva de Funciónd de verosimilitud en Función de una lambda (nº de eventos).
rango <- seq(100, 165,by = 0.01) # Vector x, donde nos devolverá 100x65 resultados y dibujaremos la gráficaAhora creamos una Función para que nos devuelva un vector de la Función de Poisson para nuestro vector x. La Función de Verosimilitud será:
\(L(\lambda,y)=\Pi\frac{e^{-\lambda} \lambda^{y_{i}}}{y_{i}!}\)
Para cada observación \(y_{i}\) de las N observaciones. As pues vamos a calcular un vector de verosimilitud.
LikelihoodFunction <- function(y, lambda){
prod(dpois(y, lambda = lambda)) #Producto para toda las y_i
}
#Cálculo del vector para dibujar, es un vector de 6501 valores -> 165-100 * 100
LikelihoodVector <- sapply(rango, FUN = LikelihoodFunction, y=obs)#Calculamos el vector de verosimilitud con un rango de nuestro vector "x"
#Para hacer que la suma de la integral sea 1 para normalizar, hacemos:
C <- 1 / (0.01 * sum(LikelihoodVector)) # Creamos la data.frame para dibujar
LikelihoodVector <- data.frame(x = rango, y = C*LikelihoodVector, class = "versemblança")
ggplot(data = LikelihoodVector) +
geom_line(aes(x = rango, y = y, color = class, linetype = class),
size = 1.2) +
scale_linetype_manual(values=c("solid"))+
scale_color_brewer(palette="Dark2") +
ylab(expression(f(lambda))) +
xlab(expression(lambda)) +
theme_light() +
theme(axis.text = element_text(size=15),
axis.title = element_text(size=18),
legend.text = element_text(size=15)) +
theme(legend.position = "bottom", legend.title = element_blank())
# Datos Previos - Nuestra creencia
Damos por buenos que los números de WhatsApps que se reciben en 1h corresponden a una distribución de Gamma con:
nº medio de whatsApps (segn nuestra creencia) -> 128
varianza (indica incerteza del dato anterior) -> 16
Para la representarción de la Función Gamma necesitamos los valores
\(E(\lambda)=\frac{a}{b}\)
\(Var(\lambda)=\frac{a}{b^2}\)
La Función dgamma en r
#dgamma(x, #El rango -> nuestro vector "rango"
# shape, #Nuestro valor "a" -> calculamos a partir de los datos media y varianza
# rate = 1, #Nuestro valor "b" -> calculamos segn la fórmula dada.
# scale = 1/rate,
# log = FALSE) # densityCalcular gamma(a,b)
mu<- 128
sigma2 <- 16
a <- mu^2/sigma2;a## [1] 1024b <- mu/sigma2;b## [1] 8#Dibujamos la Función "a priori"
Vector_Prior <- data.frame(
x = rango , #rango
y = dgamma(rango, shape = a, rate = b), #funció a priori Gamma - amb valors a,b calculats.
class = "a priori")
data_plot <- rbind(LikelihoodVector, Vector_Prior) # Dibujamos las 2 curvas, Prior y Likelihood
ggplot(data = data_plot) +
geom_line(aes(x = x, y = y, color = class, linetype = class),
size = 1.2) +
ylab(expression(f(lambda))) +
xlab(expression(lambda)) +
theme_light() +
theme(axis.text = element_text(size=15),
axis.title = element_text(size=18),
legend.text = element_text(size=15)) +
scale_linetype_manual(values=c("solid", "twodash"))+
scale_color_brewer(palette="Dark2") +
theme(legend.position = "bottom", legend.title=element_blank())
4 Distribución a Posteriori
aP <- a + sum(obs) # la nueva "a" = a a priori + la suma de las observaciones.
bP <- b + length(obs) # la nueva "b" = b a priori + N observaciones
Vector_posterior <- data.frame(
x = rango, # rango
y = dgamma(rango, shape = aP, rate = bP ),
class = "posteriori")
data_plot <- rbind(LikelihoodVector, Vector_Prior, Vector_posterior)
ggplot(data = data_plot) +
geom_line(aes(x = x, y = y, color = class, linetype = class),
size = 1.2) +
ylab(expression(f(lambda))) +
xlab(expression(lambda)) +
theme_light() +
theme(axis.text = element_text(size=15),
axis.title = element_text(size=18),
legend.text = element_text(size=15)) +
scale_linetype_manual(values=c("solid", "twodash", "dotted"))+
scale_color_brewer(palette="Dark2") +
theme(legend.position = "bottom", legend.title=element_blank())
#media posterior
muP = aP/bP;muP## [1] 1364.1 Cálculo de Probabilidades
Probabilidades a posteriori de que la media sea superior a 135 Whatsapp
q <- 135
1-pgamma(q, aP, bP)## [1] 0.6875175