Muestreador de Gibbs

Descripción e implementación del Muestreador
de Gibbs en versión bivariada
Manuel Alejandro Moreno Arévalo
Universidad Distrital Francisco Jose de Caldas

Facultad de Ciencias y Educación
Proyecto Curricular de Matemáticas
Bogotá, D.C.
Diciembre de 2016
Descripción e implementación del Muestreador
de Gibbs en versión bivariada
Manuel Alejandro Moreno Arévalo
Disertación presentada para optar al tı́tulo de

Matemático
Director
M. Sc. Luis Alejandro Másmela Caita.
Universidad Distrital Francisco Jose de Caldas

Facultad de Ciencias y Educación
Proyecto Curricular de Matemáticas
Bogotá, D.C.
Diciembre de 2016
Tı́tulo en español
Descripción e implementación del Muestreador de Gibbs en versión bivariada
Title in English
Description and implementation of Gibbs Sampler in bivariate version
Resumen: Este trabajo de grado presenta un revisión bibliográfica de los

fundamentos matemáticos del muestreador de Gibbs y una descripción especı́fica de
la mecánica de este. El muestreador de Gibbs es un algoritmo que permite generar
una muestra aleatoria a partir de la distribución de probabilidad conjunta de dos o
más variables aleatorias haciendo uso de sus distribuciones condicionales.
El trabajo se estructuró en tres capı́tulos:

• Capı́tulo 1: Se muestran los fundamentos de la teorı́a de probabilidad en el
caso multivariado, delimitando este estudio a los conceptos empleados en el
muestreador de Gibbs.
• Capı́tulo 2: Presenta la teorı́a de las Cadenas de Markov a tiempo discreto con
espacio de estados discreto. Por ser las cadenas de Markov un soporte teórico
para el muestreador de Gibbs.
• Capı́tulo 3: Expone una descripción del Muestreador de Gibbs en donde se
discute su aplicación y su eficacia.
.
Abstract: This degree work presents a literature review of the mathematical

foundations of the Gibbs sampler and a specified description of the mechanics of
this. The Gibbs sampler is an algorithm that allows to generate a random sample
from the joint probability distribution of two or more random variables making use
of their conditional distributions.
The work was divided into three chapters:

• In chapter 1 the basics of probability theory in the multivariate case are shown,
this study is limited to the concepts used in the Gibbs sampler.
• In chapter 2 it presents the theory of Markov chains with discrete time discrete
state space. The Markov chains being theoretical support for the Gibbs sampler
• In chapter 3 expose a description of the Gibbs sampler where its application
and effectiveness is discussed.
.
Palabras clave: Distribución condicional, Cadenas de Markov, Sucesión de Gibbs.

Keywords: Conditional distribution, Markov chains, Gibbs sequence.
Nota de aceptación
Trabajo de tesis
Jurado
Gabriel Córdoba Suárez
Director
Luis Alejandro Másmela Caita
Bogotá, D.C., Diciembre 7 de 2016

Dedicado a
Mi mamá Rosalia Arévalo y mi abuelita Clementina Arévalo.

Agradecimientos
Quiero agradecer a Dios por brindarme las capacidades para desarrollar mis estudios
y por todas las otras bendiciones que ha derramado sobre mi. A mi mamá Rosalia
Arévalo y mi abuelita Clementina Arévalo quienes han sido mi apoyo y guı́a en las
decisiones que he tomado en mi vida. Quiero agradecer también a los profesores del
Proyecto Curricular de Matemáticas de la Universidad Distrital por sus lecciones
académicas y consejos de vida que han compartido conmigo a lo largo de estos años.
Agradezco especialmente al profesor Luis Alejandro Másmela quien dirigió este trabajo
de grado porque sin sus conocimientos y paciencia no se hubiese podido culminar
satisfactoriamente. Por último, pero no por ello menos importante, al profesor Gabriel
Córdoba, por las correcciones y sugerencias hechas a este documento.
Índice general
Índice general I
Índice de tablas III
Índice de figuras IV
Introducción V
1. Preliminares 1
1.1. Algunos conceptos de la teorı́a de la probabilidad en el caso multivariado 1
1.2. La prueba Ji-cuadrado . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2. Cadenas de Markov 7
2.1. Comunicación . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.2. Periodo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.3. Primeras visitas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.4. Recurrencia y transitoriedad . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.5. Tiempo medio de recurrencia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.6. Clases cerradas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.7. Número de visitas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.8. Recurrencia positiva y nula . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
2.9. Evolución de distribuciones . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.10. Distribuciones estacionarias . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.11. Distribuciones lı́mite . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.12. Cadenas regulares . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.13. Cadenas reversibles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
3. Muestreador de Gibbs 20
3.1. Algoritmo del muestreador de Gibbs para el caso bivariado . . . . . . . . 21
3.2. Ilustrando el muestreador de Gibss . . . . . . . . . . . . . . . . . . . . . . . . . 21
3.2.1. Distribución Beta-binomial . . . . . . . . . . . . . . . . . . . . . . . . . 22
3.2.1.1. Distribución Beta-binomial como distribución com-
puesta . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
I
ÍNDICE GENERAL II
3.2.2. Aplicando el Muestreador de Gibbs . . . . . . . . . . . . . . . . . . . 23

3.2.3. Distribución trinomial . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
3.3. Una prueba de Convergencia simple . . . . . . . . . . . . . . . . . . . . . . . . 29
3.4. Evaluación de P n . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
Conclusiones 38
Códigos en R 39
Bibliografı́a 43
Índice de tablas
3.1. Distribución de probabilidades de (X, Y ) . . . . . . . . . . . . . . . . . . . . 30

3.2. Distribución de probabilidades conjunta de la variable aleatoria dis-
creta (X, Y ) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
3.3. Distribución de probabilidades conjunta para el ejemplo 1. . . . . . . . . 33
III
Índice de figuras
3.1. Comparación de histogramas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

3.2. Comparando el muestreador de Gibbs. . . . . . . . . . . . . . . . . . . . . . . 26
3.3. Comparación de dos Histogramas . . . . . . . . . . . . . . . . . . . . . . . . . . 29
IV
Introducción
Algunos cálculos en la teorı́a de probabilidad requieren de una cierta habilidad

matemática. Por ejemplo, si se tiene la función de densidad de probabilidad conjunta
de dos variables aleatorias y se está interesado en obtener la función de densidad
de una variable aleatoria en particular. Aquı́ presentamos un método conocido
como muestreador de Gibbs que entre sus bondades se encuentra facilitar este tipo
de operaciones. Es precisamente en el Capı́tulo 3 que el lector podrá observar la
mecánica simple y el gran alcance que tiene dicho método. Para tal fin en el citado
capı́tulo se encontrará un algoritmo que describe al muestreador de Gibbs y con
el objetivo de que se observe la bondad del método, este se ha implementado para
dos distribuciones. Se podrá apreciar las casi indistinguibles diferencias que existen
entre las distribuciones teóricas encontradas por métodos analı́ticos y las simuladas
mediante el muestreador de Gibbs. Dichas simulaciones fueron realizadas con el
software estadı́stico R.
El muestreador de Gibbs es una técnica novedosa desarrollada por los hermanos Alan
Geman y Donald Jay Geman en 1984. Al igual que los demás algoritmos conocidos
como MCMC su origen guarda una relación fuerte con la fı́sica y su desarrollo a
la computación. Como si fuera poco es aplicado en campos como las finanzas, la
economı́a y la biologı́a entre otros.
Se ha estructurado el presente trabajo en tres capı́tulos.
En el capı́tulo 1, se exponen conceptos básicos de la teorı́a de la probabilidad en el caso
multivariado, restringiendo el estudio a conceptos que se emplean para desarrollar el
En el capı́tulo 2, se presenta dichos elementos con espacio de estados discreto. Para tal
fin se ha seguido muy de cerca el texto [8]. Este capı́tulo comprende desde la definición
de Cadena de Markov a tiempo discreto con espacio de estado discreto hasta cadenas
reversibles, pasando por la distribución lı́mite que es de vital importancia para el
El capı́tulo 3, comienza presentando una muy breve revisión histórica del muestreador
de Gibbs, expone el algoritmo del muestreador para el caso bivariado, continua con
dos ejemplos para los cuales se realizó en cada caso una simulación empleando
V
INTRODUCCIÓN VI
el muestreador de Gibbs y se contrastanlos resultados con la distribución teórica.

Luego se presentan mediante ejemplos la caracterı́stica markoviana del muestreador.
Culminamos este capı́tulo con la sección llamada evaluación de la matriz P n
CAPÍTULO 1
Preliminares
El objetivo de este capı́tulo es presentar algunos conceptos básicos de la teorı́a de la

probabilidad en el caso multivariado. Al igual que como sucede en el caso univariado
donde a partir de la definición de variable aleatoria se cimientan un gran número
de definiciones, teoremas y/o proposiciones, entre otros, en el caso multivariado la
piedra angular es la definición de vector aleatorio n−dimensional.
En el campo de las aplicaciones es común estar interesado en observar fenómenos
numéricos de dos, tres o más caracterı́sticas simultáneas. Ası́, por analogı́a con
el caso univariado se presenta la definición de función de distribución conjunta
acumulada y algunas de sus propiedades. Se brinda también la definición de densidad
de probabilidad condicional, que es pieza clave para el muestreador de Gibbs.
1.1. Algunos conceptos de la teorı́a de la probabili-

dad en el caso multivariado
Un concepto básico en probabilidad tiene que ver con la definición de vector aleatorio.
Este extiende el concepto de variable aleatoria por una función de valores en Rn .
Definición 1. Sea X1 , X2 , · · · Xn n-variables aleatorias reales definidas sobre el

mismo espacio de probabilidad (Ω, =, P ) . La función X : Ω → Rn definida por
X (w) := (X1 (w) , · · · , Xn (w))
es llamada un n−dimensional vector aleatorio.

Se puede asociar a un vector aleatorio n-dimensional su función de distribución como
lo hace la siguiente definición. La distribución permite calcular la probabilidad de
que el valor del vector aleatorio, pertenezca a un conjunto.
1
CAPÍTULO 1. PRELIMINARES 2
Definición 2. Sea X un n−dimensional vector aleatorio. La medida de probabilidad

definida por:
PX (B) := P (X ∈ B) ; B ∈ Bn
Es llamada la distribución del vector aleatorio X.
Nota 1. Bn denota la σ−álgebra de Borel sobre Rn .
La naturaleza discreta de las variables aleatorias de un vector aleatorio X permite

caracterizar a este como un vector discreto. La siguiente definición evidencia esto e
incluye lo que se entenderá en adelante por función de distribución conjunta.
Definición 3. Sea X = (X1 , X2 , · · · , Xn ) un vector aleatorio n−dimensional. Si

