Introducción

La distribución \(\beta\) puede ser empleada para modelar la distribucion de medidas cuyos valores oscialan entre 0 y 1, y tambien para modelar la probabilidad de ocurrencia de algunos eventos discretos. En general, ha sido empleada para modelar el comportamiento de variables aleatorias acotadas a intervalos de tiempo de longitud finita en una amplia variedad de disciplinas. A parte de sus diversas aplicaciones en el campo de la simulacion estocastica, costos esperados de proyectos civiles o industriales, actualmente no se ha extendido y dispersado ampliamente al campo de la genómica (inferencia a partir de poblaciones estadisticas de un locus observable) (Johnson 2016). La distribución beta es frecuentemente usada como una distribucion a priori para proporciones de la binomial en analisis bayesiano (Forbes et al., 2011).

De otro lado la distribución teóricas binomial \(\beta inom\) fundada en una variable aleatoria Bernoulli, constituye el pilar de la experimentación y analisis probabilístico en el campo cientifico. Entre sus principales aplicaciones se incluyen, la estimacion de de probabilidades en cualquier conjunto de ensayos con exitos y fracasos, estimacion de probabilidades de resultados en los juegos de azar y muestreo de atributos (Forbes et al., 2011).

#Modelo

De acuerdo a lo anterior, suponga que se tiene una muestra aleatoria \(X_1\), \(X_2\), \(\ldots\), \(X_n\) proveniente de una distribución de probabilidad Binomial \(X~binom(n,p)\) con parámetros \(\n\) y \(\p\), donde \(\n\) numero de ensayos (no confundir con tamaño de muestra expuesto en los escenarios de simulación) y \(\p\) la probabilidad de exito.

\[P(x)=\frac{(1-x)^{\beta-1}X^{\alpha-1}}{B(\alpha, \beta)}\]

El parametro \(\p\) se distribuye segun la distribucion continua univariada \(~Beta(\alpha, \beta)\). La funcion de distribucion de una variable \(\x\) segun la distribucion apriori está dada por:

\[L(p|Datos) \propto \prod^m_{i=1} p^{Xi}(1-p)^{n-Xi}=p^{\sum\limits_{i=1}^mx_i}(1-p)^{mn-{\sum\limits_{i=1}^mx_i}}\]

La distribucion apriori: \[\epsilon(p) \propto p^{\alpha-1}(1-p)^{\beta-1} \]

La distribucion aposteriori esta dada por:

\[\epsilon(p|Datos) \propto p^{\alpha+\sum\limits_{i=1}^mx_i}(1-p)^{\beta+mn-\sum\limits_{i=1}^nx_i -1}\]

La funcion anterior no es una funcion de probabilidad, para que se cumpla este caso debe ser divididad por una constante, C, proveniente de:

\[C=\int\limits_{0}^{1} p^{\alpha+\sum\limits{xi-1}}(1-p)^{\beta+mn-\sum\limits{xi-1}}\]

Un estimador bayesiano esta dado por:

\[\epsilon(p|Datos) = \frac{\int\limits_{0}^{1} p^{\alpha+\sum\limits{xi-1}}(1-p)^{\beta+mn-\sum\limits_{i=1}^mx_i -1}}{C}\]

Desarrollo en R

Parta del supuesto donde la distribución de interés en este caso los datos provienen de una distribución binomial (X1, X2, ….X6), con parámetros \(p=1/2\) y \(n=6\).

A su vez, se piensa que p cumple o se asemeja a una distribución beta con parámetros \(\alpha=2\) y \(\beta=2\). Los escenarios a evaluar son:

# Datos Escenario I 
x1 <- c(3, 3, 4)
#Datos Escenario II 
x2 <- c(4, 3, 2, 2, 2)
#Datos Escenario III (n=100 y x=50 )
x3 <- c(1, 3, 2, 3, 3, 3, 4, 1, 3, 3)

Para cada uno de los escenarios evaluados se asumen valores constantes para los parametros de la distribucion \(\alpha=2\) y \(\beta=2\).

##Escenario 1 Para obtener el resultado mediante simulación el valor de \(C\) debe ser calculado:

#Número de simulaciones
m <- 100000
ps1 <- rbeta(n=m, shape1=2, shape2=2)
aux1 <- function(p)prod(dbinom(x=x1, size=6, prob=p))
aux1<-Vectorize(aux1)
C1 <- mean(aux1(ps1))

En este caso las curvas de la distribución apriori y aposteriori están dadas por:

posterior1 <- function(p) {
  aux1(p)* dbeta(x=p, shape1 = 2, shape2 = 2)/C1
}

prior1<-function(p){
  dbeta(x=p, shape1 = 2, shape2 = 2)
}

curve(posterior1(x), from = 0 , to=1,
      ylab="Density", main="Prior and Posterior Distributions", 
      xlab=expression(lambda), col="blue", las=1,lty=2)
curve(prior1(x), from=0, to=1,
      ylab="Density", col="darkgreen", las=1,lty=1,add=TRUE)
points(x=x1/6,y=c(0,0,0),cex=2,col="darkorchid",pch=19)
legend(5.5, 0.45, legend=c("Prior", "Posterior"),
       col=c("darkgreen", "blue"), lty=c(2,1), cex=1.2,border="white")

aux2 <- ps1*aux1(ps1)
num <- mean(aux2)
num/C1 
## [1] 0.5455923

##Escenario 2

En este caso las curvas de la distribución apriori y aposteriori están dadas por:

El estimador bayesiano del parámetro \(\p\) es:

aux2 <- ps1*aux1(ps1)
num <- mean(aux2)
num/C1 
## [1] 0.440995

##Escenario 3

Las curvas de la distribución apriori y aposteriori están dadas por:

El estimador bayesiano del parámetro \(\p\) es:

aux2 <- ps1*aux1(ps1)
num <- mean(aux2)
num/C1  
## [1] 0.4374796

#Referencias Robert W. Johnson/Applications of the Beta Distribution Part 1. Transformation group approach.

Forbes C, evans M, Hasting N , Peacock B. 2011. Statistical distributions. Fourth edition. John Wiley Published. Canada.p62