Tema 2: Modelos Básicos de

Estadı́stica Bayesiana

Introducción
Desde el punto de vista de la estadı́stica bayesiana se considera que hay dos tipos de
valores: conocidos y desconocidos. El objetivo es usar las cantidades o datos conocidos
mediante un modelo paramétrico dado para realizar inferencias sobre las cantidades des-
conocidas. Por cantidades desconocidas se puede entender tanto parámetros del modelo,
como observaciones missing o predicciones.
En un modelo básico se tiene un parámetro de interés θ y unos datos observados D
y se considera una distribución de probabilidad conjunta para ambos que recoge cómo se
relacionan: p(θ, D).
Aplicando la definición de probabilidad condicionada, se tiene que

p(θ, D) = p(θ) · p(D|θ)

p(θ) se denomina distribución a priori de θ. El segundo término se denomina función de

verosimilitud.
En realidad, estamos interesados en determinar p(θ|D) que es la distribución de los
parámetros desconocidos dadas las observaciones conocidas. Como p(θ, D) = p(D, θ),
entonces
p(θ) · p(D|θ) = p(D) · p(θ|D)

de modo que
p(θ) · p(D|θ)
p(θ|D) = .
p(D)
Esto se puede expresar como
p(θ) · L(θ|D)
π(θ|D) = R
θ
p(θ) · L(θ|D)dθ

donde L(θ|D) denota la verosimilitud y

Z
p(D) = p(θ) · L(θ|D)dθ
θ

1
es la constante normalizadora o la distribución predictiva a priori.
De manera alternativa, como la integral anterior no depende de θ, se puede expresar
como
π(θ|D) ∝ p(θ) · L(θ|D)

es decir,

Probabilidad a posteriori ∝ Probabilidad a priori × Verosimilitud.

El problema de inferencia
Se parte una muestra de datos, x = (x1 , . . . , xn ), de modo que la variable aleatoria X
que los genera se asume que depende de unos parámetros θ y que tiene como función de
probabilidad f :
X|θ ∼ f (·)

Se trata de realizar inferencias sobre θ, es decir, estimar su valor.

La Inferencia Clásica presenta las siguientes caracterı́sticas:

– El concepto de probabilidad está limitado a aquellos sucesos en los que se pueden

definir frecuencias relativas.

– θ es un valor fijo (pero desconocido).

– Se pueden usar estimadores de máxima verosimilitud (que se justifican asintótica-

mente) o se pueden usar estimadores insesgados.

– Se usa el concepto de intervalo de confianza:

Ejemplo: Si (1, 3) es un intervalo de confianza al 95 %, significa que si repetimos

el procedimiento muchas veces y calculamos un intervalo de confianza cada vez, el
95 % de los intervalos incluirán θ que es el verdadero valor del parámetro. Esto no
significa que la probabilidad de que θ esté dentro del intervalo sea igual a 0.95.

La definición es una consecuencia del punto de vista de la probabilidad como fre-

cuencia relativa.

– El método de muestreo es muy importante.

2
 La Inferencia Bayesiana se caracteriza por:

– Cada persona tiene sus propias creencias subjetivas previas (a priori) para cualquier
suceso: P (cruz), P (lloverá mañana),...

Nuestras probabilidades pueden ser diferentes porque son nuestras propias medidas
de incertidumbre. La única restricción es que han de ser coherentes (que cumplan
los axiomas de la probabilidad de Kolmogorov).

– θ es un variable y tiene una distribución f (θ). Modificamos nuestras creencias sobre

θ usando el teorema de Bayes.

– La estimación es un problema de decisión. En situaciones distintas se eligen estima-

dores diferentes: se usa la teorı́a de utilidad para elegirlos.

– Un intervalo de credibilidad del 95 % para θ es un intervalo tal que la probabilidad

de que contenga a θ igual a 0,95.

Comentarios y crı́ticas respecto a los métodos bayesianos

– Crı́tica: θ no está claro que tenga que ser siempre una variable.

Argumento: en la práctica, la distribución f (θ) muestra los conocimientos acerca de

θ. El conocimiento cambia y evoluciona con la recogida de los datos.

– Crı́tica: Falta de objetividad. ¿Cómo se puede elegir la distribución a priori cuando

no se tiene información?

Argumento: se pueden establecer métodos objetivos para calcular la distribución

a priori. A menudo, un análisis clásico equivale a un análisis bayesiano con una
distribución a priori no informativa.

– Los aspectos subjetivos son explı́citos en un análisis bayesiano, mientras que en los
métodos clasicos existen aunque no se muestran claramente.

– Es imprescindible hacer un análisis de sensibilidad. Si los resultados cambian cuando

se cambia la distribución inicial, entonces la verosimilitud no da mucha información
y la elección de la distribución inicial es fundamental.

3
Principios fundamentales en Inferencia
El Principio de Verosimilitud

Se trata de determinar estimadores de θ. Si se observan los datos x, toda la información

necesaria para la inferencia está contenida en la función de verosimilitud L(θ|x). Además,
dos funciones de verosimilitud tienen la misma información sobre θ si son proporcionales
entre sı́.
Los métodos bayesianos cumplen el principio de verosimilitud: Si L1 (θ|x) ∝ L2 (θ|x)
entonces, dada una distribución a priori π(θ),

f1 (θ|x) ∝ π(θ)L1 (θ|x) ∝ π(θ)L2 (θ|x) ∝ f2 (θ|x)

No obstante, se puede demostrar que los contrastes clásicos con significación a nivel
fijo (por ejemplo α = 0,05) y los intervalos de confianza no cumplen con este principio.

El Principio de Suficiencia

Un estadı́stico t es suficiente para θ si

L(θ|x) = h(t, θ)g(x).

que recibe el nombre de teorema de factorización de Neyman.

El Principio de Suficiencia afirma que si existe un estadı́stico suficiente, t, dadas dos
muestras, x1 , x2 , que cumplen t(x1 ) = t(x2 ), las conclusiones basadas en x1 y x2 deben
ser iguales.
Todos los métodos estándar de inferencia estadı́stica cumplen el principio de suficien-
cia.

El Principio de Condicionalidad

Suponiendo que se pueden hacer dos posibles experimentos E1 y E2 en relación a θ, y