las variables aleatorias Xi , con i = 1, · · · , n son todas discretas, se dice que el
vector aleatorio X es discreto. En este caso la función de masa de probabilidad de
X, también llamada función de distribución conjunta de las variables aleatorias
X1 , X2 , · · · , Xn , esta definida por:
(
P (X = x) si x pertenece a la imagen de X
pX (x) :=
0 en cualquier otro caso.
Es importante añadir que si en la definición anterior las variables aleatorias Xi ,

con i = 1, · · · , n son todas continuas, el vector aleatorio X se dice continuo. Puede
también presentarse el caso en que una variable aleatoria digamos Xj para 1 ≤ j ≤ n
pueda ser continua y las demás discretas o viceversa.
Definición 4. Sea X = (X1 , X2 , · · · , Xn ) un vector aleatorio n−dimensional. La

función definida por:
F (x1 , x2 , · · · , xn ) := P (X1 ≤ x1 , X2 ≤ x2 , · · · , Xn ≤ xn ) ,
para toda (x1 , x2 , · · · , xn ) ∈ Rn es llamada la función de distribución conjunta

acumulada de las variables aleatorias X1 , X2 , · · · , Xn , o simplemente la función de
distribución del vector aleatorio n−dimensional X.
Existe un procedimiento para encontrar la función de distribución marginal acumulada
de una variable aleatoria Xj para 1 ≤ j ≤ n de un vector aleatorio n−dimensional
X, dado este y su función de distribución conjunta acumulada digamos F. Tal
procedimiento es el que establece el siguiente teorema.
Teorema 1. Sea X = (X1 , X2 , · · · , Xn ) un vector aleatorio n−dimensional con

función de distribución conjunta acumulada F. Para cada j = 1, · · · , n, la función
de distribución acumulada de la variable aleatoria Xj está dada por
FXj (x) = lı́m · · · lı́m lı́m · · · lı́m F (x1 , · · · xn ) .

x1 →∞ xj−1 →∞ xj+1 →∞ xn →∞
La función de distribución FXj es llamada la funcion de distribución marginal

acumulada de la variable aleatoria Xj .
Dada la importancia de la función de distribución de vectores aleatorios, es de
interés presentar caracterı́sticas de los mismos. Los Teoremas 2 y 3 ilustran dichas
caracterı́sticas. Por ejemplo, cómo identificar si una función es en efecto función
de distribución conjunta acumulada en el caso bivariado. La demostración de estos
resultados puede consultarse en [4].
Teorema 2. Sea X = (X, Y ) un vector aleatorio 2−dimensional. La función de

distribución conjunta acumulada F de las variables aleatorias X, Y tiene las siguientes
propiedades:
1.
∆ba F := F (b1 , b2 ) + F (a1 , a2 ) − F (a1 , b2 ) − F (b1 , a2 ) ≥ 0
donde a = (a1 , a2 ) , b = (b1 , b2 ) ∈ R2 con a1 ≤ b1 y a2 ≤ b2 .
2.
lı́m F (x, y) = F (x0 , y) (1.1)
x&x0
3.
lı́m F (x, y) = F (x, y0 ) (1.2)
y&y0
4.
lı́m F (x, y) = 0
x→−∞
y
lı́m F (x, y) = 0
y→−∞
5.
lı́m F (x, y) = 1
(x,y)→(∞,∞)
Teorema 3. Sea X = (X1 , X2 , · · · , Xn ) un vector aleatorio n−dimensional. La fun-

ción de distribución conjunta acumulada F de las variables aleatorias X1 , X2 , · · · , Xn
tiene las siguientes propiedades:
1. ∆ba F ≥ 0 para todo a = (a1 , · · · , an ) , b = (b1 , · · · , bn ) ∈ Rn con a ≤ b, donde:
Pn
(−1)( ) F ( a + (1 − ) b , · · · , a + (1 − ) b ) .
X
i=1 i
∆ba F := 1 1 1 1 n n n n
(1 ,···n )∈{0,1}n
2. F es continua a derecha en cada componente.

3. Para todo a1 , · · · , ai−1 , ai+1 , · · · , an ∈ R con i = 1, · · · , n, tenemos:
 
lı́m F a1 , · · · , ai−1 , x , ai+1 , · · · , an  = 0

x&−∞ |{z}
i-ésima posición
4.
lı́m F (x1 , x2 , · · · xn ) = 1
(x1 ,x2 ,···xn )→(∞,∞,··· ,∞)
Definición 5. Sean X1 , X2 , · · · , Xn n variables aleatorias de valor real definidas

sobre el mismo espacio de probabilidad. Se dice que las variables aleatorias son
conjuntamente continuas, si existe una función integrable f : Rn → [0, +∞) tal que
para todo conjunto de Borel C de Rn
Z Z
P ((X1 , X2 , · · · , Xn ) ∈ C) = · · · f (x1 , x2 , · · · , xn ) dx1 dx2 · · · dxn .
C
La función f es llamada la función de densidad de probabilidad conjunta de las

variables aleatorias X1 , X2 , · · · Xn .
Nota 2. Es costumbre notar fdp a la función de densidad de probabilidad.
Los teoremas siguientes exhiben dos relaciones notables el Teorema 4 establece cómo
dado un vector aleatorio 2−dimensional y su función de distribución conjunta se
puede encontrar la función de densidad de probabilidad conjunta en las regiones en
que esta sea continua. La demostración de este hecho se puede encontrar en [4].
Teorema 4. Sean X y Y variables aleatorias continuas con funcion de distribución

conjunta F. Entonces, la función de densidad de probabilidad conjunta f es:
∂ 2 F (x, y) ∂ 2 F (x, y)
f (x, y) = =
∂x∂y ∂y∂x
para todos los puntos (x, y) donde f (x, y) es continua.

Algunas veces dado un conjunto de variables aleatorias, se está interesado en conocer
la función de densidad de una variable en particular. El teorema siguiente establece
los elementos que se deben conocer para llevar a buen fin esta tarea.
Teorema 5. Si X1 , X2 , · · · , Xn son n variables aleatorias de valor real, teniendo f

(función de densidad de probabilidad conjunta), entonces
Z ∞ Z ∞
fXj (x) = ··· f (x1 , · · · , xj−1 , x, xj+1 , · · · , xn ) dx1 · · · dxj−1 dxj+1 · · · dxn
−∞ −∞
es la función de densidad de la variable aleatoria Xj para j = 1, 2, · · · , n.

La Definición 6 y 8 son indispensables para formar la sucesión de Gibbs.
Definición 6. Sean X y Y dos variables aleatorias discretas. La función de probabi-

lidad condicional de X dado Y = y se define como
P (X = x, Y = y)
pX|Y (x|y) := P (X = x|Y = y) = .
P (Y = y)
Para todo y para el cual P (Y = y) > 0.

Definición 7. La función de distribución condicional de X dado Y = y se define
como: X
FX|Y (x|y) := P (X ≤ x|Y = y) = pX|Y (k|y) .
k≤x
Para todo y para el cual P (Y = y) > 0.

La Definición 8 caracteriza la relación entre la función de densidad de una variable
aleatoria y la función de densidad de probabilidad conjunta, llamaremos a esta la
función de probabilidad condicional. En la Definición 6 se presentó para el caso
discreto.
Definición 8. Sean X y Y variables aleatorias continuas con función de densidad

de probabilidad conjunta f. La función de densidad de probabilidad condicional de X
dado Y = y se define como
f (x, y)
fX|Y (x|y) :=
fY (y)
para todo y con fY (y) > 0.
Definición 9. La función de distribución condicional de X dado Y = y se define

como Z x
FX|Y (x|y) := P (X ≤ x|Y = y) = fX|Y (t|y) dt
−∞
Para todo y con fY (y) > 0.
1.2. La prueba Ji-cuadrado
Uno de los objetivos del presente trabajo es exponer al muestreador de Gibbs como
un método optimo para estimar la distribución marginal de una variable aleatoria. Es
ası́ que con los valores obtenidos mediante la sucesión de Gibbs, se quiere presentar
que su comportamiento se asemeja a la distribución teórica. En [5] utilizan las
gráficas para mostrar esto, nosotros además lo presentamos por medio de la prueba
ji-cuadrado implementada en R.
El test ji-cuadrado es un estadı́stico de prueba formulado en 1990 por Karl Pearson.
En este supone que desconocemos la distribución de probabilidad y la hipótesis nula
es que dicha distribución es alguna conocida. Para determinar si la hipótesis se acepta
o rechaza se emplea una muestra. El estadı́stico de prueba es
k
2
X [ni − npi ]2
X =
i=1
npi
Donde n es el tamaño de la muestra y npi son lo valores medios teóricos. Un hecho

que es conocido en estadı́stica pero no es presentado aquı́, pues su tratamiento escapa
de los objetivos del trabajo es que la distribución del estadı́stico X 2 es un distribución

χ2 con k − 1 grados de libertad, cuando la hipótesis nula es cierta.
CAPÍTULO 2
Cadenas de Markov
Debido a la naturaleza Markoviana de la iteración del muestreador de Gibbs se

hace indispensable estudiar las cadenas de Markov que es el propósito del presente
capı́tulo, para hacerlo se ha decidido seguir muy de cerca el texto [8].
Las cadenas Markov datan su aparición entre 1905 a 1907 y se deben al matemático
ruso Andrey Markov considerado por esto como precursor de los procesos estocásticos.
El estudio de cadenas de Markov se extendió ampliamente en las primeras décadas
del siglo XX debido en gran parte a que muchos fenómenos fı́sicos, financieros y
sociales encontraban su representación en este tipo de modelo.
Definición 10. Una cadena de Markov es un proceso estocástico a tiempo discreto

{Xn : n = 0, 1, . . .} , con espacio de estados discreto, y que satisface la propiedad de
Markov, esto es, para cualquier entero n ≥ 0, y para cualesquiera estados x0 , . . . , xn+1 ,
se cumple
p (xn+1 |x0 , . . . , xn ) = p (xn+1 |xn ) (2.1)
La propiedad de Markov establece que la probabilidad condicional de un evento

futuro Xn+1 = k dado el estado presente Xn = i es independiente de cualquier
evento en el pasado. Por otro lado dentro de la teorı́a de las cadenas de Markov
un concepto indispensable es la llamada probabilidad de transición siendo esta
P (Xn = k|Xn−1 = i) la probabilidad de que el proceso pase del estado i al estado
k en un solo paso o unidad de tiempo, se suele denotar por pik . Ası́ mismo se dice
que una cadena de Markov es homogénea o cadena de Markov con probabilidad
estacionaria si la probabilidad de transición es independiente del tiempo en el que se
encuentre la cadena.
7
CAPÍTULO 2. CADENAS DE MARKOV 8
Definición 11. Sea

 
p00 p01 p02 p03 · · ·
p10
 p11 p12 p13 · · ·
P = (pij ) = p20
 p21 p22 p23 · · · (2.2)
p30
 p31 p32 p33 · · ·
.. .. .. ..
. . . .
La matriz para la cual sus elementos son las probabilidades de transición entre los
estados del proceso.
La matriz que se definió se conoce como matriz estocástica o matriz de probabilidad de
transición. La Proposición 1 señala dos propiedades de esta matriz y su demostracion
puede ser consultada en [8].
Proposición 1. La matriz de probabilidad de transición P = (pij ) satisface

1. pij ≥ 0
2. ∞
X
pij = 1 i = 0, 1, 2, 3, . . .
j=0
Definición 12. Una distribución inicial para una cadena de Markov con espacio de
estados {0, 1, 2, 3, . . .} es una distribución de probabilidad sobre este conjunto.
La Proposición 2 recibe su nombre gracias a dos matemáticos de nacionalidades
distintas pero apasionados por los procesos estocásticos, esto les permitió de manera
independiente llegar a su formulación; los matemáticos aquı́ relacionados son Sydney
Chapman y Andréi Kolmogórov de nacionalidades inglesa y rusa respectivamente.
Proposición 2. (Ecuación de Chapman-Kolmogorov) Para cualesquiera en-

teros m, n tales que 0 ≤ m ≤ n, y cualesquiera estados i y j se tiene que
X
pij (n) = pik (m)pkj (n − m). (2.3)
k
La demostración de la Proposición 2 requiere del teorema de probabilidad total y

de la propiedad de Markov. La Proposición 3 es un resultado de la Ecuación de
Chapman-Kolmogorov y ambas demostraciones pueden ser consultadas en [8].
Proposición 3. La probabilidad de transición en n pasos. pij (n), está dada por la

entrada (i, j) de la n-ésima potencia de la matriz P , es decir,
pij (n) = (P n )ij .

2.1. Comunicación
La mayorı́a de textos de procesos estocásticos, afirman que la comunicación es una

relación de equivalencia, la demostración de este hecho se basa en la definición de
pij (0) como la delta Kronecker, es decir,
(
0 si i 6= j
pij (0) = δij =
1 si i = j
y la ecuación de Chapman Kolmogorov. Ası́ pues dada la relación de estados que se

comunican es posible descomponer las cadenas de Markov por los subconjuntos del
espacio de estados que están comunicados.
Definición 13. Se dice que el estado j es accesible desde el estado i si existe un

entero n ≥ 0 tal que pij (n) > 0, esto se escribe simplemente como i → j. Se dice
además que los estados i y j son comunicantes, y se escribe i ↔ j, si se cumple que
i → j y j → i.
Se puede decir que una Cadena de Markov se dice irreducible cuando la relación
de comunicación genera una única clase de equivalencia o como se presenta en la
Definición 14, ambas definiciones son equivalentes.
Definición 14. Se dice que una cadena de Markov es irreducible si todos los estados
se comunican entre sı́.
2.2. Periodo
Definición 15. El periodo de un estado i es un número entero no negativo denotado

por d (i), y definido como sigue:
d (i) = m.c.d. {n ≥ 1 : pii (n) > 0} ,
en donde m.c.d. significa máximo común divisor. Cuando pii (n) = 0 para toda n ≥ 1,
se define d (i) = 0. En particular, se dice que un estado i es aperiódico si d (i) = 1.
Cuando d (i) = k ≥ 2 se dice que i es periódico de periodo k.
Una cadena de Markov en la cual todos sus estados son aperiódicos es llamada una
cadena de Markov aperiódica. El siguiente conjunto de Proposiciones 4, 5 y 6 son
algunas propiedades del periodo. La Proposición 4 se puede interpretar de diversas
maneras, por ejemplo, se puede decir que si dos estados están comunicados entonces
tendrán el mismo periodo lo que sugiere que el periodo es una propiedad de clase.
Proposición 4. Si los estados i y j pertenecen a la misma clase de comunicación,

entonces tienen el mismo periodo.
La Proposición 5 afirma que para cada estado es posible retornar a este con múltiplos
lo suficientemente grandes del periodo d(i).
Proposición 5. Para cada estado i, existe un entero N tal que para toda n ≥ N, se
cumple pii (nd (i)) > 0.
La Proposición 6 es un resultado de las Proposiciones 5 y 2.
Proposición 6. Si pij (m) > 0 para algún entero m, entonces existe un entero N
tal que para toda n ≥ N se cumple pij (m + nd (j)) > 0.
2.3. Primeras visitas
Esta sección trata del tiempo de primera visita de un subconjunto del espacio
de estados de una cadena de Markov. Existen diversas aplicaciones, por ejemplo,
el tiempo de mantenimiento de cierta máquina, es decir, dado que la máquina
funciona correctamente en qué periodo n de tiempo la máquina presentará un
daño y necesitará mantenimiento. Esta aplicación sugiere que se puede relacionar la
Definición 16 con la Definición 17. Es decir, dado de que exista una primera visita al
subconjunto del espacio de estados A :={la máquina presenta averı́a}, con un único
estado al que se puede llamar k y suponiendo que la cadena inicia en el estado j := la
máquina funciona correctamente. Entonces según la Definición 17 fij (n) representa
la probabilidad de que τij tome el valor n. La Proposición 7 es una formula bastante
útil y la demostración requiere de la ecuación 2.1.
Definición 16. Sea A un subconjunto del espacio de estados de una cadena de

Markov {Xn : n ≥ 0} . El tiempo de primera visita al conjunto A es la variable
aleatoria
(
mı́n {n ≥ 1 : Xn ∈ A} si Xn ∈ A para algún n ≥ 1,
τA =
∞ otro caso.
Definición 17. Para cada n ≥ 1, el número fij (n) denota la probabilidad de que una
cadena que inicia en el estado i, llegue al estado j por primera vez en exactamente n
pasos, es decir,
fij (n) = P (Xn = j, Xn−1 6= j, . . . , X1 6= j|X0 = i) .
Adicionalmente se define fij (0) = 0, incluyendo el caso i = j.
Proposición 7. Para cada n ≥ 1,

n
X
pij (n) = fij (k) pij (n − k) .
k=1
2.4. Recurrencia y transitoriedad
En las secciones anteriores se han tratado las probabilidades de transición de un

estado a otro. El objetivo de esta sección es presentar una clasificación de los estados
de una cadena de Markov teniendo en cuenta su eventual probabilidad de regresar al
estado dado que se partió del mismo estado.
Definición 18. Se dice que un estado i es recurrente si la probabilidad de eventual-

mente regresar a i, partiendo de i, es uno, es decir, si
P (Xn = i para alguna n ≥ 1|X0 = i) = 1
Un estado que no es recurrente se llama transitorio, y en tal caso la probabilidad

anterior es estrictamente menor a uno.
Definición 19. Un estado i es recurrente si fii = 1, es decir, si la probabilidad de

regresar a él en un tiempo finito es uno. Análogamente, un estado i es transitorio si
fii < 1.
Proposición 8. El estado i es:

1. recurrente si, y sólo si ∞
P
n=1 pii (n) = ∞.
2. transitorio si, y sólo si ∞

P
n=1 pii (n) < ∞.
La Proposición 8 es empleada en algunos textos de procesos estocásticos como la

definición de estado recurrente y estado transitorio.
Proposición 9. La recurrencia es una propiedad de clase, es decir,

1. Si i es recurrente e i ↔ j, entonces j es recurrente.
2. Si i es transitorio e i ↔ j, entonces j es transitorio.
Proposición
P 10. Sea j un estado transitorio. Para cualquier estado inicial i, se
cumple que ∞ n=1 pij (n) < ∞. En consecuencia, lı́mn→∞ pij (n) = 0.
¿Se puede decir que todos los estados de una cadena de Markov finita son recurrentes?,
de la Proposición 11 se puede extraer respuesta.
Proposición 11. Toda cadena de Markov finita tiene por lo menos un estado
recurrente.
Las pruebas de las Proposiciones 8, 9, 10 y 11, pueden ser consultadas en [8].
2.5. Tiempo medio de recurrencia
Definición 20. El tiempo medio de recurrencia de un estado recurrente j, a partir

del estado i, se define como la esperanza de τij , y se denota por µij , es decir,
∞
X
µij = E (τij ) = nfij (n)
n=1
Notación 1. Se denota τj al tiempo de primera visita cuando el estado de inicio

y de llegada coinciden, para este caso µj denota el tiempo medio de recurrencia y
representa el promedio del número de pasos que la cadena se toma en regresar al
estado j que es recurrente.
2.6. Clases cerradas
En esta sección se considera brevemente el concepto de clase cerrada como una

colección de estados no vacı́a tales que para cualquier estado i que pertenezca a esta
colección y cualquier j que no pertenezca se tiene que pij (n) = 0 para todo n ∈ N.
Definición 21. Una colección de estados no vacı́a C es cerrada si ningún estado

fuera de C es accesible desde algún estado dentro de C , es decir, si para cualquier
i∈C yj∈ / C , i 9 j.
La Proposición 12 establece una forma de caracterizar una clase de comunicación vı́a
colección de estados que es cerrado e irreducible.
Proposición 12. Toda colección de estados que es cerrada e irreducible es una clase
de comunicación.
2.7. Número de visitas
Esta primera parte de la sección concierne a definir una nueva variable aleatoria que
representa el número de visitas que una cadena de Markov hace sobre un estado
particular para un tiempo finito n. La Proposición 13 presenta algunas propiedades
para esta variable aleatoria su demostración puede consultarse en [8].
Definición 22. La variable aleatoria que registra el número de visitas que una cadena
realiza sobre un estado j a partir del estado i, para cualquier tiempo finito n se define
por:
X n
Nij (n) = 1Xk =j .
k=1
Cuando X0 = i.
Notación 2. Se acostumbra a denotar Ni (n) cuando los estados j e i coinciden.
Proposición 13. Para cualesquiera estados i y j,
1. (
1 si k = 0,
P (Nij ≥ k) =
fij (fjj )k−1 si k ≥ 1.
2. (
1 − fij si k = 0,
P (Nij = k) =
fij (fjj )k−1 (1 − fjj ) si k ≥ 1.
3. 
∞
X 0
 si fij = 0,
fij
E (Nij ) = pij (n) = si 0 ≤ fjj < 1,
 1−fjj
n=1
∞ si fij 6= 0 y fjj = 1.

4. (
0 si j es transitorio,
P (Nij = ∞) =
fij si j es recurrente.
5. (
1 si j es transitorio,
P (Nij < ∞) =
1 − fij si j es recurrente.
La siguiente definición permite extender la propiedad de Markov a la Proposición 14.
Definición 23. Una variable aleatoria τ : Ω → {0, 1, . . . ∞} es un tiempo de paro

respecto del proceso estocástico indicado si para cada entero n ≥ 0, el evento (τ = n)
depende únicamente de X0 , X1 , . . . , Xn .
La propiedad fuerte de Markov establece que una cadena de Markov es invariante si
es observada a partir de un tiempo aleatorio τ . La demostración de la Proposición
14 se divide en dos partes verificar que las probabilidades de transición siguen siendo
estacionarias y que se satisface la ecuación 2.4.
Proposición 14. (Propiedad fuerte de Markov) Sea {Xn : n ≥ 0} una cadena

de Markov y sea τ un tiempo de paro respecto de este proceso. Condicionado al
evento (τ < ∞) , el proceso {Xτ +n : n ≥ 0} es una cadena de Markov, es decir, la
probabilidad
P (Xτ +n+1 = j|X0 = x0 , . . . , Xτ +n−1 = xn−1 , Xτ +n = i) (2.4)
Es igual a P (Xτ +n+1 = j|Xτ +n = i)

Para finalizar esta sección se presenta el teorema ergódico para cadenas de Markov
junto con un corolario. El teorema ergódico proporciona el comportamiento lı́mite
del promedio temporal del número de vistas a un estado de una cadena de Markov,
tanto el teorema ergódico como su corolario hacen parte de los teoremas lı́mites para
cadenas de Markov. La demostración de estos resultados puede ser consultada en [8].
Teorema 6. (Teorema ergódico para cadenas de Markov) Para cualesquiera

estados i y j de una cadena de Markov irreducible se cumple que
Nij (n) 1
lı́m = c.s. (2.5)
n→∞ n µj
siendo este lı́mite cero cuando µj = ∞.
Notación 3. c.s.:= casi seguramente.
Corolario 1. Si una cadena es irreducible y recurrente, entonces

n
1X 1
lı́m Pij (k) =
n→∞ n µj
k=1
2.8. Recurrencia positiva y nula
En la Sección 2.4 se presentó una primera división para toda cadena de Markov por
su espacio de estados en dos subconjuntos disjuntos a saber recurrentes y transitorios.
El objetivo de esta sección es realizar una nueva división del espacio de estados
empleando la Definición 20 y responder la pregunta ¿Existen estados recurrentes
nulos en Cadenas de Markov finitas? Para lo cual es conveniente ver la Proposición
16.
Definición 24. Se dice que un estado recurrente i es:

1. recurrente positivo si µi < ∞.
2. recurrente nulo si µi = ∞.
Proposición 15. Sea i un estado recurrente. Entonces,

1. si i es recurrente positivo e i ↔ j, entonces j es recurrente positivo.
2. si i es recurrente nulo e i ↔ j, entonces j es recurrente nulo.
Proposición 16. No existen estados recurrentes nulos en cadenas de Markov finitas.

2.9. Evolución de distribuciones
Esta sección es una adaptación de [8] de la sección del mismo nombre. El propósito es
describir cómo se puede obtener una sucesión infinita de distribuciones de probabilidad
teniendo como partida un espacio de estados finito {0, 1, . . . , n} , una distribución
inicial π 0 = (π00 , π10 , . . . , πn0 ) . y una matriz de probabilidad de transición en un
paso. La cadena se encuentre en alguno de sus posibles estados de acuerdo a la
distribución π 1 = (π01 , π11 , . . . , πn1 ) , una vez pasada la primera unidad de tiempo, para
esta distribución la j−ésima entrada de este vector es
πj1 = P (X1 = j)
Xn
= P (X1 = j|X0 = i) P (X0 = i)
i=0
n
X
= πi0 pij
i=0
Esto evidencia que mediante la multiplicación π 1 = π 0 P, se puede obtener el vector

π 1 vı́a el vector π 0 y la matriz de probabilidades de transición P , es decir
 
p00 · · · p0n
π01 , π11 , . . . , πn1 = π00 , π10 , . . . , πn0  ... .. 

. 
pn0 · · · pnn
De manera análoga se puede obtener el vector π 2 mediante π 2 = π 1 P = π 0 P 2 , si se

continuara con este procedimiento, se llegarı́a que para m ≥ 1,
π m = π m−1 P = π 0 P m (2.6)
Ası́ se obtendrı́a sucesión infinita de distribuciones de probabilidad π 0 , π 1 , π 2 , . . . ,

en donde a excepción de la primera, todas se obtuvieron de su inmediata anterior
multiplicando a derecha por la matriz de probabilidades de transición en un paso.
2.10. Distribuciones estacionarias
Definición 25. Una distribución de probabilidad π = (π0 , π1 , . . .) es estacionaria o

invariante para una cadena de Markov con matriz de probabilidades de transición
P = (pij ) si X
πj = πi pij .
i
Al igual que sucede en álgebra lineal donde un sistema de ecuaciones lineales puede
tener infinitas soluciones o única solución o no existir solución una distribución
estacionaria para una cadena de Markov puede también presentarse de manera única,
múltiple o no existir. En libros de procesos estocásticos pueden encontrarse un buen
número de ejemplos.
Por lo dicho anteriormente es óptimo enunciar una definición alterna de distribución
estacionaria a saber: una distribución de probabilidad π = (π0 , π1 , . . .) es estacionaria
o invariante para una cadena de Markov con matriz de probabilidades de transición
P = (pij ) si,
π = πP. (2.7)
Producto de esta definición alterna se facilita obtener algunos resultados, por ejemplo,
si se tuviera una matriz de transición P para una cadena de Markov y se quisiera
conseguir una probable distribución P estacionaria se podrı́a resolver el sistema de
ecuaciones π = πP con la condición j πj = 1. una interesante propiedad matemática
de las distribuciones estacionaria es que forma un conjunto convexo.
Se conoce como distribución estacionaria debido a la propiedad de ser invariante por
el paso del tiempo, esto es una consecuencia de la ecuación 2.7 y se enuncia como:
una distribución π estacionaria satisface que para todo m ∈ N se tiene que π = πP m .
Es preciso destacar que aunque el vector nulo o vector de ceros cumple la ecuación
2.7 este no es una distribución de probabilidad.
Proposición 17. (Soporte de una distribución estacionaria) Sea π una

distribución estacionaria para una cadena de Markov. Si j es un estado transitorio o
recurrente nulo, entonces πj = 0.
La demostración de las Proposiciones 17 y 18 pueden ser consultadas en [8]. La
importancia de la Proposición 18 es que concibe una fórmula de la distribución
estacionaria para una cadena de Markov irreducible y recurrente positiva, tomando
el inverso del tiempo medio de recurrencia para el correspondiente estado. Como
también establece que toda cadena finita e irreducible tiene una única distribución
estacionaria.
Proposición 18. (Existencia y unicidad de la distribución estacionaria)

Toda cadena de Markov que es irreducible y recurrente positiva tiene una única
distribución estacionaria dada por
1
πj = > 0,
µj
en donde µj es el tiempo medio de recurrencia del estado j. En particular, toda

cadena finita e irreducible tiene una única distribución estacionaria.
2.11. Distribuciones lı́mite
Definición 26. Considere una cadena de Markov con matriz de probabiliades de

transición P = (pij ) y distribución inicial π 0 . Se le llama distribución lı́mite de esta
cadena al vector X
π = lı́m π 0 P n = lı́m πi0 pij (n)
n→∞ n→∞
i
La definición de distribución lı́mite para una cadena de Markov establece que esta se
determina por el lı́mite de las potencia de P y por tanto no depende de la distribución
inicial. Además es una distribución estacionaria.
No siempre una distribución estacionaria es una distribución lı́mite, es natural
entonces preguntarse qué condiciones debe satisfacer una distribución lı́mite para
ser estacionaria, para responder esta pregunta se presenta la Proposición 19 que es
válida tanto para espacios finitos como infinitos.
Obsérvese también que la Proposición 19 no hace la salvedad en que los lı́mites
puedan ser todos nulos pero finaliza recalcando que se obtiene una distribución de
probabilidad verdadera cuando el espacio de estados es finito. La importancia de
este argumento es que existen vectores que satisfacen la Definición 26 pero no son
distribuciones de probabilidad.
Proposición 19. Considere una cadena de Markov con probabilidades de transición

pij tales que los lı́mites πj = lı́mn→∞ pij (n) existen para cada j, y no dependen del
estado i. Entonces
1. X
πj ≤ 1.
j
2. X
πj = πi pij .
i
Cuando el espacio de estados es finito se cumple la igualdad en el primer resultado

obteniéndose una distribución de probabilidad verdadera.
La Proposición 20 amplia el conocimiento de la distribución lı́mite pues enuncia
qué condiciones se deben tener para la existencia de los lı́mites de probabilidad. Suele
pensarse debido a la condición 3 que la Proposición 20 aparenta ser un reciproco de
la Proposición 19.
Proposición 20. (Convergencia a la distribución estacionaria) Considere

una cadena de Markov que es:
1. irreducible,
2. aperiódica, y
3. con distribución estacionaria π

entonces para cualesquiera estados i y j, lı́mn→∞ pij (n) = πj .
La Proposición 21 juega un papel importante dentro de las cadenas de Markov, pues
establece las condiciones para la existencia de probabilidades limite y suministra
además una forma explı́cita para hallarlas. Por otro lado es un argumento que
relaciona una buena cantidad de elementos matemáticos de las cadenas de Markov
y fruto de esto su demostración es bastante corta. La prueba de las proposiciones
presentadas en esta sección puede consultarse en [8].
Proposición 21. (Convergencia para cadenas de Markov) Considere una

cadena de Markov que es:
1. irreducible,
2. recurrente positiva, y
3. aperiódica.
Entonces las probabilidades lı́mite πj = lı́mn→∞ pij (n) existen, están dadas por
πj = 1/µj , y constituyen la única solución al sistema de ecuaciones
X
πj = πi pij (2.8)
i
P
sujeto a las condiciones πj ≥ 0, y j πj = 1.
2.12. Cadenas regulares
Existen cadenas de Markov para las cuales alguna potencia n de la correspondiente

matriz de transición tiene todos sus elementos no nulos, es decir, P n es una matriz
con todas sus entrada mayores que cero para algún n ∈ N. A este tipo de cadenas
de Markov se le conocen como cadenas regulares; un método no muy eficiente para
conocer si una cadena es regular consiste en calcular las potencias de la matriz de
transición asociada a la cadena hasta encontrar un n tal que P n tenga todas sus
entradas estrictamente positivas.
La Proposición 22 es equivalente a la Definición 27. La Proposición 23 enuncia el

comportamiento lı́mite de algunas cadenas finitas.
Definición 27. Se dice que una cadena finita o su matriz de probabilidades de

transición es regular si existe un natural n tal que pij (n) > 0, para cualesquiera
estados i y j.
Proposición 22. Una matriz estocástica es regular si, y sólo si, es finita, irreducible
y aperiódica.
Proposición 23. Toda cadena finita que además es:

1. regular, tiene como distribución lı́mite la única solución no negativa del sistema
de ecuaciones 2.8.
2. regular y doblemente estocástica, tiene como distribución lı́mite la distribución
uniforme.
3. irreducible, aperiódica y doblemente estocástica, tiene como distribución lı́mite
la distribución uniforme.
2.13. Cadenas reversibles
Se cierra este capı́tulo con un muy breve estudio de un tipo particular de cadenas
irreducibles conocidas como cadenas reversibles. La condición de reversibilidad suele
describirse afirmando que para cualesquiera dos estados j e i de una cadena de
Markov la tasa de transición de i a j coincide con la tasa de transición de j a i. La
Proposición 24 describe como caracterizar una cadena de Markov reversible para la
cual exista una distribución que satisfaga la identidad 2.9.
Definición 28. Se dice que una cadena de Markov irreducible con probabilidades
de transición pij y con distribución estacionaria π es reversible en el tiempo si para
cualesquiera estados i y j,
πi pij = πj pji . (2.9)
Proposición 24. Considere una cadena irreducible para la cual existe una distribu-
ción π que cumple la identidad 2.9. Entonces π es una distribución estacionaria y la
cadena es reversible.
CAPÍTULO 3
Muestreador de Gibbs
El muestreador de Gibbs es uno de los algoritmos markovianos, estos son conocidos

como métodos MCMC (proviene del inglés Monte Carlo Markov Chain). Los cuales se
basan en emplear el método de Monte Carlo vı́a cadenas de Markov. A continuación,
presentaremos un muy breve desarrollo histórico de los métodos MCMC. Para esto,
es conveniente iniciar por el método de Montecarlo, ya que como se ha dicho es la
raı́z de los algoritmos MCMC. El método de Montecarlo surge entre 1944 y 1946, su
descubrimiento lo comparten los matemáticos Stanislaw Ulam y John von Neumann.
Es preciso destacar que la idea original del método se debe a Ulam, pero es Neumann
quien realiza mejoras y un aporte significativo para que el método de Monte Carlo
se hubiese implementado en 1948, para obtener los valores singulares de la ecuación
de Schrödinger aunque también se recibió aportes del fı́sico Enrico Fermi.
Más tarde, en 1953 el fı́sico-matemático Nicholas Constantine Metrópolis propone el
primer algoritmo que lleva su nombre, para generar muestras de la distribución de
Boltzmann. Dicho algoritmo es generalizado en 1970 por el estadı́stico canadiense
W.K. Hastings
Vale la pena señalar, la estrecha relación que ha existido entre la fı́sica y los algoritmos
MCMC. Del cual el muestreador de Gibbs no es la excepción, pues surge en el
campo de la fı́sica estadı́stica con el propósito de poder simular campos aleatorios
markovianos, en especial los campos aleatorios de Gibbs.
El muestreador de Gibbs fue propuesto en 1984 por los hermanos y matemáticos
Stuart Alan Geman y Donald Jay Geman en el paper “Stochastic Relaxation, Gibbs
Distribution, and the Bayesian Restoration of Images”. Pero, es hasta en 1990
con la publicación de Alan E. Gelfan y Adrian F.M. Smith de “Sampling-Based
Approaches to Calculating Marginal Densities ”y por sus conferencias impartidas que
el muestreador de Gibbs se considera como una revelación en estadı́stica. Además
debido a los recientes avances de la computación el muestreador de Gibbs se aplica
ampliamente en campos como las finanzas, la economı́a y la biologı́a entre otros.
20
CAPÍTULO 3. MUESTREADOR DE GIBBS 21
3.1. Algoritmo del muestreador de Gibbs para el

caso bivariado
Esta sección es un primer acercamiento a la mecánica del muestreador de Gibbs. Pues

aquı́ se presenta un algoritmo para implementarlo en el caso bivariado. El objetivo
de dicho algoritmo es que el lector comprenda cómo surge la sucesión de Gibbs para
el caso que estamos considerando y que observe también, el mecanismo sencillo que
describe al muestreador de Gibbs.
Algoritmo 1.
Entrada: • Distribuciones condicionales, fX|Y (x | y) y fY |X (y | x) .
• Número máximo de iteraciones k.
0 0
• Valor inicial de Y0 = y0 .
Salida: Sucesión de Gibbs
0 0
Paso 1. Sea i = 0, Y0 = y0 .
Paso 2. Mientras i + 1 ≤ k realizar pasos 3 y 4.
0 0 0
Paso 3. Obtener Xi mediante f (x | Yi = yi ).
0 0 0
Paso 4. Obtener Yi+1 mediante f (y | Xi = xi ).
0 0
• Salida Xi , Yi+1 (Algoritmo culminado satisfactoriamente).
• PARAR.
Paso 5. Salida (Se pidió la implementación del método como máximo para k
iteraciones).
• PARAR.
Los valores
0 0 0 0 0 0 0 0
Y0 , X0 , Y1 , X1 , Y2 , X2 , . . . , Yk , Xk , (3.1)
obtenidos al aplicar el Algoritmo 1 se conocen como “Sucesión de Gibbs”.
El éxito del muestreador de Gibss se basa en que haciendo k lo suficientemente
0 0
grande, la observación de Xk = xk es un punto muestral de f (x), es decir, cuando
0
k → ∞ la distribución de Xk → f (x) la verdadera marginal de X.
3.2. Ilustrando el muestreador de Gibss
Ahora ya familiarizados con el muestreador de Gibbs podemos implementarlo. Esta

sección lo hace, su objetivo es comparar los resultados obtenidos mediante métodos
analı́ticos y el muestreador de Gibbs. Para esto se consideran dos distribuciones con-
juntas, a saber, las que establecen las ecuaciones 3.6 y 3.9. Para dichas distribuciones
se encuentra por medios analı́ticos sus distribuciones marginales y condicionales.

Luego se emplea el muestreador de Gibbs para simular la distribución marginal de
X, para las distribuciones en 3.6 y 3.9. Por último se comparan los histogramas
entre la distribución teórica y la simulada. Esto se desarrolla en la subsecciones
3.2.2 y 3.2.3. Pues en la primera se quiso realizar una introducción a la distribución
marginal que se simulará mediante Gibbs en la subsección 3.2.2. Además, porque
dicha distribución es frecuente en la estadı́stica Bayesiana y se presenta una forma
para obtenerla con métodos de estadı́stica clásica.
3.2.1. Distribución Beta-binomial
En estadı́stica clásica los parámetros de una distribución son siempre constantes.

Ahora bien, en estadı́stica Bayesiana se consideran distribuciones cuyos parámetros
no son constantes sino aleatorios. Un ejemplo, de tal situación es la familia de
distribuciones discretas de probabilidad conocida como Beta-binomial. Esta se obtiene
considerando que la probabilidad de éxito p, en una distribución Bin(n, p), por cada
ensayo no es fija sino que se distribuye Beta(α, β) .
Notación 4. La función Beta es,
Z 1
B (α, β) = y α−1 (1 − y)β−1 dy.
0
La notaremos por B (α, β) .
3.2.1.1. Distribución Beta-binomial como distribución compuesta
Sea X una variable aleatoria con parámetro P que se distribuye Beta(α, β). Sabemos
que,
n k
pX|P (x | p) = p (1 − p)n−x , (3.2)
x
entonces
pα−1 (1 − p)β−1

n x
f (p, x) = pX|P (x | p) = p (1 − p)n−x . (3.3)
x B (α, β)
Ası́ que la distribución compuesta está dada por,
pα−1 (1 − p)β−1
Z 1
n x
pX (x) = p (1 − p)n−x dp
0 x B (α, β)
Z 1
n 1
= px (1 − p)n−x pα−1 (1 − p)β−1 dp
x B (α, β) 0
(3.4)
px+α−1 (1 − p)n−x+β−1
Z 1
n 1
= B (x + α, n − x + β) dp
x B (α, β) 0 B (x + α, n − x + β)

n B (x + α, n − x + β)
= .
x B (α, β)
Empleando la identidad B (α, β) = Γ(α)Γ(β)

Γ(α+β)
, se puede reescribir la ecuación 3.4 como
sigue,

n B (x + α, n − x + β)
pX (x) =
x B (α, β)
Γ(x+α)Γ(n−x+β)
n Γ(n+α+β)
= Γ(α)Γ(β)
(3.5)
x
Γ(α+β)

n Γ (α + β) Γ (x + α) Γ (n − x + β)
=
x Γ (α) Γ (β) Γ (α + β + n)
la ecuación 3.5, corresponde a la función de probabilidad de la distribución beta-

binomial.
3.2.2. Aplicando el Muestreador de Gibbs
Se busca exponer que el muestreador de Gibbs es un método óptimo. Por ejemplo,

para generar una muy buena aproximación a la densidad marginal f (x) de una
distribución conjunta f (x, y). Dicho lo anterior, dada una función de distribución
conjunta se obtendrá de manera analı́tica y haciendo uso del muesteador de Gibbs la
distribución marginal para una variable aleatoria distribuida beta-binomial.
A continuación, se presenta el desarrollo analı́tico para las variables aleatorias X y
Y con función de probabilidad conjunta dada por:

n Γ (α + β) x+α−1
fX,Y (x, y) = y (1 − y)n−x+β−1 x = 0, 1, 2, 3, . . . , n 0 ≤ y ≤ 1.
x Γ (α) Γ (β)
(3.6)
Con respecto a 3.6, es posible encontrar las distribuciones fX (x) , fY (y) , fX|Y (x, y)
y fY |X (y|x) como sigue

Z 1
n Γ (α + β) x+α−1
fX (x) = y (1 − y)n−x+β−1 dy
x Γ (α) Γ (β)
0 Z 1
n Γ (α + β)
= y x+α−1 (1 − y)n−x+β−1 dy (3.7)
x Γ (α) Γ (β) 0

n Γ (α + β) Γ (x + α) Γ (n − x + β)
= .
x Γ (α) Γ (β) Γ (α + n + β)
Ası́ que fX (x) es beta-binomial de parámetros n, α y β. A continuación,

n
X n Γ (α + β) x+α−1
fY (y) = y (1 − y)n−x+β−1
x=0
x Γ (α) Γ (β)
n
Γ (α + β) α−1 β−1
X n x
= y (1 − y) y (1 − y)n−x
Γ (α) Γ (β) x=0
x
Γ (α + β) α−1
= y (1 − y)β−1 .
Γ (α) Γ (β)
Es decir, fY (y) es Beta (α, β) . Luego,

n Γ(α+β) x+α−1
(1 − y)n−x+β−1

x Γ(α)Γ(β)
y
fY |X (y | x) = n Γ(α+β) Γ(x+α)Γ(n−x+β)

x Γ(α)Γ(β) Γ(α+n+β)
Γ (α + n + β)
= y x+α−1 (1 − y)n−x+β−1
Γ (x + α) Γ (n − x + β)
−1
Γ (x + α) Γ (n − x + β)
= y x+α−1 (1 − y)n−x+β−1
Γ (n + α + β)
1
= y x+α−1 (1 − y)n−x+β−1 ,
B (x + α, n − x + β)
Γ(α)Γ(β)
donde se hizo uso de B (α, β) = Γ(α+β)
, para obtener que fY |X (y | x) es Beta(x +
α, n − x + β). Por último,
n Γ(α+β) x+α−1
(1 − y)n−x+β−1

x Γ(α)Γ(β)
y
fX|Y (x | y) = Γ(α+β) α−1
Γ(α)Γ(β)
y (1 − y)β−1

n x
= y (1 − y)n−x ,
x
es decir fX|Y (x | y) es Binomial(n, y).

0
Hemos realizado una simulación en R con k = 100000 iteraciones, de valor inicial Y0 =
0.4, para la distribución marginal de X de la distribución conjunta 3.6, implementando
el muestreador de Gibbs. Esto con el fin de contrastar los resultados analı́ticos y el
muestreador de Gibbs. Ahora, el lector puede apreciar el objetivo de haber encontrado

por metodos analı́ticos la ecuación 3.7. Pues esta indica que la distribución marginal
de X es beta-binomial y que la gráfica producida por la marginal simulada se debe
parecer al comportamiento de dicha distribución. En efecto, esto lo muestra la
figura 3.1, donde el histograma de color azul representa el obtenido empleando el
muestreador de Gibbs y el de color rojo se obtuvo directamente de la distribución
beta-binomial con parámetros n = 10, α = 2 y β = 5. El lector puede apreciar
también que ambos histogramas son similares, hecho que evidencia lo eficaz del
muestreador de Gibss.
Figura 3.1. Comparación de histogramas
Hemos también querido implementar la prueba ji-cuadrado. Con el fin de comparar

los valores obtenidos mediante la sucesión de Gibbs y los teóricos. El tamaño de la
muestra es 100000, se construyeron k = 11 celdas. La hipótesis nula es H0 : Los datos
provienen de una distribución Beta-binomial con parámetros n = 10, α = 2 y β = 5.
La hipótesis alterna es H1 : Los datos no provienen de una distribución Beta-binomial
con parámetros n = 10, α = 2 y β = 5. Con el código 3.4 se obtuvó un p-valor de
1 y un valor χ2 de 0.0001129436 dado que este es menor que χ20.05 = 18.3070 para
10 grados de libertad. Se acepta H0 , es decir, no hay evidencia para decir que estos
datos no provienen de una distribución Beta-binomial con parámetros n = 10, α = 2
y β = 5.
También se ha simulado la distribución beta-binomial con el muestreador de Gibbs,
con el ánimo de observar que para valores pequeños la distribución simulada se asemeja
rápidamente a la teórica. La figura 3.2 permite apreciar esto. Los histogramas en
azul y verde representan la distribución simuladas empleando Gibbs con k = 1000 y
k = 10000 respectivamente y el histograma de color rojo la distribución teórica, es
decir, la distribución beta-binomial con parámetros n = 10, α = 2 y β = 5.
Figura 3.2. Comparando el muestreador de Gibbs.
3.2.3. Distribución trinomial
La distribución trinomial es un caso particular de la distribución multinomial, que

es la análoga de la distribución binomial en el caso multivariado. La distribución
trinomial cuantifica la probabilidad de que en n ensayos independientes se obtenga
precisamente xi veces un evento Ei , donde i = 1, 2, 3. Con la condición de que en
un único ensayo se obtenga uno de los Ei , Además, pi denota la correspondiente
probabilidad. La distribución trinomial tiene función de probabilidad dada por
(
n!
p x1 p x2 p x3 si 3i=1 xi = n,
P
x1 !x2 !x3 ! 1 2 3
f (x1 , x2 , x3 ; n, p1 , p2 , p3 ) = (3.8)
0 en otro caso.
También 3.8 satisface 3i=1 pi = 1, de donde se sigue que p3 = 1 − p1 − p2 . Por otro

P
lado x3 = n − x1 − x2 . Ası́ que reemplazando estos valores en 3.8, se logra suprimir
la dependencia de x3 y p3 , es decir, podemos reescribir 3.8 como sige,
n!
f (x1 , x2 ; n, p1 , p2 ) = px1 px2 (1 − p1 − p2 )n−x1 −x2 . (3.9)
x1 !x2 ! (n − x1 − x2 )! 1 2
Encontremos ahora la distribución marginal de X1 y X2 .

n−x
X1 n!
fX1 (x1 ) = px1 px2 (1 − p1 − p2 )n−x1 −x2
x2 =0
x1 !x2 ! (n − x1 − x2 )! 1 2
n−x
(3.10)
px1 1 X1 n!
= px2 (1 − p1 − p2 )n−x1 −x2
x1 ! x =0 x2 ! (n − x1 − x2 )! 2
2
antes de realizar el siguiente paso observe que x1 + x2 + x3 = n, de donde n ≥ x1 ,

ası́ que
n! = (1)(2)(3) · · · (n − x1 )(n − x1 + 1)(n − x1 + 2) · · · (n)
| {z }
(n−x1 )!
regresando a 3.10
n−x
px1 1 X1 (n − x1 )!
= (n − x1 + 1)(n − x1 + 1) · · · (n) px2 2 (1 − p1 − p2 )n−x1 −x2
x1 ! x2 =0 2
x ! (n − x 1 − x 2 )!
n−x
px1 1 X1 n − x1

= (n − x1 + 1) · · · (n) px2 2 (1 − p1 − p2 )n−x1 −x2
x1 ! x =0
x2
2
px1 1
(n − x1 )!(n − x1 + 1) · · · (n)
= (p2 + (1 − p1 − p2 ))n−x1
x1 ! (n − x1 )!
n!
= px1 1 (1 − p1 )n−x1
x ! (n − x1 )!
1
n x1
= p (1 − p1 )n−x1 .
x1 1
(3.11)
Es ası́ que fX1 (x1 ) es Binomial(n, p1 ). En seguida para fX2 (x2 ), se tiene que
n−x
X2 n!
fX2 (x2 ) = px1 1 px2 2 (1 − p1 − p2 )n−x1 −x2
x1 =0
x1 !x2 ! (n − x1 − x2 )!
n−x
px2 2 n! X2 (n − x2 )!
= px1 (1 − p1 − p2 )n−x1 −x2
x2 ! x =0 (n − x2 )!x1 !(n − x2 − x1 )! 1
1
n−x
px2 2 n! X2 (n − x2 )!
= px1 (1 − p1 − p2 )n−x2 −x1
(n − x2 )!x2 ! x =0 x1 !(n − x2 − x1 )! 1
1
n−x 2
n x2 X n − x 2 x1
= p2 p1 (1 − p1 − p2 )n−x2 −x1
x2 x1 =0
x1

n x2
= p (p1 + (1 − p1 − p2 ))n−x2
x2 2

n x2
= p (1 − p2 )n−x2 .
x2 2
De donde fX2 (x2 ) es Binomial(n, p2 ). Encontremos a continuación la distribución

condicional, fX1 |X2 (x1 | x2 ),
n!
px1 px2 (1 − p1 − p2 )n−x1 −x2
x1 !x2 !(n−x1 −x2 )! 1 2
fX1 |X2 (x1 | x2 ) = n!
px2 (1 − p2 )n−x2
x2 !(n−x2 )! 2
n!x2 ! (n − x2 )! px2 2 px1 1 (1 − p1 − p2 )n−x1 −x2 (1 − p2 )x1

=
n!x1 !x2 ! (n − x1 − x2 )! px2 2 (1 − p2 )n−x2 (1 − p2 )x1
x1 n−x1 −x2
(n − x2 )! p1 (1 − p1 − p2 )
= x1
x1 ! (n − x2 − x1 !) (1 − p2 ) (1 − p2 )−x1 (1 − p2 )n−x2
x1 n−x2 −x1
n − x2 p1 1 − p1 − p 2
=
x1 1 − p2 1 − p2
x1 n−x2 −x1
n − x2 p1 1 − p2 p1
= −
x1 1 − p2 1 − p2 1 − p2
x1 n−x2 −x1
n − x2 p1 p1
= 1− .
x1 1 − p2 1 − p2
p1
Por lo tanto fX1 |X2 (x1 | x2 ) se distribuye binomial de parametros n − x2 y 1−p2
. Para
la distribución condicional, fX2 |X1 (x2 | x1 ), se tiene que,
n!
px1 px2 (1 − p1 − p2 )n−x1 −x2
x1 !x2 !(n−x1 −x2 )! 1 2
fX2 |X1 (x2 | x1 ) = n!
px1 (1 − p1 )n−x1
x1 !(n−x1 )! 1
n!x1 ! (n − x1 )! px1 1 px2 2 (1 − p1 − p2 )n−x1 −x2 (1 − p1 )x2

=
n!x1 !x2 ! (n − x1 − x2 )! px1 1 (1 − p1 )n−x1 (1 − p1 )x2
x2 n−x1 −x2
(n − x1 )! p2 (1 − p1 − p2 )
= x2
x2 ! (n − x1 − x2 )! (1 − p1 ) (1 − p1 )n−x1 (1 − p1 )−x2
x2 n−x1 −x2
n − x1 p2 1 − p 1 − p2
=
x2 1 − p1 1 − p1
x2 n−x1 −x2
n − x1 p2 1 − p1 p2
= −
x2 1 − p1 1 − p 1 1 − p1
x2 n−x1 −x2
n − x1 p2 p2
= 1− .
x2 1 − p1 1 − p1
p2
Es decir, fX2 |X1 (x2 | x1 ) se distribuye binomial de parametros n − x1 y 1−p1
.
La figura 3.3 contrasta la distribución teórica en 3.11 y la simulada empleando
el muestreador de Gibbs, que ha sido creada en R con k=200000 iteraciones. En
dicha figura puede apreciarse lo casi indistinguible que es el comportamiento de
ambas distribuciones, en donde el histograma en color azul representa la distribución
simulada mediante Gibbs y el histograma en color rojo la distribución teórica en

3.11.
Figura 3.3. Comparación de dos Histogramas
A continuación se presenta la prueba ji-cuadrado con el fin de comparar los valores

obtenidos mediante la sucesión de Gibbs y los teóricos. El tamaño de la muestra es
200000, se construyeron k = 11 celdas. La hipótesis nula es H0 : Los datos provienen
de una distribución trinomial con parámetros n = 10 p1 = 0.3, p2 = 0.2 y p3 = 0.5.
La hipótesis alterna es H1 : Los datos no provienen de una distribución trinomial
con parámetros n = 10 p1 = 0.3, p2 = 0.2 y p3 = 0.5. En R se obtuvó un p-valor de
0.9999997 y un valor χ2 de 0.2569682 como este valor es menor que χ20.05 = 18.3070
con 10 grados de libertad. Se acepta H0 , es decir, no hay evidencia para decir que
estos datos no provienen de una distribución distribución trinomial con parámetros
n = 10 p1 = 0.3, p2 = 0.2 y p3 = 0.5
3.3. Una prueba de Convergencia simple
En la sección anterior se obtuvo de manera analı́tica y por medio del muestreador

de Gibbs la distribución marginal para la variable aleatoria X de dos distribuciones
conjuntas f (X, Y ). Ahora, esta nueva sección presenta una matriz que notaremos
por AX|X y posee la propiedad de que su distribución estacionaria es fX .Si bien
su desarrollo se hace en el caso discreto también presentamos un paralelo entre la
versión discreta y la continua.
El objetivo de esta sección es que el lector se percate de la naturaleza markoviana
del muestreador de Gibbs, para ello se han hecho comentarios y consideraciones que
permitan apreciarlo.
Para comenzar consideremos primero un ejemplo, Sean X y Y dos variables aleatorias
cada una Bernoulli, con función de distribución de probabilidad conjunta como lo
indica la tabla siguiente:
HH
X
H 0 1 Suma
Y HH
H
0 0.2 0.1 0.3
1 0.4 0.3 0.7
Suma 0.6 0.4 1
Tabla 3.1. Distribución de probabilidades de (X, Y ).
podemos expresar la tabla 3.1 en términos de la función de probabilidad conjunta, es

decir,
fX,Y (0, 0) fX,Y (1, 0) 0.2 0.1
=
fX,Y (0, 1) fX,Y (1, 1) 0.4 0.3
ası́ que las distribuciones marginales de X y Y están dadas por
fX = [fX (0), fX (1)] = [0.6, 0.4]
fY = [fY (0), fY (1)] = [0.3, 0.7] ,

que se distribuyen Bernoulli con probabilidad de éxito 0.4 y 0.7. Luego podemos
calcular la distribución marginal de X | Y = y,
P (0, 0) 0.2
fX|Y (0 | 0) = =
P (0) 0.3
P (1, 0) 0.1
fX|Y (1 | 0) = =
P (0) 0.3
P (0, 1) 0.4
fX|Y (0 | 1) = =
P (1) 0.7
P (1, 1) 0.3
fX|Y (1 | 1) = = .
P (1) 0.7
Para construir una matriz que contenga las probabilidades condicionales de X dado
Y =y 0.2 0.1
AX|Y = 0.40.3 0.3
0.3
0.7 0.7
y de manera análoga podemos calcular la distribución marginal de Y | X = x, y

obtener la siguiente matriz 0.2 0.4
0.6
AY |X = 0.1 0.6 .
0.3
0.4 0.4
Si consideramos el espacio de estados {0, 1}, es natural pensar en las matrices

AX|Y y AY |X como matrices de transición. Donde por ejemplo AX|Y contiene las
probabilidades de ir del estado X al estado Y y AY |X posee las probabilidades de
acceder del estado Y al estado X.
Con esta información uno estarı́a interesado en obtener la matriz de probabilidad
de transición de la distribución marginal de X (Denotada por AX|X ). Para tal fin
0
el interés recae sobre la subsucesión X de la sucesión 3.1. Dado que 3.1 produce
0 0 0 0 0 0
X0 → Y1 → X1 → Y2 y ya que X0 → X1 forma una cadena de Markov e implica tener
0
que pasar por Y1 . Podemos escribir esto en términos de la probabilidad condicional
gracias a la ecuacion 2 para obtener
0 0
X 0 0
0 0

P X1 = x1 |X0 = x0 = P X1 = x1 |Y1 = y × P Y1 = y|X0 = x0 ,
y
con esto AX|X viene dado por

0.2 0.4
0.2 0.1
38 25

AX|X = AY |X AX|Y = 0.6 0.6 0.3 0.3 = 63 63 . (3.12)
0.1 0.3 0.4 0.3 25 17
0.4 0.4 0.7 0.7 42 42
0 0
En resumen hemos construidos la matriz de transición para ir de X0 a X1 . Esto
es fundamental para el presente trabajo pues se puso en evidencia la naturaleza
0 0
Markoviana de la iteración del Muestreador de Gibbs al considerar que ir de X0 a X1
forma una Cadena de Markov y poder escribir AX|X como un producto de matrices de
transición. En general la idea fundamental de los algoritmos MCMC, es el diseño de
una cadena de Markov cuya distribución estacionaria es precisamente la distribución
que se quiere muestrear. Para el caso del muestreador de Gibbs presentamos la
generalización del ejemplo anterior y uno empleando la distribución trinomial.
En este orden de ideas consideremos consideremos la tabla 3.2, con muestreo multi-
nomial dado que X y Y son variables aleatorias con distribución conjunta Bernoulli.
H
HH X
0 1
Y H
HH
0 p1 p2
1 p3 p4
Tabla 3.2. Distribución de probabilidades conjunta de la variable aleatoria discreta (X, Y ).
Teniendo en cuenta que pi ≥ 0, para i ∈ {1, 2, 3, 4} y es tal que p1 + p2 + p3 + p4 = 1.

Luego expresando la tabla 3.2 en términos de la función de probabilidad conjunta,
tenemos que
fX,Y (0, 0) fX,Y (1, 0) p1 p 2
=
fX,Y (0, 1) fX,Y (1, 1) p3 p 4
se sigue también que la distribución marginal de X es una distribución Bernoulli con
probabilidad de éxito p2 + p4 .
fX = [fX (0), fX (1)] = [p1 + p3 p2 + p4 ] (3.13)

de manera análoga a como se hizo en el ejemplo anterior se puede construir dos

matrices con todas las distribuciones condicionales de X|Y y Y |X a saber:
p1 p2
AX|Y = p1p+p
3
2 p1 +p2
p4
p3 +p4 p3 +p4
y p1 p3
p1 +p3 p1 +p3
AY |X = p2 p4 .
p2 +p4 p2 +p4
Si consideramos el espacio de estados {0, 1}, es natural pensar en las matrices

AX|Y y AY |X como matrices de transición. Donde por ejemplo AX|Y contiene las
probabilidades de ir del estado X al estado Y y AY |X posee las probabilidades de
acceder del estado Y al estado X. Además sı́ aplicamos el esquema iterativo (3.1) se
obtendrı́a una sucesión de ceros y unos.
Nuevamente nos preguntamos ¿cómo obtener la distribución marginal de X?, para
hacerlo centraremos nuestra atención en dos hechos, primero en la subsucesión de
0
X de la sucesión de Gibbs y segundo en que los valores en 3.1 se obtuvieron fijando
un valor inicial y evaluando en las condicionales respectivas su inmediato anterior.
0 0 0
Ası́ que si consideramos la sucesión de Gibbs como X0 → Y1 → X1 esta nos permite
0 0 0 0 0
acceder de X0 a X1 pasando por Y1 . Además ya que X0 → X1 forma una cadena
de Markov podemos aplicar la ecuación 2 para expresar esto en términos de la
probabilidad de transición, es decir,
0 0
X 0 0
0 0

P X1 = x1 |X0 = x0 = P X1 = x1 |Y1 = y × P Y1 = y|X0 = x0 (3.14)
y
de donde resulta que la matriz de probabilidad de transición AX|X de la sucesión

0
X , viene dada por
 2 
1 p1 p23 1 p1 p2 p3 p4
p +p p +p
+ p3 +p4
+
AX|X = AY |X AX|Y =  1 3 p1 p 2 p1 +p3 p1 +p 2 p3 +p4
 . (3.15)
1 2 1 p4 p3 1 p22 p24
p2 +p4 p1 +p2
+ p3 +p4 p2 +p4 p1 +p2
+ p3 +p4
0 0 k
De manera que la matriz de transición para P Xk = xk |X0 = x0 es AX|X . Esto
implica que podemos calcular fácilmente la distribución de probabilidad de cualquier
0 0
Xk en la sucesión 3.1. En efecto, si denotamos la distribución marginal de Xk por
fk = [fk (0) , fk (1)]
entonces para cualquier k, tendremos que,

fk = f0 AkX|X = f0 Ak−1
X|X AX|X = fk−1 AX|X . (3.16)
El lector en este punto puede llegarse a preguntar ¿que relación tiene esto con el
muestreador de Gibbs? la respuesta la daremos de inmediato. En efecto, ya que como
se ha indicado desde el inició de este capı́tulo el muestreador de Gibbs es uno de
los algoritmos MCMC, los cuales tienen como idea fundamental el diseño de una
cadena de Markov cuya distribución estacionaria es precisamente la distribución que
se quiere muestrear. Esto y el hecho de que cuando la matriz AX|X tiene todas sus
entradas positivas, la ecuación 3.16 garantiza que cuando k → ∞, fk converge a la
distribución f que es en un punto estacionario de 3.16. Sin importar la probabilidad
inicial f0 . Mas aún el hecho de que se satisface,
f AX|X = f. (3.17)
Responde al lector la pregunta planteada, es decir, la relación de este desarrollo y el
muestreador de Gibbs es que hemos visto la naturaleza markoviana del muestreador
de Gibbs. Ası́ mismo de la ecuación 3.17 para el caso descrito, podemos entender
porque el muestreador de Gibbs es uno de los algoritmos MCMC.
Veamos ahora que para fX en 3.13 se verifica 3.17,
 2 
1 p1 p23 1 p1 p2 p3 p4
p +p p +p
+ p3 +p4 p1 +p3
+
fX AX|X = [p1 + p3 p2 + p4 ]  1 3 p1 p 2 p1 +p 2 p3 +p4

1 2 1 p4 p3 1 p22 p24
p2 +p4 p1 +p2
+ p3 +p4 p2 +p4 p1 +p2
+ p3 +p4
= [p1 + p3 p2 + p4 ]
= fX
Además, si la sucesión 3.1 converge, la f que satisface 3.17 es la distribución marginal

de X. Para indicar este hecho se ha implementado el siguiente ejemplo con el software
estadı́stico R.
Ejemplo 1.
Consideremos la distribución dada por la tabla 3.3
HH
X
HH 0 1 2 3 suma
Y HH
0 0.008 0.036 0.054 0.027 0.125
1 0.06 0.18 0.135 0 0.375
2 0.150 0.225 0 0 0.375
3 0.125 0 0 0 0.125
suma 0.343 0.441 0.189 0.027 1
Tabla 3.3. Distribución de probabilidades conjunta para el ejemplo 1.
Sea P la matriz obtenida de expresar la tabla 3.3 en términos de la función de

probabilidad conjunta,es decir,
   
fX,Y (0, 0) fX,Y (1, 0) fX,Y (2, 0) fX,Y (3, 0) 0.008 0.036 0.054 0.027
fX,Y (0, 1) fX,Y (1, 1) fX,Y (2, 1) fX,Y (3, 1)  0.06 0.18 0.135 0 
 =  = P.
fX,Y (0, 2) fX,Y (1, 2) fX,Y (2, 2) fX,Y (3, 2) 0.150 0.225 0 0 
fX,Y (0, 3) fX,Y (1, 3) fX,Y (2, 3) fX,Y (3, 3) 0.125 0 0 0
Los renglones en la matriz AY |X se obtienen dividiendo los respectivos renglones

de P T por la suma de la correspondiente columna en P . Por ejemplo calculemos la
entrada (2, 1) de AY |X esta es 0.036/0.441 = 0.08163265. Ası́ obtenemos
 
0.02332362 0.1749271 0.4373178 0.3644315
0.08163265 0.4081633 0.5102041 0 
AY |X = 0.28571429 0.7142857

0 0 
1 0 0 0
para obtener AX|Y dividimos cada renglon de P por la suma de todos los elementos
en este renglón, por ejemplo calculemos la entrada (1, 1) en la matriz AX|Y , esto es,
0.008/0.125 = 0.064.
 
0.064 0.288 0.432 0.216
0.160 0.480 0.360 0 
AX|Y =  0.40

0.6 0 0 
1 0 0 0
0
de manera que la matriz de transición de la subsucesión X en la sucesión de Gibbs
viene dada por
 
0.5688397 0.3530729 0.07304956 0.005037901
0.2746122 0.5255510 0.18220408 0.017632653
AX|X = AY |X AX|Y = 
0.1325714 0.4251429 0.38057143 0.061714286 (3.18)

0.064 0.2880 0.4320 0.2160
0
la cadena de Markov de las X es irreducible pues todos sus estados se comunican,
es aperiodica ya que todos sus estados tienen periodo 1, como lo muestra la matriz
3.18. Además también posee distribución estacionaria f = [0.343, 0.441, 0.189, 0.027]
(que se corresponde con la marginal de X en la tabla 3.3). Ası́ que por el teorema
20, se tiene que para cualesquiera estados i y j, lı́mn→∞ pij (n) = fj , en efecto, para
n = 21 se tiene que:
 
0.343 0.441 0.189 0.027
0.343 0.441 0.189 0.027
A21
X|X = 
 .
0.343 0.441 0.189 0.027
0.343 0.441 0.189 0.027
Como un primer resultado se puede afirmar que para valores suficientemente grandes
0
de k la distribución marginal de Xk tiende a fX . Ası́ si detenemos la iteración (3.1)
0
a un valor suficientemente grande de k, podemos suponer que la distribución de Xk
es aproximadamente fX .
Casella y George (1992) argumentan que es posible considerar X, Y alguna o ambas
continuas, si bien no presentan qué condiciones se deben imponer para extender los
argumentos de dimensión finita, si afirman que la teorı́a continúa sirviendo y que
por tanto el muestreador de Gibbs determina una muy buena aproximación a la
distribución marginal de X. Es en este contexto que presentaremos un paralelo entre
la versión discreta que se ha desarrollado y el caso continuo. Para ello es pertinente
presentar pimero una versión para este último caso de la ecuación 3.14, es decir, una
0 0
formula que represente la densidad condicional de X1 dado X0 , esta es,
Z
f 0 0
X1 |X0 (x1 | x0 ) = fX 0 |Y 0 (x1 | y) fY 0 |X 0 (y | x0 ) dy.
1 1 1 0
Una vez hecha esta precisión, nuestro interés recae sobre tratar de describir una
k
matriz análoga a Ax|x , es decir, la matriz de transición de k−pasos. Teniendo
presente que X o Y pueden ser continuas, entonces el elemento buscado es “una
matriz de transición infinita ”. Cuyas entradas satisfacen la siguiente ecuación:
Z
fX 0 |X 0 (x | x0 ) = fX 0 |X 0 (x | t) fX 0 (t | x0 ) dt, (3.19)
k 0 k k−1 k−1 |X0
Por ultimo vale la pena mencionar que la ecuación 3.19 representa la versión continua
3.16, ası́ que la densidad fX 0 |X 0 representa una transición de un solo paso, y las
k k−1
otras dos densidades hacen referencia a fk y fk−1 . También en [5] se describe que
cuando k → ∞ se deduce que el punto estacionario de (3.19) es la densidad marginal
de X, la densidad para la cual fX 0 |X 0 converge.
k k−1
El siguientes aspecto a tratar es revelar el argumento por el que funciona el mues-

treador de Gibbs además establecemos la forma limite de la ecuación 3.19 todo esto
lo presentamos aquı́ en versión bivariada. Para comenzar sabemos que conocidas
las distribuciones condicionales podemos determinar la distribución conjunta si esta
existe, Lo que parece poner al muestreador de Gibbs como una aplicación de este
hecho. Pero es conveniente señalar que no siempre conocidas las distribuciones condi-
cionales estas determinen una distribución marginal adecuada. Para tal situación el
lector interesado puede consultar [5] donde se presenta un ejemplo de tal situación.
Hay que recordarle al lector que el muestreador de Gibbs es capaz de simular la
distribución marginal mediante las distribuciones condicionales correspondientes.
Esta idea es clave para el siguiente desarrollo. Para comenzar sean X y Y dos
variables aleatorias para las cuales conocemos las densidades condicionales fX|Y (x | y),
fY |X (y | x) y desconocemos fXY (X, Y ). Encontraremos ahora la densidad marginal
de X, fX (x), para hacerlo por definición sabemos que

Z
fX (x) = fXY (x, y) dy,
donde fXY (x, y) es la densidad conjunta desconocida, pero sabemos que fXY (x, y) =
fX|Y (x | y) fY (y), se deduce que
Z
fX (x) = fX|Y (x|y) fY (y) dy, (3.20)
R
ahora reemplazando fY (y) = fY |X (y|t) fX (t) dt, en 3.20 tenemos que
Z Z
fX (x) = fX|Y (x|y) fY |X (y|t) fX (t) dtdy
 
Z Z 
=  fX|Y (x|y) fY |X (y|t) dy  fX (t) dt (3.21)
 
 
| {z }
h(x,t)
Z
= h (x, t) fX (t) dt.
Es decir, escribimos la distribución marginal de X en términos de las condicionales.

Razón por la cual el lector debe apreciar el porque del exito del muestreador de
Gibbbs. Además la ecuación (3.21) define una ecuación integral de punto fijo para la
cual fX (x) es solución. La unicidad de la solución se explica en [?].
Además de que la ecuación 3.21 coloca a la distribución marginal de X en términos
de las condicionales. Un hecho que no se mostrara aquı́ pues su desarrollo excede el
objetivo del presente trabajo, es que cuando k → ∞ en (3.19), se tiene que
fX 0 |X 0 (x|x0 ) → fX (x) y fX 0 |X 0 (x|t) → h (x, t) (3.22)

k 0 k k−1
muestran que 3.21 es la forma limite de 3.19.
3.4. Evaluación de P n
La ecuación 3.16 en la sección anterior motiva trabajar con la k-ésima potencia de la

matriz de transición. Por otro lado, evaluar una matriz cuadrada digamos n × n con
un exponente k ≥ 5, es una trabajo bastante tedioso en especial si n ≥ 4. Ası́ pues,
el objetivo de esta sección es presentar la matriz de transición para una cadena de
Markov con espacio de estados finito, cuando es evaluada para alguna potencia en
los enteros mayores que cero. Se ha decidido seguir el texto [2], para el presente
desarrollo.
Sea P la matriz de transición para una cadena de Markov con k estados. Supongamos
que los autovalores de P son distintos entre sı́. Además la ecuación
0 0
P x = λx y y P = λy
tiene solución x y y distinta de cero si y sólo si λ es una autovalor de P . Se sigue

que para cada λj , podemos encontrar los correspondientes autovectores a derecha y
0 0
a izquierda xj y yj .
Ahora si definimos H = [x1 , x2 , · · · , xk ] y K = [y1 , y2 , · · · , yk ], es posible mediante
algunos cálculos con estos elementos, llegar a caracterizar la matriz P n como
k
X
Pn = λnj Aj ,
j=1
0 0
dado que yj xj = 1 y Aj definido como xj yj . Sin embargo en [2], también se considera
0
el caso donde no se ajusta la condición yj xj = 1 sino yj xj = cj , con el fin de obtener
k
X
n
P = c−1 n
j λj Aj .
j=1
También es posible presentar la matriz P n cuando P tiene multiples autovalores,

dada por, rj −1 n
n
X 1 d λ adj (λI − P )
P =
j
(rj − 1)! dλrj −1 ψj (λ) λ=λj
donde ψj (λ) = i6=j (λ − λi )ri y rj es la multiplicidad de λj (Frazer, R.A., Duncan,

Q
W.J., AND Collar, A.R., 1946).
Conclusiones
• El muestreador de Gibbs es un algoritmo versátil para simular las distribuciones

marginales de una distribución conjunta. Además dado que la mecánica que
lo describe es sencilla su aplicación en software como R no presenta mayor
dificultad y permite apreciar el buen comportamiento que tiene la distribución
simulada mediante el muestreador de Gibbs a la distribución teórica. Por otro
lado la simulación desarrollada también permite apreciar hechos como lo rápido
que converge la matriz de transición a AX|X a la distribución marginal de la
variable aleatoria X.
• Las cadenas de Markov hacen parte del soporte teórico del muestreador de
Gibbs. Pues permiten construir la matriz de transición en k pasos para la
variable aleatoria que se quiere simular, considerando los valores de la variable
en la sucesión de Gibbs como una cadena de Markov. La importancia de
construir dicha matriz es que su distribución estacionaria coincide con la que
se quiere simular, como se mencionó y se ilustró por medio de ejemplos en el
trabajo.
• El trabajo desarrollado también permite apreciar que el éxito del muestreador
de Gibbs se basa en que la función marginal de X se puede escribir en términos
de las distribuciones condicionales X | Y y Y | X.
38
APÉNDICE
Códigos en R
Empleando el Muestreador de Gibbs para la distri-

bución Beta-binomial
x<-c()
y<-c(0.4) ##Valor inicial para la sucesión de Gibbs##
##Algoritmo del Muestreador de Gibbs##
for(i in 1:100000){
x[i] <- rbinom(1,10,y[i])
y[i+1] <-rbeta(1,x[i]+2,15-x[i])
}
###Sucesión de Gibbs####
datos<-table(x)/100000
####Distribución Teórica Beta-binomial

con parámetros n=10, alfa=2 y beta=5#######
dbb <- function(x){
(choose(10,x)*30/factorial(16))*gamma(x+2)*gamma(15-x)
}
#### Evaluando la distribución Beta-binomial
a <- 0:10
t <- dbb(a)
39
CÓDIGOS EN R 40
Prueba ji-cuadrado para la distribución Beta-

binomial
x<-c()
y<-c(0.4) ##Valor inicial para la sucesión de Gibbs##
for(i in 1:100000){
x[i] <- rbinom(1,10,y[i])
y[i+1] <-rbeta(1,x[i]+2,15-x[i])
}
####Distribución Teórica Beta-binomial

dbb <- function(x){
(choose(10,x)*30/factorial(16))*gamma(x+2)*gamma(15-x)
}
#### Evaluando la distribución Beta-binomial
a <- 0:10
t <- dbb(a)
###Cuadrado de bondad de ajuste####

frecbb <- as.vector(table(x))/100000
estad_chi <- sum((frecbb-t)^2/t)
pvalor <- p_valor<- pchisq(estad_chi,10,lower.tail=FALSE)
Empleando el Muestreador de Gibbs para la distri-

bución Binomial
x<-c()
y<-c(2) ##Valor inicial para la sucesión de Gibbs##
for(i in 1:200000){
x[i] <- rbinom(1,10-y[i],0.3/(1-0.2))
y[i+1] <-rbinom(1,10-x[i],0.2/(1-0.3))
}
###Sucesión de Gibbs####
datos2<-table(x)/200000
#### Evaluando la distribución Binomial #######
db <- function(x){
choose(10,x)*((0.3)^(x))*((0.7)^(10-x))
CÓDIGOS EN R 41
}
l<-db(0:10)
Prueba ji-cuadrado para la distribución Binomial
x<-c()
y<-c(2) ##Valor inicial para la sucesión de Gibbs##
for(i in 1:200000){
x[i] <- rbinom(1,10-y[i],0.3/(1-0.2))
y[i+1] <-rbinom(1,10-x[i],0.2/(1-0.3))
}
#### Evaluando la distribución Binomial #######
db <- function(x){
choose(10,x)*((0.3)^(x))*((0.7)^(10-x))
}
l<-db(0:10)
###Cuadrado de bondad de ajuste####
frecb <- as.vector(table(x))/200000
estad_chi <- sum((frecb-l)^2/l)
pvalor <- p_valor<- pchisq(estad_chi,10,lower.tail=FALSE)
Código Ejemplo 1
install.packages("Matrix",dependencies=TRUE)
install.packages("expm",dependencies=TRUE)
#################################
#Distribución conjunta trinomial#
#################################
trinomial<- function(x,y,n,p,q){
if(x+y<=n)
(factorial(n)/(factorial(x)*factorial(y)*factorial(n-x-y)))
*p^x*q^y*(1-p-q)^(n-x-y)
else 0 }
#################################
#Matriz probabilidades trinomial#
#################################
conjunta <- function(n,p,q){
T <- matrix(nrow=n+1,ncol=n+1)
for(i in 0:n){
for(j in 0:n){
T[i+1,j+1] <- trinomial(i,j,n,p,q)
CÓDIGOS EN R 42
}
}
print(T)}
########################################################
#Obtención de la matriz de probabilidades de transición#
########################################################
##########
#Conjunta#
##########
fxy <- conjunta(3,0.5,0.3)
###############
#Marginal de x#
###############
fx <- apply(fxy,2,sum)
###############
#Marginal de y#
###############
fy <- apply(fxy,1,sum)
###################################################
#Marginal de X a partir de la matriz de transición#
###################################################
Ayx <- t(fxy)/fx
Axy <- fxy/fy
Ax <- Ayx%*%Axy
PAx<- Ax %^% 21
###################################################
#Marginal de Y a partir de la matriz de transición#
###################################################
Ay <- Axy%*%Ayx
PAy <- Ay %^% 21
Bibliografı́a
[1] Alan, G AND Smith, A. Sampling-Based Approaches to Calculating Marginal

Densities. Journal of the American Statistical Association, 85(40), 398–409, 1990.
[2] Bailey, N. The elements of Stochastic Processes with applications to the natural
sciences. United States of America: John Willey & Sons, Inc., 1964.
[3] Blanco, L. Probabilidad. Bogotá, Colombia: Unibiblos, 2004.
[4] Blanco, L., Arunachalam W., AND Dharmaraja, D. Introduction to
Probability and Stochastic Processes with Applications. New Jersey, USA: John
Wiley & Sons, Inc., 2012.
[5] Casella, G AND Edward, G. Explaining the Gibbs Sampler. The American
Statistical Association, 46(3), 167–174, 1992.
[6] Frazer, R.A., Duncan, W.J., AND Collar, A.R. Elementary Matrices.
Cambridge University Press, 1946.
[7] Rincón, L. Curso intermedio de Probabilidad. Mexicó DF, Mexicó: Departa-
mento de Matemáticas Facultad de Ciencias UNAM, 2007.
[8] Rincón, L. Introducción a los procesos estocásticos. Mexicó DF, Mexicó: Depar-
tamento de Matemáticas Facultad de Ciencias UNAM, 2012.
[9] Ross, S. Introduction to probability Models. San Diego, United States of America:
John Willey & Sons, Inc., 2010.
[10] Prieto, V. VII Coloqio colombiano de matemáticas, Cadenas de Markov.
Bogotá, Colombia: Universidad Nacional de Colombia, 1977
43

Muestreador de Gibbs

Cargado por

Copyright:

Formatos disponibles

Muestreador de Gibbs

Cargado por

Información del documento

Derechos de autor

Formatos disponibles

Compartir este documento

Compartir o incrustar documentos

Opciones para compartir

¿Le pareció útil este documento?

¿Este contenido es inapropiado?

Copyright:

Formatos disponibles

Muestreador de Gibbs

Cargado por

Copyright:

Formatos disponibles

Descripción e implementación del Muestreador

de Gibbs en versión bivariada

Manuel Alejandro Moreno Arévalo

Universidad Distrital Francisco Jose de Caldas

Manuel Alejandro Moreno Arévalo

Disertación presentada para optar al tı́tulo de

Universidad Distrital Francisco Jose de Caldas

Resumen: Este trabajo de grado presenta un revisión bibliográfica de los

El trabajo se estructuró en tres capı́tulos:

Abstract: This degree work presents a literature review of the mathematical

The work was divided into three chapters:

Palabras clave: Distribución condicional, Cadenas de Markov, Sucesión de Gibbs.

Bogotá, D.C., Diciembre 7 de 2016

Mi mamá Rosalia Arévalo y mi abuelita Clementina Arévalo.

Índice de tablas III

3.2.2. Aplicando el Muestreador de Gibbs . . . . . . . . . . . . . . . . . . . 23

3.1. Distribución de probabilidades de (X, Y ) . . . . . . . . . . . . . . . . . . . . 30

3.1. Comparación de histogramas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

Algunos cálculos en la teorı́a de probabilidad requieren de una cierta habilidad

el muestreador de Gibbs y se contrastanlos resultados con la distribución teórica.

El objetivo de este capı́tulo es presentar algunos conceptos básicos de la teorı́a de la

1.1. Algunos conceptos de la teorı́a de la probabili-

Definición 1. Sea X1 , X2 , · · · Xn n-variables aleatorias reales definidas sobre el

X (w) := (X1 (w) , · · · , Xn (w))

es llamada un n−dimensional vector aleatorio.

Definición 2. Sea X un n−dimensional vector aleatorio. La medida de probabilidad

Nota 1. Bn denota la σ−álgebra de Borel sobre Rn .

La naturaleza discreta de las variables aleatorias de un vector aleatorio X permite

Definición 3. Sea X = (X1 , X2 , · · · , Xn ) un vector aleatorio n−dimensional. Si

Es importante añadir que si en la definición anterior las variables aleatorias Xi ,

Definición 4. Sea X = (X1 , X2 , · · · , Xn ) un vector aleatorio n−dimensional. La

para toda (x1 , x2 , · · · , xn ) ∈ Rn es llamada la función de distribución conjunta

Teorema 1. Sea X = (X1 , X2 , · · · , Xn ) un vector aleatorio n−dimensional con

FXj (x) = lı́m · · · lı́m lı́m · · · lı́m F (x1 , · · · xn ) .

La función de distribución FXj es llamada la funcion de distribución marginal

Teorema 2. Sea X = (X, Y ) un vector aleatorio 2−dimensional. La función de

Teorema 3. Sea X = (X1 , X2 , · · · , Xn ) un vector aleatorio n−dimensional. La fun-

2. F es continua a derecha en cada componente.

lı́m F a1 , · · · , ai−1 , x , ai+1 , · · · , an  = 0

Definición 5. Sean X1 , X2 , · · · , Xn n variables aleatorias de valor real definidas

La función f es llamada la función de densidad de probabilidad conjunta de las

Nota 2. Es costumbre notar fdp a la función de densidad de probabilidad.

Teorema 4. Sean X y Y variables aleatorias continuas con funcion de distribución

para todos los puntos (x, y) donde f (x, y) es continua.

Teorema 5. Si X1 , X2 , · · · , Xn son n variables aleatorias de valor real, teniendo f

es la función de densidad de la variable aleatoria Xj para j = 1, 2, · · · , n.

Definición 6. Sean X y Y dos variables aleatorias discretas. La función de probabi-

Para todo y para el cual P (Y = y) > 0.

Para todo y para el cual P (Y = y) > 0.

Definición 8. Sean X y Y variables aleatorias continuas con función de densidad

para todo y con fY (y) > 0.

Definición 9. La función de distribución condicional de X dado Y = y se define

Para todo y con fY (y) > 0.

1.2. La prueba Ji-cuadrado

Donde n es el tamaño de la muestra y npi son lo valores medios teóricos. Un hecho

de los objetivos del trabajo es que la distribución del estadı́stico X 2 es un distribución

Debido a la naturaleza Markoviana de la iteración del muestreador de Gibbs se

Definición 10. Una cadena de Markov es un proceso estocástico a tiempo discreto

La propiedad de Markov establece que la probabilidad condicional de un evento

Definición 11. Sea