Likelihood exponential prior Gamma
Marcos Lacasa
13/11/2020
library(dplyr)
library(ggplot2)1 Us Likelihood Exponential - prior Gamma
Se utiliza para modelizar el tiempo que transcurre entre 2 eventos.
2 Exemple
Vamos a analizar la variable “tiempo en minutos que esperamos” antes que nos atiendan. Supongamos que tenemos una creencia, basada en estudios previos, que la espera media es de 10 minutos, pero no nos fiamos mucho, y aceptamos un margen de error de 10 (o sea que estaría entre 0-20). Indica poca fiabilidad en estos estudios.
3 Distribució Gamma
La distribución conjugada para este tipo de análisis es la distribución gamma, y la función de verosimilitud se ajusta a la distribución exponencial.
\(f_\Theta(\theta) = \frac{\beta^\alpha\theta^{\alpha - 1}e^{-\beta\theta}}{\Gamma(\alpha)}\)
4 Likelihood - Exponential
Funció de densitat:
\(f_{X_i|\Theta}(x_i|\theta) = \theta e^{-\theta x_i}\)
Likelihood Function:
\(f(X|\theta) = \theta^n e^{-\theta \sum x_i}\)
on:
N -> número d’observacions
\(\sum x_i =\) suma d’observacions
5 Posterior
\(\text{Posterior} \propto \text{Likelihood} \times \text{Prior}\)
\(\Rightarrow f_{\Theta | X}(\theta | x) \propto f_{x|\Theta}(x|\theta)\times f_\Theta(\theta)\)
\(f_{\Theta | X}(\theta | x) \propto \theta^n e^{-\theta \sum x_i} \times \frac{\beta^\alpha\theta^{\alpha - 1}e^{-\beta\theta}}{\Gamma(\alpha)}\)
6 Implementamos en R
Creamos la función de Verosimilitud donde:
lambda -> rango
n_obs -> nº de observaciones
sum_obs -> suma de todas las observaciones
Las observaciones las vamos a crear mediante una función aleatoria utilizando la función rgamma, con los siguientes datos:
media = 10 (minuts) -> \(\mu\) (de la observación)
desviación tipo \(\sigma\) = 10
De la distribución exponencial sabemos que:
\(\theta = 1 / \mu\)
\(\theta^{2} = 1 / \sigma^{2}\)
# Para simular una serie exponencial sólo necesitamos saber el nº obs y valor de lambda.
set.seed(10)
mu<- 10
lambda <- 1/mu
obs <- rexp(10, lambda)
N <- 10
hist(obs,
breaks = 5,
main = "")
obs## [1] 0.1495641 9.2022120 7.5215894 15.7504185 2.3165862 10.8667300
## [7] 23.2762287 7.2912382 12.8831010 6.7226829- Sobre ese lambda, parámetro desconocido no observable, es sobre el que hacemos inferencia y empezamos por una distribución a priori gamma.
Si pensamos, a priori, que lambda tiene un valor de 0.1 (lo que supondría una media de 10 minutos de espera) podemos establecer un rango entre 0.001 (media de 1000) y 0.25 (media de 4)
rango <- seq(0.001,0.25, length.out = 10000)Como varianza podemos pensar en una varianza de 1 que ya sería grande en este caso.
Para la distribución gamma tenemos que calcular los parámetros (a,b) que tienen la siguiente relación con la variable a estudio \(\theta\).
\(E(\theta) = a/b\)
\(Var(\theta) = a/b^{2}\)
Pero recordemos las relaciones entre \(\theta\) y los datos iniciales.
Dato inicial, media de 10 minutos y margen (desviación tipo) de 10:
\(\mu(\theta) = 1/10 = 0.1\)
\(\sigma(\theta) = 1/100 = 0.01\)
#Dibujamos la función "a priori" amb els valors a,b de la distribució gamma
mu <- 0.1
sigma2 <-0.01
#sigma2 <- sigma^2
a <- mu^2/sigma2;a## [1] 1b <- mu/sigma2;b## [1] 10Vector_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")
ggplot(data = Vector_Prior) +
geom_line(aes(x = x, y = y, color = class, linetype = class),
size = 1.2) +
ylab(expression(f(theta))) +
xlab(expression(theta)) +
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())
Ahora calculamos/dibujamos la verosimilitud. Vemos que esta centrada en la media de las observaciones. Por otra parte, podemos fijarnos en que los valores del eje de la y son muy pequeños incluso si intentamos normalizarlos por lo que para pintarlo de forma que sea comparable con la posterior tendremos que multiplicar por alguna constante. Pero ojo, esto es solo a efectos de pintarlo, en el fondo, en nuestro análisis sobre lambda no nos hace falta.
LikelihoodFunction <- function(obs, theta){
(theta^N)*exp(-theta*obs)
}
#rango theta
LikelihoodFunction(obs, 0.1)## [1] 9.851549e-11 3.984309e-11 4.713478e-11 2.069989e-11 7.932169e-11
## [6] 3.373369e-11 9.752731e-12 4.823314e-11 2.757364e-11 5.105492e-11length(rango)## [1] 10000#Likelihood Vector
LikelihoodVector <- sapply(rango, obs = obs, FUN = LikelihoodFunction)sum(LikelihoodVector)## [1] 0.002123164LikelihoodFunction <- function(obs, theta){
(theta^N)*exp(-theta*sum(obs))
}
LikelihoodVector <- sapply(rango, obs = obs, FUN = LikelihoodFunction)
data <- data.frame(x = rango, y = LikelihoodVector/sum(LikelihoodVector), class = "versemblança")
ggplot(data = data) +
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(theta))) +
xlab(expression(theta)) +
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())
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 según la fórmula dada.
# scale = 1/rate,
# log = FALSE) # density7 Distribución a Posteriori
aP <- a + length(obs);aP # la nueva "a" = a a priori + Numero de observaciones.## [1] 11bP <- b + sum(obs);bP # la nueva "b" = b a priori + la suma de las observaciones## [1] 105.9804Vector_posterior <- data.frame(
x = rango, # rango
y = dgamma(rango, shape = aP, rate = bP ),
class = "posteriori")
LikelihoodVector_norm<- data.frame(
x = rango, # rango
y = 30000*LikelihoodVector/sum(LikelihoodVector), #Multiplicamos por 30000 para que sea comparable a la previa y la posterior
class = "versemblança")
data_plot <- rbind(LikelihoodVector_norm, 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(theta))) +
xlab(expression(theta)) +
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())
La media de la posterior
muP <- aP/bP;muP## [1] 0.10379287.1 Cálculo de Probabilidades
Probabilidades a posteriori de que la media sea menos de 10 min -> lambda = 0.1
q <- 0.1
1-pgamma(q, aP, bP)## [1] 0.50865048 y qué pasaría si la varianza fuera mucho más baja?
Partimos de la base que esos 10 minutos que teníamos de creencia con un margen de error de 10 minutos, nos han llegado estudios fiables que hacen que ese margen de eror se reduzca a la mitad!!
mu<- 1/10
sigma2 <- 1/25
a <- mu^2/sigma2;a## [1] 0.25b <- mu/sigma2;b## [1] 2.5#simulem les nostres observacions
set.seed(10)
mu<- 10
theta <- 1/mu
obs <- rexp(10, theta)
N <- 10
hist(obs,
breaks = 5,
main = "")
obs## [1] 0.1495641 9.2022120 7.5215894 15.7504185 2.3165862 10.8667300
## [7] 23.2762287 7.2912382 12.8831010 6.7226829#Ampliem el rang segons la nova distribució de les observacions
rango <- seq(0.001,0.35, length.out = 10000)
#Likelihood Vector
LikelihoodVector <- sapply(rango, obs = obs, FUN = LikelihoodFunction)LikelihoodFunction <- function(obs, theta){
(theta^N)*exp(-theta*sum(obs))
}
LikelihoodVector <- sapply(rango, obs = obs, FUN = LikelihoodFunction)
data <- data.frame(x = rango, y = LikelihoodVector/sum(LikelihoodVector), class = "versemblança")
ggplot(data = data) +
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(theta))) +
xlab(expression(theta)) +
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())
aP <- a + length(obs);aP # la nueva "a" = a a priori + Numero de observaciones.## [1] 10.25bP <- b + sum(obs);bP # la nueva "b" = b a priori + la suma de las observaciones## [1] 98.48035Vector_posterior <- data.frame(
x = rango, # rango
y = dgamma(rango, shape = aP, rate = bP ),
class = "posteriori")
LikelihoodVector_norm<- data.frame(
x = rango, # rango
y = 30000*LikelihoodVector/sum(LikelihoodVector), #Multiplicamos por 30000 para que sea comparable a la previa y la posterior
class = "versemblança")
data_plot <- rbind(LikelihoodVector_norm, 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(theta))) +
xlab(expression(theta)) +
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())
Probabilidades a posteriori de que la media sea menos de 10 min
q <- 0.1
1-pgamma(q, aP, bP)## [1] 0.5089927