se elige uno con probabilidad 0,5, entonces la inferencia sobre θ depende sólo del resultado
del experimento seleccionado. La inferencia bayesiana cumple el principio, pero la clásica
no.
El principio de verosimilitud equivale al principio de suficiencia más el principio de
condicionalidad (ver demostración en Lee (2012)).

4
Modelo Beta-Binomial
Este esquema se va a aplicar a problemas relativos a proporciones. Nos interesarán,
por tanto, experimentos con dos posibles resultados, uno de los cuales se designará éxito
y otro fracaso. Existen muchos ejemplos de este tipo de experimentos, por ejemplo:

Estamos diseñando un sistema de control de una lı́nea de producción y consideramos

que los productos pueden ser defectuosos o no.

Al diseñar un sistema de decisión inteligente, por ejemplo en Medicina, podemos

considerar que los pacientes pueden tener o no una enfermedad.

Para diseñar una polı́tica de mantenimiento de un mecanismo de control registramos

si el mecanismo indica situación de peligro o no.

Estamos interesados en evaluar la efectividad de un programa de educación para la

infancia, de modo que consideramos la proporción de niños y niñas que han mejorado
su rendimiento.

Nos interesa obtener información acerca de una proporción determinada, por ejemplo,
sus valores tı́picos o si es menor que un valor prefijado. También se puede considerar el
problema de comparar dos proporciones: si cierta proporción es mayor que otra, o si la
diferencia entre dos proporciones es menor que 0.2, por ejemplo.
La situación que consideramos es la siguiente: se desea recoger y proporcionar infor-
mación sobre p, la proporción de casos en que se produce cierto fenómeno, pudiéndose
dar sólo dos resultados. Disponemos de unas creencias iniciales sobre p, que puede tomar
valores entre 0 y 1.
Se tiene que determinar una distribución suficientemente flexible que modelice estas
creencias sobre proporciones. Una posible distribución es la distribución beta.

NOTA:

Una v.a. tiene distribución beta de parámetros α, β en [0, 1] (y se representa como

X ∼ Be (α, β) ), con α, β > 0 si su función de densidad es

Γ (α + β) α−1
f (x|α, β) = x (1 − x)β−1 I0<x<1
Γ (α) Γ (β)
Recuerda que la función gamma se define como
Z ∞
Γ (α) = xα−1 e−x dx
0

5
y sus propiedades básicas son

Γ (α + 1) = αΓ (α)
Γ (n + 1) = n!

para n ∈ N
Los momentos de la distribución son

Z 1
Γ (α + β) α−1
µ= x x (1 − x)β−1 dx =
0 Γ (α) Γ (β)

Z 1
Γ (α + β) Γ (α + 1) Γ (β) Γ (α + β + 1) α
= x (1 − x)β−1 dx =
Γ (α) Γ (β) Γ (α + β + 1) 0 Γ (α + 1) Γ (β)

Γ (α + β) Γ (α + 1) Γ (β) α
= ·1= .
Γ (α) Γ (β) Γ (α + β + 1) α+β
Del mismo modo se demuestra que

αβ
σ2 = 2 .
(α + β) (α + β + 1)

Cuando α = β = 1, se obtiene la distribución uniforme y, en general, cuando α = β

la distribución es simétrica alrededor de 1/2. Cuando α < β se concentra la probabilidad
hacia la izquierda, y a la inversa cuando α > β. Cuanto mayores son α y β las densidades
están más concentradas en algún punto.
En la gráfica siguiente se muestran ejemplos sobre distintas formas que presenta va-
riando sus parámetros:

6
7
 Las gráficas anteriores se han dibujado con R:

z = seq (0 ,1 , length =1000)

par ( mfrow = c (3 ,3))

# Parametros : alfa =5 beta =5

plot (z , dbeta (z ,5 ,5) , ylab = " Beta (5 ,5) " , type = " l " )

# Parametros : alfa =5 beta =1

plot (z , dbeta (z ,5 ,1) , ylab = " Beta (5 ,1) " , type = " l " )

# Parametros : alfa =5 beta =3

plot (z , dbeta (z ,2 ,3) , ylab = " Beta (2 ,3) " , type = " l " )

# Parametros : alfa =1 beta =5

plot (z , dbeta (z ,1 ,5) , ylab = " Beta (1 ,5) " , type = " l " )

# Parametros : alfa =1 beta =1

plot (z , dbeta (z ,1 ,1) , ylim = c (0 ,1 .5 ) , ylab = " Beta (1 ,1) " , type = " l " )

# Parametros : alfa =1 beta =3

plot (z , dbeta (z ,1 ,3) , ylab = " Beta (1 ,3) " , type = " l " )

# Parametros : alfa =3 beta =5

plot (z , dbeta (z ,3 ,5) , ylab = " Beta (3 ,5) " , type = " l " )

# Parametros : alfa =3 beta =1

plot (z , dbeta (z ,3 ,1) , ylab = " Beta (3 ,1) " , type = " l " )

# Parametros : alfa =3 beta =3

plot (z , dbeta (z ,3 ,3) , ylab = " Beta (3 ,3) " , type = " l " )

En nuestro caso, se tiene que


Γ (α + β) α−1
 Γ (α) Γ (β) p (1 − p)β−1 para 0≤p≤1



f (p|α, β) =



 0 resto

α
E (p) = ,
α+β
αβ
V ar (p) = 2 ,
(α + β) (α + β + 1)
α−1
M oda (p) = .
α+β−2
Por ejemplo, para α = β = 1 se tiene la distribución uniforme que modeliza la igno-
rancia acerca de p.

8
Planteamiento del modelo beta-binomial
El problema que nos planteamos es el siguiente. Suponemos un experimento que con-
siste en observar n casos independientes, registrándose el número de casos favorables que
se presentan. La verosimilitud (o el modelo) es, en este caso, binomial, teniéndose
 
n x
P (X = x|p) = p (1 − p)n−x ,
x
para x = 0, 1, . . . n.
Se realiza, por tanto, el experimento y supongamos que se producen x éxitos. Nuestro
interés se centra en estimar la proporción p de éxitos dada la muestra.
Desde el punto de vista bayesiano asumimos que tenemos una información previa sobre
p que se modeliza mediante la distribución a priori de p que hemos supuesto que se recoge
mediante una densidad beta. Actualizamos, entonces, nuestras creencias sobre p aplicando
el teorema de Bayes. Se tiene
 
