library(rmarkdown)
library(dplyr)
library(ggplot2)

1 Uso de la función de verosimilitud Binomial con una distribución previa Beta

Variables dicotómicas (valores 1,0 o Si-No). Lanzamientos de una moneda, compro o no compro, llueve o no …

1.1 Función de probabilidad de Bernouilli

\(P(Y = y | \theta)= \theta^{y}(1- \theta)^{1-y}\)

Donde \(\theta\) es la probabilidad de la variable interés, y suele puede tener un valor entre [0,1], decimos que su espacio paramétrico.

\(\Theta [0,1]\).

\(E(Y)= \theta\)

\(Var (Y)= \theta (1- \theta)\)

1.2 Prèvia conjugada

Una distribución se dice conjugada si la distribución “posterior” tiene la misma forma que la “prior”. en este caso, nuestra creencia produciría una distribución “a priori” Beta, igual que la “posterior”.

\((\theta |a,b) \sim beta(a,b)\)

Beta Distribution

\(p(\theta | a,b) = beta (\theta|a,b)\)
\(= \frac{\theta^{a-1}(1- \theta)^{b-1}}{B(a,b)}\)

Siendo \(B(a,b) =\int_0^1 \theta^{a-1}(1- \theta)^{b-1} d\theta\)

1.3 En R

\(beta (\theta|a,b)\) es la función \(dbeta(\theta, a, b)\)

\(B(a,b) = \int_0^1 \theta^{a-1}(1- \theta)^{b-1} d\theta\) es la función beta(a,b).

Vamos a ver cómo se comporta la función dbeta en función de los parámetros a,b.
Recordemos que los parámetros a y b corresponden a la distribución previa según la información o creencia que disponemos.
Vamos a suponer que tenemos información de una liga de baloncesto de un país que apenas tenemos información. Nos indican que el porcentaje de acierto de tiros libres es de un 75% (0.75). Entendemos que hay un margen de (0.15) de error.

Hacemos una representación de las formas que tiene una distribución Beta en función de sus valores (a,b):

i = 1
df_list = list()
for(a in c(1, 2, 3, 4)){ #Combinamos los valores del vector a con el vector b. 5x5 = 25 combinaciones
    for(b in c(1, 2, 3, 4)){
        x = seq(0, 1, 0.01) #El rango, que ha de ser entre 0 y 1.
        y = dbeta(x, shape1=a, shape2=b) # La función que vamos a dibujar
        df_list[[i]] = data.frame(x=x, y=y, 
                                  group = paste("b =", b, ", a =", a))
        i = i + 1
    }
}

df = bind_rows(df_list)
ggplot(df, aes(x, y)) + 
  geom_area(fill="lightblue") + 
  facet_wrap(~ group)+
  coord_cartesian(ylim=c(0,4))+
  theme_minimal()+
  theme(axis.text.x = element_text(angle = 90, hjust = 1),
        axis.text.y = element_blank())+
  xlab("Theta")+
  ylab("Probability density")

Fijémonos en la diagonal como se aproxima a una distribución normal(no puede ser normal ya que está acotada entre 0 y 1) a medida que aumentan los valores a,b proporcionalmente reduciendo su varianza.

\(\mu = a /(a+b)\) la media.

\(\omega=(a-1)/(a+b-2)\) la moda.

Pero lo interesante es averiguar las relaciones entre a,b sabiendo la media y varianza, que es lo que se suele conocer.

\(a=\mu \left(\frac{\mu(1-\mu)} {\sigma ^2}-1\right)\)

\(b=(1-\mu)\left(\frac{\mu(1-\mu)}{\sigma ^2}-1 \right)\)

Bien, ahora ya sabemos calcular tanto media y varianza, sabiendo a,b y al revés.

2 Ejemplo

Siguiendo con nuestro ejemplo, estamos haciendo un seguimiento de equipos de un país que apenas conocemos y nuestra creencia es que tienen un % de acierto de 0.75 con una desviación media de 0.15.

Varianza = a * b / {[(a + b)^2]*(a+b+1)}

media = a/c

y tenemos que:

c = media (1 - media) / sigma2 - 1 a = media * c b = (1 - media) * c

#Hallamos a,b sabiendo la media y sigma

mu <- 0.75 # seria la media del % de acierto que entendemos por nuestro conocimiento
sigma <- 0.15 # es el margen por la falta de conocimiento que tenemos. la varibilidad de la media.

c <- mu * (1 - mu) / (sigma^2) -1
a <- mu * c;a

## [1] 5.5

b <- (1 - mu)*c;b

## [1] 1.833333

# los valores a,b definen una distribución Beta para esa media y varianza


x = seq(0, 1, 0.001) # el rango, siempre entre 0,1
y = dbeta(x, shape1=a, shape2=b) #Dibujamos la función
VectorPrior = data.frame(x=x, y=y) #dibujamos el vector
mean=a/(a+b)
mode=(a-1)/(a+b-2);mode

## [1] 0.84375

ggplot(VectorPrior, aes(x, y)) + 
  geom_area(fill="lightblue") + 
  theme_minimal()+
  theme(axis.text.x = element_text(angle = 90, hjust = 1),
        axis.text.y = element_blank())+
  geom_vline(xintercept = mean, colour="red")+  #En rojo la media
  geom_vline(xintercept = mode, colour="darkblue")+#En azul la moda
  xlab("Theta")+
  ylab("Probability density")+
  ggtitle("Distribucion Prior Beta")

2.1 Likelihood

Vamos a ver en directo un partido de un equipo de esa liga en cuestión. Vamos a recoger nuestra observación de ese partido para poder evaluar, juntamente con nuestra creencia, una distribución “posterior”.

Recordemos que la función es:

\(P(Y = y | \theta)= \theta^{y}(1- \theta)^{1-y}\)

Cuando \(\theta\) debe tomar valores entre 0,1, o sea, que será el rango. Vamos a generar un ramdom con estas condiciones:

Random con distribución Binomial Pongamos por caso que ahora hacemos una observación y resulta que la estadística mejora, en 10 lanzamientos, anota 8. Donde N = 10 y a = 8

num_random <- 1
prob <- 0.8    
size <- 10
obs <-rbinom(num_random, prob = prob, size = size);obs

## [1] 7

La función de Verosimilitud es:

LikelihoodFunction <- function(prob, size, data){
  LL <- prod(dbinom(x = data, prob = prob, size = size)) 
  return(LL)
}

Ahora calculamos un vector para esta función:

rango <- seq(0, 1, 0.001)
LikelihoodVector <- sapply(rango, FUN = LikelihoodFunction, size = size, data = obs)

Dibujamos la función Likelihood, teniendo en cuenta la distribución de nuestra función de Verosimilitud con los datos recogidos:

LikelihoodVector = LikelihoodVector/sum(LikelihoodVector)  #Normalizamos para poder dibujar junto las otras distribuciones
LikelihoodVectorPlot <- data.frame(x = rango, y = LikelihoodVector, class = "Verosimilitud")
ggplot(data = LikelihoodVectorPlot) +
  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())

El cálculo de la distribución “posterior” sabiendo la distribución “prior” se realiza actualizando los valores (a,b) de la distribución beta mediante la suma de:

aP -> parámetro a, distribución posterior bP -> parámetro b, distribución posterior

Posterior Beta: Beta (a + y, b + N - y) -> siendo N (num obs), y (aciertos)

Observaciones -> 8 anotadas y 2 fallos

aP = a + 8 bP = b + 2

aP= a + 8;aP # Suma suceso éxito - anota el TL

## [1] 13.5

bP= b + 2;bP  # Auma falla el TL de los 20 que damos por bueno

## [1] 3.833333