f (p) P (x|p) Γ (α + β) α−1 β−1 n x
f (p|x) = ∝ p (1 − p) · p (1 − p)n−x
P (x) Γ (α) Γ (β) x

∝ px+α−1 (1 − p)n−x+β−1 ,

siendo p ∈ (0, 1) y f (p|x) = 0 para p ∈

/ (0, 1) .
Resulta, por tanto, que la distribución a posteriori de p sigue una beta de parámetros
(x + α) y (n − x + β) .
Se puede emplear esa distribución para resolver los problemas tı́picos que se presentan
en inferencia.

Estimación puntual
Se trata de dar un valor-resumen representativo de la distribución a posteriori. Las
más habituales son:

1. Media a posteriori,
x+α
n+α+β
2. Moda a posteriori,
x+α−1
n+α+β−2
3. Mediana a posteriori, que es la solución p∗ de la ecuación
Z p∗
Γ (n + α + β) 1
px+α−1 (1 − p)n−x+β−1 dp = .
0 Γ (α + x) Γ (n − x + β) 2
En este caso, se hace necesario usar cálculo numérico.

9
Ejemplo en R:
# Se calcula el cuantil del 50 % de una beta 10 , 15
qbeta (0 .5 ,10 ,15)

[1] 0 .3972924

Estimación por intervalos

De manera informal, se trata de dar un intervalo que, con cierta probabilidad r, cubra
el valor de p, esto es, un intervalo con extremos a y b que satisfaga lo siguiente
Z b
f (p|x) dp = r.
a

1−r
Normalmente, se escogen a y b de manera que a deja probabilidad 2
a su izquierda
1−r
y 2
a su derecha.
Por ejemplo, si r = 0,90, a y b satisfacen
Z a
f (p|x) dp = 0,05
0
Z 1
f (p|x) dp = 0,05.
b

Por ejemplo, en R:
a = qbeta (0 .05 ,4 ,5)
b = qbeta (0 .05 ,4 ,5 , lower.tail = F )

Contraste de hipótesis
En muchas ocasiones se está interesado en un modelo especı́fico o subconjuntos de
modelos que denominamos hipótesis nula (H0 ) frente al resto de modelos que se denominan
hipótesis alternativa (H1 ). Por ejemplo, podemos contrastar

H0 : p = 0,5

frente a
H1 : p 6= 0,5;

o podemos contrastar
H0 : p ≤ 0,6

frente a
H1 : p > 0,6.

Contrastar H0 frente a H1 , significa calcular la probabilidad a posteriori de cada una

de las hipótesis y quedarse con la hipótesis más probable.

10
 En el segundo ejemplo propuesto se calcuları́a
Z 0,6
P (H0 |x) = f (p|x) dp,
0
Z 1
P (H1 |x) = f (p|x) dp,
0,6

y se dirı́a que los datos apoyan H0 (frente a H1 ) si

P (H0 |x) > P (H1 |x)

(o equivalentemente, P (H0 |x) > 0,5). El primer ejemplo es algo más delicado, pues
P (H0 |x) = 0, al ser una distribución continua. Una solución rigurosa requiere cálcu-
los más sofisticados. Alternativamente, se puede calcular un intervalo I de probabilidad r,
centrado en la media a posteriori de p, y si 0,5 ∈ I, se dice que hay evidencia a favor de la
hipótesis nula. En otro caso, se dice que hay evidencia a favor de la hipótesis alternativa.

Ejemplo

Supongamos que estamos interesados en estudiar los hábitos de sueño de los estudian-
tes de un cierto centro. Parece ser que los médicos recomiendan un mı́nimo de 8 horas de
sueño para una persona adulta, con lo cual el estudio se plantea en términos de averiguar
la proporción de de estudiantes que duermen al menos 8 horas. Llamaremos p a dicha
proporción.
Se toma una muestra de 27 estudiantes de modo que 11 de ellos duermen al menos
8 horas y el resto no. Nos planteamos hacer inferencias sobre la proporción p teniendo
en cuenta la información previa que se tiene. Además interesa predecir el número de
estudiantes que duermen al menos 8 horas si se toma una nueva muestra de 20 estudiantes.
Supongamos que consideramos una distribución a priori beta:
p = seq (0 ,1 , length =500)
s =11
f =16
a =3 .4
b =7 .4

prior = dbeta (p , a , b )
like = dbeta (p , s +1 , f +1)
post = dbeta (p , a +s , b + f )

plot (p , post , type = " l " , ylab = " Density " , lty =2 , lwd =3 , col =4)
lines (p , like , lty =1 , lwd =3 , col =2)
lines (p , prior , lty =3 , lwd =3 , col =6)
legend (0 .7 , 4 , c ( " A priori " ," Verosimilitud " ," A posteriori " ) ,
lty = c (3 ,1 ,2) , lwd = c (3 ,3 ,3) , col = c (6 ,2 ,4))

11
 Se pueden considerar resúmenes de la distribución: ¿Cuál es P (p ≥ 0,5|datos)? ¿Cuál
es el intervalo de confianza al 90 % para el verdadero valor de p?
s = 11
f = 16
a = 3 .4
b = 7 .4

1 - pbeta (0 .5 , a +s , b + f )

[1] 0 .06842569

qbeta ( c (0 .05 , 0 .95 ) , a +s , b + f )

[1] 0 .2562364 0 .5129274

Otra alternativa es hacerlo mediante simulación:

s = 11
f = 16
a = 3 .4
b = 7 .4

ps = rbeta (1000 , a +s , b + f )

sum ( ps >= 0 .5 ) / 1000

12
[1] 0 .072

quantile ( ps , c (0 .05 , 0 .95 ))

5% 95 %
0 .2538683 0 .5107199

lattice :: histogram ( ps , xlab = " p " , main = " " )

13
Resolución con MCMCpack

library ( MCMCpack )

posterior = MCbinomialbeta (11 , 27 , mc =10000)

summary ( posterior )

X11 ()
plot ( posterior )

grid = seq (0 ,1 ,0 .01 )

X11 ()
plot ( grid , dbeta ( grid , 1 , 1) , type = " l " , col = " red " ,
lwd =3 , ylim = c (0 ,4 .5 ) , xlab = expression ( pi ) , ylab = " Densidad " )
lines ( density ( posterior ) , col = " blue " , lwd =3)
legend ( .75 , 3 .6 , c ( " A priori " , " A posteriori " ) ,
lwd =3 , col = c ( " red " , " blue " ))

Se obtiene

Iterations = 1:10000
Thinning interval = 1
Number of chains = 1
Sample size per chain = 10000

1 . Empirical mean and standard deviation for each variable ,

plus standard error of the mean :

Mean SD Naive SE Time - series SE

0 .4130949 0 .0901232 0 .0009012 0 .0009012

2 . Quantiles for each variable :

2 .5 % 25 % 50 % 75 % 97 .5 %
0 .2407 0 .3501 0 .4111 0 .4750 0 .5923

14
15
16
Resolución con SAS

Entrar en
https://odamid.oda.sas.com/SASStudio

OPTIONS nodate ;
TITLE ’ Ejemplo BetaBinomial Bayes ’;

DATA cosas ;
INPUT sopa @@ ;
DATALINES ;
1 1 1 1 1 1 1 1 1 1 1 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0
;

PROC mcmc data = cosas nmc =10000 outpost = sale ;

/* Valor i n i c i a l de la ca den a */
parm p 0.2;

/* A priori no i n f o r m a t i v a :
prior p ∼ u n i f o r m (0 ,1) ; */

/* A priori c o n j u g a d a */
prior p ∼ beta (4 ,4) ;

model sopa ∼ Binary ( p ) ;

RUN ;

17
Ejemplo

Un distribuidor anuncia que el 95 % de sus ordenadores no se tiene que reparar durante

el año de garantı́a. Desde el punto de vista bayesiano se puede modelizar esta creencia
mediante una Be (α = 4,75, β = 0,25) , ya que si

p ∼ Be (4,75, 0,25) ,

entonces
4,75
E (p) = = 0,95,
4,75 + 0,25
V ar (p) = 0,0016.

Debido a la aparente calidad, compramos 20 ordenadores, de los cuales 12 requieren re-

paración durante el año de garantı́a. Nuestra distribución a posteriori sobre la proporción
es entonces ahora
Be (4,75 + 8, 0,25 + 12)

La figura muestra las distribuciones a priori y a posteriori.

18
Estimadores puntuales.
Algunos valores que resumen la distribución a posteriori son:
12,75
Media a posteriori: E (p|x) = 25
= 0,51
12,75−1
Moda a posteriori: 25−2
= 0,51087
Mediana a posteriori: 0.5102709 (en R: qbeta(0.5,12.75,12.25)).

Estimación por intervalos:

Un intervalo de probabilidad 0.95 es
En R:
p = c (0 .025 , 0 .975 )
qbeta (p , 12 .75 , 12 .25 )

Obteniéndose [0,32; 0,70]

Contraste de Hipótesis:
Contrastamos la hipótesis H0 : p ≥ 0,95 frente a H1 : p < 0,95, que se corresponderı́a
con contrastar el anuncio.
Como, al considerar la orden en R, pbeta(0.95,12.75,12.25), se obtiene 1, entonces
se tiene que

P (H1 |x) ≈ 1
P (H0 |x) ≈ 0,

por lo que los datos sugieren que el distribuidor está exagerando.

Predicción
Dados los datos x = (x1 , . . . , xn )0 , supongamos que se quiere predecir el valor de una
futura observación Xn+1 . Para ello, hay que calcular la distribución predictiva. Se puede
considerar tanto la distribución predictiva a priori como a posteriori dependiendo de si
se usa una u otra distribución de los parámetros.
Sus expresiones son
Z
f (xn+1 |x) = f (xn+1 |θ, x)f (θ) dθ
Z
f (xn+1 |x) = f (xn+1 |θ, x)f (θ|x) dθ.

En nuestra situación, los Xi |θ son intercambiables, es decir, su distribución es igual

cuando se les cambia de orden y lo anterior se puede simplificar como
Z
f (xn+1 |x) = f (xn+1 |θ)f (θ) dθ
Z
f (xn+1 |x) = f (xn+1 |θ)f (θ|x) dθ.

19
Ejemplo:
Supongamos que se tira una moneda equilibrada y que salen 9 cruces y 3 caras (con-
sidero como éxito ≡ cruz ). Suponiendo una distribución a priori uniforme, se tiene que
la distribución a posteriori de p sigue una distribución beta de parámetros (x + α) y
(n − x + β): Beta(10, 4).
A continuación, se quiere predecir el número de cruces que se obtendrı́an en diez
tiradas más.
Se tiene que, ahora,

Xf ut |θ ∼ Bin(10, θ),

entonces
!
Z 1
10
f (xf ut |x) = θxf ut (1 − θ)10−xf ut ×
0 xf ut
Γ(14)
× θ10−1 (1 − θ)4−1 dθ
Γ(10)Γ(4)

!
10 Γ(14)
= ×
xf ut Γ(10)Γ(4)
Z 1
× θ10+xf ut −1 (1 − θ)14−xf ut −1 dθ
0

!
10 Be(10 + xf ut , 14 − xf ut )
=
xf ut Be(10, 4)

donde se ha usado la notación

Γ(α)Γ(β)
Be(α, β) =
Γ(α + β)

En general, si estamos interesados en r éxitos en los m siguientes ensayos haremos

Z 1 
m r
P [r éxitos en m ensayos | x] = p (1 − p)m−r f (p|x) dp =
0 r

Z 1  
m r Γ (n + α + β)
= p (1 − p)m−r px+α−1 (1 − p)n−x+β−1 dp =
0 r Γ (α + x) Γ (n − x + β)

 
m Γ (n + α + β) Γ (r + x + α) Γ (m − r + n − x − β)
= .
r Γ (α + x) Γ (n − x + β) Γ (m + n + α + β)

Se puede evaluar la media de Xf ut |x sin tener que evaluar la distribución predictiva.

20
NOTA: se tenı́a que,
Ez [Z] = Ey [Ez [Z|Y ]]

para cualquier variable Z y Y .

Ası́, respecto al ejemplo anterior,

Ex [Xf ut |x] = Eθ [Ex [Xf ut |θ]]

como Xf ut |θ ∼ Bin(10, θ) entonces Ex [Xf ut |θ] = 10θ, y como θ|x ∼ Be(10, 4) entonces
10
Eθ [θ|x] = 14
, por tanto
10
Ex [Xf ut |x] = Eθ [Ex [Xf ut |θ]] = 10 × ≈ 7,141
14
Para evaluar la varianza predictiva, se usa la fórmula

Vz [Z] = Ey [Vz [Z|Y ]] + Vy [Ez [Z|Y ]]

= Eθ [Vx [Xf ut |θ]] + Vθ [Ex [Xf ut |θ]]

Evaluamos ambas partes.

Vx [Xf ut |θ] = 10θ(1 − θ) = 10(θ − θ2 )

Eθ [Vx [Xf ut |θ]] = 10 Eθ [θ|x] − Eθ [θ2 |x]

= 10 Eθ [θ|x] − (Eθ [θ|x]2 + Vθ [θ|x])

 2
10 10
= 10 − −
10 + 4 10 + 4

10 · 4
− =
(10 + 4)2 (10 + 4 + 1)
40
=
21
y como Ex [Xf ut |θ] = 10θ,

Vθ [Ex [Xf ut |θ]] = 100Vθ [θ|x]

10 · 4
= 100
(10 + 4)2 (10 + 4 + 1)
200
=
147
40 200
Entonces tenemos Vx [Xf ut |x] = 21
+ 147
≈ 3,3.
Predicción en el ejemplo de los ordenadores:
La probabilidad de que el siguiente ordenador que adquiramos no deba repararse en
el periodo de garantı́a es igual a la media a posteriori ya que es una Ber(θ), luego
Z 1
12,75
f (xf ut |x) = θf (θ|x)dθ = E(θ|x) = = 0,51.
0 25

21
Asignación de distribuciones beta
Se trata, ahora de encontrar procedimientos para asignar la distribución a priori,
supuesto que es Be (α, β) . Esto implica escoger α y β, lo que implica extraer dos juicios
de un experto.
Una forma de proceder es la siguiente. Primero pedimos al experto que indique la
α
probabilidad r de éxito en el primer ensayo. Sabemos que ésta es α+β
.
Después, se pide al experto que suponga que el primer ensayo fue un éxito y que
proporcione la probabilidad r2 de éxito en el segundo ensayo; la densidad actualizada por
α+1
el primer éxito es Be (α + 1, β) , por lo que la probabilidad r2 será α+β+1
.
Se resuelve, entonces, el sistema
α
r1 = ,
α+β
α+1
r2 =
α+β+1
y se obtiene,
r1 (1 − r2 )
α= ,
r2 − r1
(1 − r1 ) (1 − r2 )
β= .
r2 − r1

Nota:

Si se emplea la Be (1, 1) se corresponde a la U (0, 1) modeliza la ignorancia o total

incertidumbre. En este sentido, cuando tengamos mucha incertidumbre deberı́amos favo-
recer distribuciones a priori con valores pequeños de α y β que facilitan la absorción de
evidencia.

Ejemplo

Supongo que el experto afirma que r1 = 0,7 y r2 = 0,75. Entonces,

0,7 (1 − 0,75)
α= = 3,5,
0,75 − 0,7
(1 − 0,7) (1 − 0,75)
β= = 1,5.
0,75 − 0,7
Observemos que todos los cálculos que hagamos quedarán afectados por r1 y r2 , a
través de α y β, especialmente con pocos datos.

22
Alternativas para la distribución a priori
En una muestra aleatoria tomada de una clase de 27 alumnos, se considera la variable
aleatoria X que es igual a 1 si los alumnos duermen más de 8 horas al dı́a ó 0 en caso
contario. La verosimilitud es por tanto binomial de parámetro desconocido p.
Supongamos que la distribución de probabilidad a priori para p se denota como g(p).
Si denominamos como éxito al hecho de dormir más de 8 horas, la verosimilitud es

L(p) ∝ ps (1 − p)f

de modo que la densidad a posteriori es

g (p|datos) ∝ g(p) · ps (1 − p)f

Luego interesa determinar la distribución a priori.

Alternativa 1: A priori discreta.

Supongamos que se piensa que los valores posibles para p son

0,05; 0,15; 0,25; 0,35; 0,45; 0,55; 0,65; 0,75; 0,85; 0,95

basado en su información previa le asigna los siguientes pesos:

2, 4, 8, 8, 4, 2, 1, 1, 1, 1

que se pueden convertir en probabilidades a priori normalizando los pesos para que sumen
1.
En R esto se puede hacer ası́:
p = seq (0 .05 , 0 .95 , by =0 .1 )
prior = c (2 , 4 , 8 , 8 , 4 , 2 , 1 , 1 , 1 , 1)
prior = prior / sum ( prior )
plot (p , prior , type = " h " , ylab = " Probabilidad a priori " , col =4)

23
 En nuestro caso hay 11 estudiantes que duermen bien, por lo que

L(p) ∝ p11 (1 − p)16

Se puede observar que es una distribución beta: Be(12, 17)

Para calcular las probabilidades a posteriori usamos la librerı́a LearnBayes:
library ( LearnBayes )
data = c (11 , 16)
post = pdisc (p , prior , data )
cbind (p , prior , post )
plot (p , post , type = " h " , ylab = " Probabilidad a posteriori " , col =2)

Se obtiene

24
 p prior post
[1 ,] 0 .05 0 .06250 2 .882642e -08
[2 ,] 0 .15 0 .12500 1 .722978e -03
[3 ,] 0 .25 0 .25000 1 .282104e -01
[4 ,] 0 .35 0 .25000 5 .259751e -01
[5 ,] 0 .45 0 .12500 2 .882131e -01
[6 ,] 0 .55 0 .06250 5 .283635e -02
[7 ,] 0 .65 0 .03125 2 .976107e -03
[8 ,] 0 .75 0 .03125 6 .595185e -05
[9 ,] 0 .85 0 .03125 7 .371932e -08
[10 ,] 0 .95 0 .03125 5 .820934e -15

Si dibujamos ambas para comparar cómo ha variado la distribución de p después de

obtener los datos:
par ( mfrow = c (2 ,1))
plot (p , prior , type = " h " , ylab = " Probabilidad a priori " , col =4)
plot (p , post , type = " h " , ylab = " Probabilidad a posteriori " , col =2)

25
Alternativa 2: A priori por un histograma. Se tendrı́a una versión continua de la anterior.
Supongamos que se puede asignar una distribución a priori sobre intervalos. Asumimos
que nos basta con 10 para el rango de p: (0, 0.1), (0.1, 0.2),... (0.9, 1) de modo que
asignamos como pesos:
2, 4, 8, 8, 4, 2, 1, 1, 1, 1

library ( LearnBayes )
midpt = seq (0 .05 , 0 .95 , by = 0 .1 )
prior = c (2 , 4 , 8 , 8 , 4 , 2 , 1 , 1 , 1 , 1)
prior = prior / sum ( prior )
p = seq (0 , 1 , length = 500)

plot (p , histprior (p , midpt , prior ) , type = " l " ,

ylab = " Densidad a priori " , ylim = c (0 , .25 ))

26
s = 11
f = 16
like = dbeta (p , s +1 , f +1)
post = like * histprior (p , midpt , prior )
plot (p , post , type = " l " , ylab = " Densidad a posteriori " , col =4)

27
 Se puede obtener una muestra simulada de la densidad a posteriori, convirtiendo los
productos que aparecen en la rejilla anterior en probabilidades. Se toma una muestra con
reemplazamiento de de la rejilla:

post = post / sum ( post )

ps = sample (p , replace = TRUE , prob = post )
hist ( ps , xlab = " p " )

28
Ejemplo

En el ejemplo de las horas de sueño se puede considerar la distribución predictiva a

posteriori:
library ( LearnBayes )
s = 11
f = 16
a = 3 .4
b = 7 .4
ab = c ( a +s , b + f )
m = 20
ys = 0:20
pred = pbetap ( ab , m , ys )
cbind ( ys , pred )

Se obtiene

29
 ys pred
[1 ,] 0 0 .0006
[2 ,] 1 0 .0040
[3 ,] 2 0 .0143
[4 ,] 3 0 .0348
[5 ,] 4 0 .0653
[6 ,] 5 0 .1002
[7 ,] 6 0 .1299
[8 ,] 7 0 .1456
[9 ,] 8 0 .1431
[10 ,] 9 0 .1242
[11 ,] 10 0 .0957
[12 ,] 11 0 .0655
[13 ,] 12 0 .0398
[14 ,] 13 0 .0212
[15 ,] 14 0 .0099
[16 ,] 15 0 .0040
[17 ,] 16 0 .0013
[18 ,] 17 0 .0004
[19 ,] 18 0 .0001
[20 ,] 19 0 .0000
[21 ,] 20 0 .0000

De manera alternativa, para obtener las predicciones se puede simular p∗ de la distri-

bución a posteriori, obtiéndose las predicciones a partir de una distribución binomial con
parámetro igual a p∗ :
p = rbeta (1000 , a +s , b + f )
y = rbinom (1000 , 20 , p )
table ( y )

y
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
1 3 12 37 60 113 127 143 147 124 110 57 39 15 7 2 3

Se puede dibujar la distribución de las predicciones

freq = table ( y )
if ( names ( freq )[1]== " 0 " ){
# Si aparece alguna observacion 0
ys = c (0: max ( y ))
} else {
# Si NO aparece ninguna observacion 0
ys = c (1: max ( y ))}
predprob = freq / sum ( freq )
plot ( ys , predprob , type = " h " , xlab = " y " , ylab = " Probabilidad predictiva " )

30
 Si se quiere resumir la distribución predictiva mediante un intervalo que cubra el 90 %
de la función, se puede usar la función discint.
Se obtienen los valores que incluye el intervalo y la probabilidad exacta de los mismos.
# Highest probability interval for a
# discrete distribution
covprob = 0 .9
dist = cbind ( ys , as.vector ( predprob ))
discint ( dist , covprob )

$ prob
[1] 0 .92

$ set
[1] 4 5 6 7 8 9 10 11 12

31
Comportamiento del modelo con muestras grandes
Ejemplo

Dos ingenieros se enfrentan a un problema de control de calidad; para ello se fijan en la

proporción de piezas defectuosas. El primero es inexperto sobre el proceso de producción
de modo que se modeliza sus creencias sobre dicha proporción mediante una Be (1, 1) (es
decir, con una U (0, 1)). El segundo tiene más experiencia y modeliza sus creencias sobre
la proporción mediante una Be (10, 70) . Ambos observan tres piezas, siendo una de ellas
defectuosa. Las creencias del primer ingeniero pasan a ser Be (2, 3) , cambiando radical-
mente, mientras que las del segundo pasan a ser Be (11, 72) , sin cambiar prácticamente
nada.
Después de observar 220 piezas en buen estado y 160 defectuosas, las creencias del
primero pasan a ser Be (161, 221) y las del segundo pasan a ser Be (170, 290) , pareciéndose
ahora las creencias de ambos ingenieros.

El ejemplo describe una situación caracterı́stica de un modelo beta-binomial que se da,

también, en otros modelos: al aumentar la evidencia aportada por los datos, la influencia
de la distribución a priori disminuye.

32
 En efecto, como la distribución a posteriori es Be (x + α, n − x + β) , cuando n es muy
grande, la media a posteriori es
x
x+α n
+ αn x
= α+β
≈
α+β+n 1+ n n

y la varianza a posteriori
(x + α) (β + n − x) x(n − x) x x 1
≈ ≈ 1 −
(α + β + n)2 (α + β + n + 1) n3 n n n

Por tanto, para n suficientemente grande, al tender la varianza hacia 0, la distribución

se concentrará más y más en torno a x/n. De hecho, el resultado es más preciso puesto
que se puede probar que para n suficientemente grande, la distribución a posteriori es
 q

x x x 1
aproximadamente normal N n ; n
1 − n n . Este resultado simplifica las cuestiones
computacionales.
Ası́, por ejemplo, se puede considerar:

x
(i ) Estimación puntual con la media a posteriori:
n
(ii ) Estimación por intervalos
" r  r  #
x x x 1 x
 x x 1

I= − z α2 1− ; + z α2 1−
n n n n n n n n

(iii ) Contraste de hipótesis (por ejemplo, H0 : p = 0,5): Se rechaza H0 si 0,5 ∈

/ I.

El resultado anterior es bastante general y se describe técnicamente como la normali-

dad asintótica de la distribución a posteriori.
Una alternativa consiste en emplear una simulación de la distribución a posteriori.
Este caso es muy sencillo: basta simular de una Be (α + x, β + n − x) y sustituir mo-
mentos poblacionales por momentos muestrales.
Por ejemplo, podemos simular 1000 observaciones de tal distribución (x1 , x2 , . . . , x1000 )
y
P
xi
1. Aproximar la media a posteriori por n
i

2. Aproximar el intervalo de probabilidad 1−α (pongamos, para fijar ideas que α = 0,1)
 
mediante X(50) , X(950)

Ejemplo

Generamos 1000 observaciones de una Be (12,75, 12,25) . En R se harı́a:

33
 x = rbeta (1000 ,12 .75 ,12 .25 )
mean ( x )
y = sort ( x )
y [975]
y [25]

La media muestral es 0,5110.

Al ordenar la muestra se obtiene

X(975) = 0,7017
X(25) = 0,3197

con lo que un intervalo de probabilidad a posteriori de 0,95 es [0,3197; 0,7017].

Ejemplo con R

Se necesita instalar previamente la librerı́a Bolsdtad de R.

library ( Bolstad )
# Ejemplo con 6 exitos en 8 ensayos con una a priori
# uniforme Beta (1 ,1)
x11 ()
binobp (6 ,8)

34
# 6 exitos en 8 ensayos con una a priori Beta (0 .5 ,6)
x11 ()
binobp (6 ,8 ,0 .5 ,6)

# 4 exitos observados in 12 ensayos con una a priori Beta (3 ,3)

# Se dibujan todas las graficas
x11 ()
results = binobp (4 ,12 ,3 ,3)

35
x11 ()
par ( mfrow = c (3 ,1))

y.lims = c (0 ,1 .1 * max ( results $ posterior , results $ prior ))

plot ( results $ pi , results $ prior , ylim = y.lims , type = " l " ,

xlab = expression ( pi ) , ylab = " Density " , main = " Prior " )

polygon ( results $ pi , results $ prior , col = " red " )

plot ( results $ pi , results $ likelihood , ylim = c (0 ,0 .25 ) , type = " l " ,

xlab = expression ( pi ) , ylab = " Density " , main = " Likelihood " )

polygon ( results $ pi , results $ likelihood , col = " green " )

plot ( results $ pi , results $ posterior , ylim = y.lims , type = " l " ,

xlab = expression ( pi ) , ylab = " Density " , main = " Posterior " )

polygon ( results $ pi , results $ posterior , col = " blue " )

36
37
Comparación de proporciones
Se estudia, ahora, el problema de la comparación entre dos proporciones. Para ello se
consideran dos poblaciones que dan lugar a dos experimentos, cada uno con una propor-
ción p1 y p2 de éxito.
Deseamos comparar p1 y p2 . Por ejemplo, si p1 ≥ p2 o viceversa.
Suponemos que se tiene una muestra de tamaño n1 de la primera población con x1
éxitos. La verosimilitud es, en este caso,
 
n 1 x1
P (X1 = x1 |p1 ) = p1 (1 − p1 )n1 −x1 .
x1
Además, se tiene otra muestra de tamaño n2 de la segunda población con x2 éxitos. La
verosimilitud es, en este caso,
 
n 2 x2
P (X2 = x2 |p2 ) = p (1 − p2 )n2 −x2 .
x2 2

Supondremos que hay independencia entre ambas muestras (a veces esta hipótesis no
será adecuada, como ocurre, por ejemplo, en muestras apareadas), con lo cual la verosi-
militud conjunta será el producto de las verosimilitudes marginales
   
n1 x1 n1 −x1 n2
P (X1 = x1 , X2 = x2 |p1 , p2 ) = p1 (1 − p1 ) px2 2 (1 − p2 )n2 −x2 .
x1 x2
Para cada pi disponemos de información a priori que modelizamos mediante una dis-
tribución Be (αi , βi ) . Supondremos que las pi son independientes. De nuevo esta hipótesis
no será siempre adecuada; por ejemplo, podemos creer que si una proporción es grande,
es probable que la otra también lo será, o, a la inversa, al ser un proporción grande, la
otra podrı́a ser pequeña. Ası́, la distribución a priori conjunta es

f (p1 , p2 ) ∝ pα1 1 −1 (1 − p1 )β1 −1 pα2 2 −1 (1 − p2 )β2 −1 I[0,1] (p1 ) I[0,1] (p2 ) .

Calculamos la distribución a posteriori de p1 y p2 multiplicando por la verosimilitud,

f (p1 , p2 |x1 , x2 ) ∝ pα1 1 +x1 −1 (1 − p1 )n1 −x1 +β1 −1 pα2 2 +x2 −1 (1 − p2 )n2 −x2 +β2 −1 I[0,1] (p1 ) I[0,1] (p2 ) .

Se comprueba, además, que

f (p1 |x1 , x2 ) ∝ pα1 1 +x1 −1 (1 − p1 )n1 −x1 +β1 −1 I[0,1] (p1 ) ,

f (p2 |x1 , x2 ) ∝ pα2 2 +x2 −1 (1 − p2 )n2 −x2 +β2 −1 I[0,1] (p2 ) ,

con lo que p1 y p2 son también independientes a posteriori. Este resultado es, de hecho,
general.

38
 La distribución a posteriori se emplea en sentido similar a como lo hemos hecho en el
caso de una sola proporción. Por ejemplo, podemos emplearla para el cálculo de probabi-
lidades de conjuntos interesantes:

Conjuntos rectangulares
[c1 , d1 ] × [c2 , d2 ]
Utilizando la independencia de la distribución a posteriori, se tiene
Z d1 Z d2
f (p1 , p2 |x1 , x2 ) dp1 dp2 =
c1 c2
Z d1  Z d2 
= f (p1 |x1 , x2 ) dp1 × f (p2 |x1 , x2 ) dp2 ,
c1 c2

y para calcular cada una de las probabilidades aplicamos los métodos ya vistos.

Conjuntos triangulares
Se está más interesado en calcular probabilidades a posteriori de conjuntos trian-
gulares, y más precisamente, de conjuntos triangulares de la forma (p1 − p2 ) ≥ c. La
probabilidad se puede calcular numéricamente y es igual a
Z Z
f (p1 , p2 |x1 , x2 ) dp1 dp2 .
{p1 −p2 ≥c}

Alternativamente, se puede emplear la normalidad asintótica: si n1 y n2 son suficien-

temente grandes
   
x1 x 1 x1 1
p1 |x1 , x2 ∼ N , 1−
n1 n1 n1 n1
   
x2 x 2 x2 1
p2 |x1 , x2 ∼ N , 1−
n2 n2 n2 n2
y son independientes. Se tiene, por las propiedades de las distribuciones normales, que
     
x1 x2 x1 x1 1 x2 x2 1
p1 − p2 |x1 , x2 ∼ N − , 1− + 1−
n1 n2 n1 n1 n1 n2 n2 n2
y el cálculo de la probabilidad pedida es
   
x1 x2
 c− n1
− n2 
1 − φ
r    
.

x1 x1 1 x2 x2 1
n1
1− n1 n1
+ n2
1− n2 n2

Ası́, un intervalo de probabilidad a posteriori (1 − α) serı́a

  s    
x1 x 2 x1 x1 1 x2 x2 1
− ± z α2 1− + 1− .
n1 n2 n1 n1 n1 n2 n2 n2

39


Otra posibilidad es emplear simulación. Generamos k observaciones p1i , p2i , . . . pki de
pi |x1 , x2 , para i = 1, 2 y se aproxima la probabilidad mediante

cardinal {pi1 − pi2 ≥ c}

.
k
Mediante simulación podemos generar, por ejemplo, 1000 observaciones de p1 y p2 (pij ,
j ∈ {1, 2} , i ∈ {1, 2, . . . , 1000}), ordenar los valores ri = pi1 − pi2 y emplear el intervalo
 
r(50) , r(950) para dar un intervalo aproximado de probabilidad de 0.9.
Por último, se puede contrastar hipótesis, por ejemplo,

H0 : p1 − p2 = 0

frente a
H1 : p1 − p2 6= 0.

Para ello, calculamos un intervalo de probabilidad a posteriori igual 0.95; si incluye a

0, decimos que los datos justifican la hipótesis de que ambas proporciones son iguales.

Ejemplo

Un grupo de investigadores evaluó si un grupo de personas en la república del Congo

tenı́a el virus del sida y registró la frecuencia con que habı́an empleado preservativos
el año anterior. La tabla recoge el número de personas que pasaron de seronegativos a
seropositivos, en función del uso de preservativos.

Uso Nunca <25 % 26-50 % 50-75 % >75 %

Sero+ 74 19 7 0 0
Sero- 214 36 15 2 6
Total 288 55 22 2 6

Agrupamos las observaciones entre personas que usan muy frecuentemente preserva-
tivos (≥75 %) y las que no lo usan (<75 %)

<75 % ≥75 %
Sero+ 100 0
Sero- 267 6
367 6

La cuestión que nos planteamos es si el uso alto de preservativos previene la transmisión

del virus del sida. Consideramos distribuciones a priori Be (1, 1) independientes para
ambas proporciones pSI , pN O y calculamos un intervalo de probabilidad a posteriori 0.95
para d = pN O − pSI .

40
 La densidad a posteriori para pN O es Be (101, 268) y para pSI es Be (1, 7) . Las medias
y varianzas a posteriori son, respectivamente,
101
E (pN O |datos) = = 0,2737
369
V ar (pN O |datos) = 0,0005
1
E (pSI |datos) = = 0,125
8
V ar (pSI |datos) = 0,0121,

con lo que

E (pN O − pSI |datos) = 0,149

V ar (pN O − pSI |datos) = 0,0127.

Suponiendo válido el argumento asintótico resulta que, aproximadamente,

pN O − pSI |datos ∼ N (0,149, 0,0127) .

Como, para la N (0, 1) , P (−1,96, 1,96) = 0,95, un intervalo de probabilidad 0.95 es

h p p i
0,149 − 1,96 0,0127, 0,149 + 1,96 0,0127 =
= [−0,07, 0,37]

y 0 ∈ [−0,07, 0,37] .
Con esto ya se ha respondido a la pregunta formulada, pero podemos completar el
argumento.
Como P (−1, 1) = 0,68, un intervalo de probabilidad 0.68 es
h p p i
0,149 − 1 0,0127, 0,149 + 1 0,0127 =
= [0,04, 0,26]

y0∈
/ [0,04, 0,26]
 
0−0,149
Además, P (pN O ≤ pSI ) = P (pN O − pSI ≤ 0) = Φ √
0,0127
= Φ (−1,318) = 1 −
Φ (1,318) = 0,1.
Ası́, la hipótesis de que usar los preservativos es efectivo tiene algo de apoyo, pero no
mucho.
Procedemos en forma similar mediante simulación. Generamos 1000 observaciones de
una Be (101, 268) y de una Be (1, 7) . Se obtienen las diferencias, se ordenan de menor a
mayor y se calcula el intervalo correspondiente. Para una probabilidad del 0.95 se tiene
 
z(25) , z(975) = [−0,0981, 0,2875] .

41
 Para probabilidad 0.68, se tiene
 
z(160) , z(840) = [0,0436, 0,2197]

y
cardinal {z ≤ 0} 106
P (pN O − pSI ≤ 0|datos) = = = 0,106.
1000 1000
# Ejemplo
# Se tiene que p1 sigue una Be (101 ,268)
# Se tiene que p2 sigue una Be (1 ,7)

# Simulo dos distribuciones a posteriori

p1 = rbeta (1000 ,101 ,268)
p2 = rbeta (1000 ,1 ,7)

# Calculo la diferencia entre ambas muestras

dif = p1 - p2

# Calculo la probabilidad de los valores que son >=0

difmas0 = dif [ dif >= 0]
( prob = length ( difmas0 ) / 1000)

[1] 0 .893

# Calculo los cuantiles del 5 % y el 95 %

# Esto equivale a un intervalo del 0 .90
quantile ( dif , c (0 .05 ,0 .95 ))

5% 95 %
-0 .09852904 0 .27579033

# Calculo los cuantiles del 2 .5 % y el 97 .5 %

# Esto equivale a un intervalo del 0 .95
quantile ( dif , c (0 .025 ,0 .975 ))

2 .5 % 97 .5 %
-0 .1813385 0 .2886177

# O bien ordeno valores y calculo un intervalo de 0 .90

ordedif = sort ( dif )

ordedif [50]; ordedif [950]

[1] -0 .1007258
[1] 0 .2757859

ordedif [100]; ordedif [900]

[1] -0 .008156479
[1] 0 .2637232

ordedif [150]; ordedif [850]

[1] 0 .03205061
[1] 0 .2512118

42