x = rango # el rango, siempre entre 0,1
y = dbeta(x, shape1=aP, shape2=bP) #Dibujamos la función
VectorPosterior = data.frame(x=x, y=y) #dibujamos el vector

meanP=aP/(aP+bP);meanP

## [1] 0.7788462

modeP=(aP-1)/(aP+bP-2)

ggplot(VectorPosterior, aes(x, y)) + 
  geom_area(fill="lightblue") + 
  theme_minimal()+
  theme(axis.text.x = element_text(angle = 90, hjust = 1),
        axis.text.y = element_blank())+
  geom_vline(xintercept = meanP, colour="red")+  #En rojo la media
  geom_vline(xintercept = modeP, colour="darkblue")+#En azul la moda
  xlab("Theta")+
  ylab("Probability density")

3 Dibujamos las 3 juntas

#Normalizamos el vector de verosimilitud cuya integral sería 1, para compararlo con las distribuciones "prior" y "posterior"


VectorPrior <- data.frame(
        x = rango , #rango
        y = dbeta(x, shape1=a, shape2=b), #funció a priori Gamma - amb valors a,b calculats.
        class = "a priori")

LikelihoodVectorNorm<- data.frame(
        x = rango, # rango
        y = LikelihoodVector*1200, #Multiplicamos por 1200 para que sea comparable a la previa y la posterior
        class = "Verosimilitud")

VectorPosterior <- data.frame(
       x = rango , #rango
       y = dbeta(x, shape1=aP, shape2=bP), #funció a priori Gamma - amb valors a,b calculats.
       class = "posterior")


data_plot <- rbind(LikelihoodVectorNorm, VectorPrior, VectorPosterior)
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())

4 Comentario a la gráfica

A priori teníamos una creencia que se ha visto contrastada con una observación, que hace que nuestra distribución posterior recoja tanto la información previa que teníamos sobre la liga en cuestión y lo observado.

5 Cálculo de probabilidades

Qué probabilidad hay que un equipo de esa liga en cuestión tenga una media superior del 80% de tiros libres?

q = 0.8
1- pbeta(q, aP, bP)

## [1] 0.4581547

6 Qué pasaría si tuviéramos mayor confianza en los datos? -> menos varianza en nuestra creencia.

En un principio pensábamos que teníamos una creencia de mu 0.75 y sigma 0.15. Según nuestros estudios actualizados, reducimos esta varianza al 0.05. Con los mismos datos, ¿Cómo quedarían las curvas?

#Recalculamos los valores prior
mu <- 0.75 # seria la media del % de acierto que entendemos por nuestro conocimiento
sigma <- 0.05 # sigma modificada gracias a nuestro nuevo conicimiento.

c <- mu * (1 - mu) / (sigma^2) -1
a <- mu * c;a

## [1] 55.5

b <- (1 - mu)*c;b

## [1] 18.5

#Actualizamos también los parámetros a posterior

aP= a + 8;aP # Suma suceso éxito - anota el TL

## [1] 63.5

bP= b + 2;bP  # Auma falla el TL de los 20 que damos por bueno

## [1] 20.5

VectorPrior <- data.frame(
        x = rango , #rango
        y = dbeta(x, shape1=a, shape2=b), #funció a priori Gamma - amb valors a,b calculats.
        class = "a priori")

LikelihoodVectorNorm<- data.frame(
        x = rango, # rango
        y = LikelihoodVector*2200,
        class = "versemblança")

VectorPosterior <- data.frame(
       x = rango , #rango
       y = dbeta(x, shape1=aP, shape2=bP), #funció a priori Gamma - amb valors a,b calculats.
       class = "posterior")


data_plot <- rbind(LikelihoodVectorNorm, VectorPrior, VectorPosterior)
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())

Al estar más seguros en un inicio, los valores de la observación no han tenido tanta influencia en nuestra distribución “posterior”. Comprobamos la probabilidades ante la misma pregunta: