TFG: Mé odos Mon e
Ca lo basados en
cadenas de Ma ko
José Jiménez
JIMENEZLUNA.COM
Licensed unde he C ea i e Commons A ibu ion-NonComme cial 4.0 Unpo ed License ( he
“License”). You may no use his ile excep in compliance wi h he License. You may ob ain a
copy o he License a
h p://c ea i ecommons.o g/licenses/by-nc/4.0
. Unless equi-
ed by applicable law o ag eed o in w i ing, so wa e dis ibu ed unde he License is dis ibu ed
on an “AS IS”BASIS,WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, ei he exp ess o
implied. See he License o he speci ic language go e ning pe missions and limi a ions unde
he License.
P ime a edición, Ene o 2015
Abs ac :
Ma ko Chain Mon e Ca lo (o sho ly MCMC) is a powe ul me hod o sampling
om high dimensional p obabili y dis ibu ions. He e we p esen a swi in oduc ion o he
heo y behind hese me hods as well as se e al applica ions o he Bayesian MCMC amewo k.
These include, among o he s Bayesian Mix u e Models, Bayesian Image Analysis o Tex Mining.
Implemen a ions o solu ions ega ding hese p oblems can be ound in se e al p og amming
languages.
Índice gene al
1Mo i ación del p oblema .................................... 7
1.1 P oblemá ica en in e encia bayesiana 7
1.2 Cálculo de espe anzas 8
2Concep os p e ios .......................................... 9
2.1 Mues eo po echazo 9
2.2 Cocien e de uni o mes 10
2.3 In eg ación Mon e Ca lo 10
2.4 Mues eo po impo ancia 11
2.5 Cadenas de Ma ko 12
2.6 P opiedades de las Cadenas de Ma ko 13
3Ma ko Chain Mon e Ca lo ................................. 15
3.1 Algo i mo de Me opolis-Has ings 15
3.1.1 Dis ibucionesp oposición....................................... 16
3.2 Algo i mos de mues eo 17
3.2.1 Algo i modeMe opolis ........................................ 17
3.2.2 Paseoalea o ioMe opolis ...................................... 17
3.2.3 Mues eoconindependencia ................................... 17
3.2.4 Ac ualizandoenbloques ....................................... 18
3.3 Mues eo de Gibbs 19
3.4 Mues eo po echazo adap a i o (ARS) 20
3.5 O as conside aciones 21
3.5.1 Valo esiniciales ............................................... 21
3.5.2 Calen amien o ............................................... 21
3.5.3 Análisisdelasalida ............................................ 22
3.5.4 Tasadecon e gencia.......................................... 22
3.5.5 Es imacióndela a ianza ....................................... 22
4Modelos g a icos y DAGs ................................... 25
4.1 Dig a os acíclicos 25
4.2 G a o de independencia condicional 26
4.3 Ejemplos de aplicación 27
5Mix u as gaussianas ........................................ 33
5.1 Modelos de mix u as ini os 33
5.2 Usando MCMC 35
5.3 Di icul ad en easignación 39
5.4 De e minando el núme o de subpoblaciones 42
6Análisis de imagen ......................................... 45
6.1 Recons ucción de imágenes 45
6.1.1 Camposalea o iosdeMa ko ................................... 46
6.1.2 ModelodeIsing............................................... 46
6.1.3 ModelodePo s .............................................. 51
6.1.4 Sob e la cons an e de in eg ación . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
6.1.5 Segmen acióndeimágenes..................................... 56
7Impu ación múl iple ........................................ 63
7.1 Impu ación usando ecuaciones encadenadas 63
7.1.1 Modelos de impu ación uni a ian e . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
7.2 P obando el algo i mo 66
8Mine ía de ex o ............................................ 69
8.1 Descub imien o de emá icas 69
8.2 Asignación la en e de Di ichle 69
8.2.1 Fo malización de la gene ación . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
8.3 In e encia y es imación usando mues eo de Gibbs colapsado 70
8.3.1 Fo malización del ap endizaje . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
8.4 Ejemplos de aplicación 72
8.4.1 Análisisdepe ilenTwi e ....................................... 73
8.4.2 Temá ica según au o es clásicos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
Bibliog a ía ................................................. 77
P oblemá ica en in e encia bayesiana
Cálculo de espe anzas
1. Mo i ación del p oblema
Los mé odos MCMC (Ma ko Chain Mon e Ca lo) su gen de la necesidad de simula el
compo amien o de a iables alea o ias y de es ima pa áme os de las unciones de densi-
dad/p obabilidad de las mismas. El g an impulso a es as écnicas se las da mayo men e (pe o
no únicamen e) el en oque es adís ico bayesiano, donde la in e encia se ealiza sob e lo que se
denomina una unción a pos e io i, que denomina emos π(θ|x).
Las siglas MCMC ienen ma cadas po las cadenas de Ma ko y po la in eg ación Mon e
Ca lo. La in e encia bayesiana, en mul i ud de ocasiones necesi a in eg a sob e dis ibuciones
de dimensión muy ele ada (en muchas ocasiones, con cien os de pa áme os). Exis en mé odos
numé icos ap oximados que p oducen buenas soluciones, pe o que no escalan bien con la di-
mensión, siendo en muchos casos, compu acionalmen e in a ables. De mane a poco o mal, la
aplicación de écnicas MCMC cons a de dos pasos:
1.
Gene a una mues a
X1,...,Xn
median e una cadena de Ma ko cuya dis ibución es acio-
na ia sea la buscada.
2.
Toma medias mues ales (in eg ación Mon e Ca lo) y ealiza in e encias sob e la mues a
an e io men e ci ada
Algunos é minos no quedan explicados a es as al u as. No obs an e, es a es uc u a gene al
que queda á acla ada en los siguien es capí ulos pe mi e esol e muchos p oblemas.
1.1 P oblemá ica en in e encia bayesiana
Como hemos mencionado, en in e encia bayesiana, se ealiza in e encia sob e una dis i-
bución a pos e io i
π(θ|x)
. Bajo un en oque bayesiano, no exis e di e encia concep ual en e
pa áme os y alo es obse ables, es deci , son en su o alidad can idades alea o ias. Denomine-
mos po
x
a los da os obse ados, y
θ
al conjun o de pa áme os (nó ese que pueden se ambos
mul idimensionales). En es adís ica bayesiana, la unción a pos e io i puede exp esa se de la
siguien e mane a, aplicando el eo ema de Bayes:
π(θ|x) = L(x|θ)π(θ)
RL(x|θ)π(θ)(1.1)
Donde
L(x|θ)
no es más que la unción de e osimili ud de los da os asumiendo que siguen
8Capí ulo 1. Mo i ación del p oblema
una dis ibución pa ame izada y
π(θ)
es lo que se denomina una dis ibución a p io i. És a
úl ima exp esa en in e encia bayesiana oda in o mación p e ia o c eencia que se iene sob e el
compo amien o de los pa áme os del modelo. O lo que es lo mismo, en in e encia bayesiana
los pa áme os son ambién a iables alea o ias, cuya in o mación se incluye en un análisis a
pos e io i, omando ambién como e e encia la mues a xensayada.
En muchísimas ocasiones (po no deci casi odas), es a dis ibución a pos e io i se conoce
únicamen e sal o cons an e mul iplica i a (que ue za a que la dis ibución in eg e
1
sob e su
dominio Ω). Po es a azón, gene almen e se suele no a de la siguien e mane a:
π(θ|x)∝L(x|θ)π(θ)(1.2)
Es deci , se dice que la unción a pos e io i es p opo cional al nume ado de 1.1. Una ez que
enemos de e minada es a dis ibución, el en oque bayesiano ealiza in e encia sob e espe anzas
de unciones de la misma.
1.2 Cálculo de espe anzas
Como acabamos de menciona , la p oblemá ica aho a su ge de es ima
E[ (x)]
sob e la
dis ibución a pos e io i. En la mayo ía de ocasiones abaja emos sob e espacios pa amé icos
ele ados, donde las soluciones analí icas pasan a se un dolo de cabeza, y donde las numé icas
no escalan bien. (De hecho, se pueden ealiza simulaciones donde se comp ueba que en espacios
pa amé icos con más de 20 dimensiones es os mé odos suelen alla es epi ósamen e.)
En el es o del abajo seguimos una hoja de u a encaminada a ap ende cómo esol e
es os p oblemas mencionados de acue do con los pasos ci ados en el pun o an e io . An es de
ello, ealiza emos un muy b e e epaso a algunos mé odos de simulación de a iables alea o ias
básicos (que se án ú iles a lo la go del abajo), luego eco da emos algunas p opiedades básicas
de las cadenas de Ma ko que nos da án sopo e eó ico sob e lo que que emos hace , y pa a
inaliza los concep os p e ios habla emos muy b e emen e de in eg ación Mon e Ca lo.
Una ez inalizada la in oducción, explicamos los mé odos MCMC en gene al, dando a ios
algo i mos gene ales pa a esolución de dichos p oblemas (con ejemplos sencillos), algunos
undamen os eó icos pa a sus en a alidación de esul ados y algunos mé odos de simulación
más a anzados.
La úl ima pa e (y la más pesada) del abajo co esponde a a ias aplicaciones p ác icas
de es os mé odos. Algunos ejemplos son a amien o de imágenes o mix u as. En es a pa e del
abajo especialmen e se ha á un uso in ensi o de lenguajes de p og amación, como pueden se
R o Py hon 2.7. Las lib e ías u ilizadas en cada pa e se de allan pa a pos e io ep oducibilidad.
Todo el código del abajo queda disponible en el eposi o io
gi hub.com/hawk31/MCMC_ g
bajo licencia MIT.
Mues eo po echazo
Cocien e de uni o mes
In eg ación Mon e Ca lo
Mues eo po impo ancia
Cadenas de Ma ko
P opiedades de las Cadenas de Ma ko
2. Concep os p e ios
En es e capí ulo epasamos algunas ideas p e ias pa a la mejo comp ensión de las aplicacio-
nes pos e io es. Es as ideas básicas se i án pa a esboza el ma co de ac uación bayesiano an e
un p oblema de e minado.
2.1 Mues eo po echazo
El mues eo po echazo p opo ciona una mane a muy e icien e de simula alo es de una
a iable alea o ia
(x)
. Suponemos en es e sen ido que no enemos un algo i mo sencillo (como
puede se median e in e sión) pa a gene a alo es alea o ios de dicha dis ibución. La idea es u i-
liza una dis ibución ins umen al
g(x)
, que aco e absolu amen e a
(x)
, es deci
(x)<Mg(x)
,
con M>1 y de la que sepamos gene a alo es alea o ios de una mane a ácil.
El algo i mo unciona de la siguien e mane a:
Mues ea un alo xde g(x)yupe enecien e a una U(0,1).
Si
u< (x)
Mg(x)
, acep a
x
como alo alea o io de
(x)
. En caso con a io echaza el alo
y ol e al p ime paso.
Nó ese que es e mé odo supone que an o
como
g
son e aluables y que se conoce la
cons an e mul iplica i a
M
has a cie o pun o. Es e mé odo, aunque e icien e, iene algunas
des en ajas, como que la dis ibución
g
debe se pa ecida en cu osis y sime ía a
, de mane a
que en el segundo paso del algo i mo se acep e una p opo ción de eces acep able. En caso de
que es o no ocu a, desechamos muchas mues as. En es e sen ido, exis en mé odos un poco
más complejos (Mues eo po echazo adap a i o) que a an es e p oblema de una mane a más
a ac i a.
El mues eo po echazo es p ecu so di ec o del algo i mo de Me ópolis-Hash ings, que
además iene la en aja de escala mucho mejo con la dimensionalidad de la simulación.
Podemos ealiza una pequeña simulación de alo es alea o ios de p oceden es de una N(0,1)
u ilizando como dis ibución ins umen al una Cauchy(0,2), con cons an e
M=3
. (No a: si se
desea e la adecuación de es a dis ibución ins umen al bas a ía con ep esen a la)
16 Capí ulo 3. Ma ko Chain Mon e Ca lo
La dis ibución
q(.|.)
puede ene cualquie o ma y la cadena p opo cionada po el algo i mo
con e ge á a la dis ibución es aciona ia
π(x)
. Si bien es a condición es su icien e, la dis ibución
p oposición debe ene la misma dimensión que la es aciona ia, y se capaz de gene a alo es
que se acep en. En o o caso la cadena puede pasa pe iodos la gos de iempo en un mismo
es ado. Nó ese que en el momen o que
X
ya pe enezca a la dis ibución es aciona ia, odos los
alo es siguien es X +1,X +2,... pe enece án igualmen e a la misma dis ibución.
Eje cicio 3.1
Pa a ilus a cómo unciona el algo i mo de Me opolis-Has ings, u iliza emos
un ejemplo de jugue e. Supongamos que que emos mues ea alo es de una dis ibución
N(0,1)
,u ilizando como dis ibución p opues a
N(X ,0,5)
. Una posible implemen ación po-
d ía se la siguien e.
n = 1000
alea = nume ic(n)
alea[1]=0 #mu
o (i in 2:n)
{y = no m(1,alea[i-1],0.5)
u = uni (1)
alpha = min(1,(dno m(y)*dno m(alea[i-1],y,0.5))/(dno m(alea[i-1])*
dno m(y,alea[i-1],0.5)))
i (u<alpha) alea[i]=y
else alea[i]=alea[i-1]
}
En la p ác ica, se suele ealiza lo que se denomina ’ aza’ de la simulación, que no es
más que un g á ico de línea pa a comp oba la con e gencia y el mezclado de la cadena.
plo (alea, ype="l")
3.1.1 Dis ibuciones p oposición
Como hemos no ado an es, cualquie dis ibución p oposición
q(.|.)
es su icien e pa a que
la cadena con e ja a su dis ibución es aciona ia. No obs an e, la elección adecuada de es a
dis ibución y su elación con la es aciona ia ga an iza á una con e gencia más ápida hacia es a
3.2 Algo i mos de mues eo 17
úl ima. Es más, aún ga an izada con e gencia, la cadena pod ía ’mezcla ’ de mane a len a (es
deci , que se mue a len amen e po el sopo e de la dis ibución obje i o). Po ello, escoge una
dis ubución p opues a adecuada se con ie e en un p oblema. Exis en a ias o mas canónicas
que de allamos a con inuación.
3.2 Algo i mos de mues eo
T adicionalmen e, la li e a u a sob e MCMC ha hablado de mues eado es y algo i mos. No
obs an e, aunque aquí seguimos esas con enciones, con iene no ol ida que unos mues eado es
no excluyen a o os. Como e emos un poco más adelan e, es común combina los pa a cons ui
una cadena de Ma ko con mejo endimien o. Es po ello que iene más sen ido denomina los
como ’ac ualizado es MCMC’ más que mues eado es. Es os concep os queda án más asen ados
cuando expliquemos el mues eado de Gibbs.
3.2.1 Algo i mo de Me opolis
El o iginal p opues o po Me opolis en 1950 supone que u ilizamos dis ibuciones p oposi-
ción simé icas (es o es q(X|Y) = q(Y|X)). Con es a condicion 3.1, pasa a se :
α(X ,Y) = min1,π(y)
π(x)(3.2)
3.2.2 Paseo alea o io Me opolis
Si
q(x,y) = (y−x)
pa a una densidad pa icula
, el algo i mo se denomina de paseo
alea o io Me opolis. Se suelen usa densidades
simé icas, pa a acep a con p obabilidades
idén icas a las del apa ado an e io .
3.2.3 Mues eo con independencia
Supongamos a con inuación que
q(x,y) = (y)
, y po an o los posibles candida os a la
cadena se gene an independien emen e del es ado
X
de la cadena. En es e caso, la p obabilidad
de acep ación puede esc ibi se como:
α(x,y) = min1,w(y)
w(x)(3.3)
donde
w(x) = π(x)/ (x)
es lo que se denomina el peso de impo ancia. Es e mé odo gua da
g andes simili udes con las ideas mos adas en el mues eo po impo ancia. La di e encia esencial
en e los dos mé odos es que el mues eado po impo ancia gene a la masa de p obabilidad
sob e pun os con pesos al os, escogiéndolos ecuen emen e. En con as e, el mues eo de
independencia cons uye masa de p obabilidad en pun os con al os pesos, pe maneciendo en
ellos la gos pe iodos de iempo.
Eje cicio 3.2
Supongamos que que emos mues ea de una dis ibución Gamma de pa áme-
os
a,b
cualesquie a usando el mues eado de independencia. U iliza emos una dis ibución
N(a/b,a/b2)como p opues a. Una posible implemen ación pod ía se :
#! /us /bin/en RSc ip
gammaSample <- unc ion (n, a, b)
{mu <- a/b
sig <- sq (a/b^2)
18 Capí ulo 3. Ma ko Chain Mon e Ca lo
ec <- nume ic(n)
x <- mu
ec[1] <- x
o (i in 2:n) {
can <- no m(1, mu, sig)
ap ob <- min(1, (dgamma(can, a, b)/dgamma(x,
a, b))/(dno m(can, mu, sig)/dno m(x,
mu, sig)))
u <- uni (1)
i (u < ap ob) ec[i] <- can
else ec[i]<- ec[i-1]
}
ec
}
Nó ese que hemos u ilizado la media de la dis ibución obje i o como p ime alo de
la cadena pa a in en a consegui una con e gencia ápida. En la p ác ica es o no suele se
posible.
En la p ác ica, el mues eado de independencia puede unciona muy bien o muy mal. Pa a
que uncione bien,
q(.)
debe ía se una buena ap oximación a la dis ibución obje i o, aunque
gene almen e bas a con que enga colas más pesadas. En e ec o, si
q(.)
no cumple es a condición,
puede pasa mucho iempo a ascado en las colas de la dis ibución obje i o, p opo cionando mal
endimien o.
3.2.4 Ac ualizando en bloques
En los ejemplos que hemos es ado ealizando has a aho a, hemos u ilizado dis ibuciones con
un núme o bajo de pa áme os, pa a acili a la comp ensión. En la p ác ica desa o unadamen e
nos encon amos con si uaciones en las que podemos llega a ene cien os de pa áme os. Es po
ello que se han desa ollado p ocedimien os pa a ac ualiza cada uno de los componen es ’en
bloque’ o secuencialmen e. El mé odo o iginal p opues o po Me opolis se denomina single-
componen Me opolis-Has ings.
Supongamos que enemos que ac ualiza
X,1,....,X.h
. Es deci , enemos que ac ualiza h
elemen os, secuencialmen e. Deno emos po
X.−i
al mismo ec o sin el elemen o
i
-ésimo. Pa a
cada
i
de la i e ación
+1
, ac ualizamos cada componen e u ilizando Me opolis-Has ings. El
candida o
Y.i
se gene a median e dis ibución p opues a
qi(Y.i|X .i,X .−i)
. Es o quie e deci que
la dis ibución p opues a puede depende del es o de componen es del ec o y de sí misma
en la i e ación an e io . (Teniendo en cuen a que si se han ac ualizado algunos an es en la
misma i e ación
+1
pueden ac ualiza se los siguien es eniendo es o en cuen a). Po an o, pa a
cada componen e del ec o end emos una p obabilidad de acep ación al es ilo an e io , con
modi icaciones pe inen es:
α(X.−i,X.i,Y.i) = min1,π(Y.i|X.−i)qi(X.−i|Y.i,X.−i)
π(X.i|X.−i)qi(Y.−i|X.i,X.−i)(3.4)
En es a exp esión
π(.|X.−i)
se denomina dis ibución o almen e condicionada. Es o no es más
que la exp esión ijando odos los pa áme os menos el de in e és, a ándolos como cons an es.
El uso de dis ibuciones o almen e condicionadas juega un papel muy impo an e en los mé odos
MCMC y en especial in e encia bayesiana, donde se aplican en modelos condicionales de
3.3 Mues eo de Gibbs 19
dependencia (ya los e emos en aplicaciones), y que en la mayo ía de si uaciones o ecen una
simpli icación de 3.4.
3.3 Mues eo de Gibbs
El mues eado de Gibbs cob a especial impo ancia jus o después de es e úl imo apa ado.
El mismo nos p opo ciona una mane a gene al y sencilla de ac ualiza los componen es del
ec o de pa áme os mencionado u ilizando dis ibuciones o almen e condicionadas. U iliza
Me opolis-Has ings pa a cada componen e iene la des en aja e iden e de que enemos que
busca
h
densidades p opues a, una pa a cada uno de los componen es. El mues eado de Gibbs
p opone lo siguien e: bajo las mismas condiciones que an es, la dis ibución p opues a de cada
uno de los componen es se á:
qi(Y.i|X.i,X.−i) = π(Y.i|X.−i)(3.5)
Es deci , se p opone u iliza la dis ibución o almen e condicionada de un pa áme o al es o
de pa áme os ( alga la edundancia) como dis ibución p opues a. Es o, como hemos dicho
an es pasa po ija el es o de pa áme os y a a los como cons an es ( alo es iniciales o los de
la i e ación an e io ). Finalmen e, ac ualizamos la densidad esul an e (uni a ian e, es a ez) con
el mé odo que más con enga, ya sea Me opolis-Has ings u o o.
Es e mé odo iene g an aplicabilidad en odos los ejemplos p ác icos que e emos al inal del
abajo, ya que p opo ciona una mane a di ec a de ac ualiza los componen es.
Eje cicio 3.3
Pa a ejempli ica cómo unciona el mues eado de Gibbs, supongamos que
que emos gene a alo es alea o ios de la siguien e dis ibución con 3 a iables.
(y1,y2,y3)∝exp(−(y1+y2+y3+θ12y1y2+θ13y1y3+θ23y2y3)(3.6)
Pa a aplica el mues eado de Gibbs es necesa io de e mina las dis ibuciones o almen e
condicionadas de y1,y2,y3. De es a mane a:
(y1|y2,y3)∝exp(−y1(1+θ12y2+θ13y3)) (3.7)
Nó ese que podemos igno a aquellos é minos en los que no apa ece
y1
po que son cons an es
mul iplica i as. Con un poco de ojo, y eniendo en cuen a que aquí
y2
e
y3
son cons an es
nos damos cuen a que
y1|y2,y3∼Exp(1+θ12y2+θ13y3)
. Po sime ía y de mane a idén ica,
y2|y1,y3∼Exp(1+θ12y1+θ23y3)
y
y3|y1,y2∼Exp(1+θ12y1+θ23y2)
. El ejemplo no es a ía
comple o sin una simulación pe inen e, pa a la cual supond emos unos θi j cualesquie a.
ni e = 10^4
he a = c(2,3,1/2)
sim = ma ix(NA,n ow=ni e ,ncol=3)
sim[1,] = c(0.4,1.2,0.3) ## Valo es iniciales (cualesquie a)
o (i in 2:ni e ){
sim[i,1] = exp(1,1+ he a[1]*sim[i-1,2]+ he a[2]*sim[i-1,3])
sim[i,2] = exp(1,1+ he a[1]*sim[i,1]+ he a[3]*sim[i-1,3])
## Nó ese que aquí ya hemos ac ualizado el p ime
## componen e y lo usamos pa a acele a con e gencia.
sim[i,3] = exp(1,1+ he a[2]*sim[i,1]+ he a[3]*sim[i,2]) ## Idem
20 Capí ulo 3. Ma ko Chain Mon e Ca lo
}
En es e codigo hacemos uso de
exp
, la unción ya implemen ada en R pa a gene a
alo es según una exponencial de e minada. No obs an e, no enemos po qué ene un
algo i mo e icien e en alguna de es as dis ibuciones condicionadas, po lo que pod iamos
usa Me opolis-Has ings o alguna de las écnicas ya desc i as.
3.4 Mues eo po echazo adap a i o (ARS)
Ap o echa emos que ya hemos explicado el mues eado de Gibbs pa a explica una écnica
bas an e usada en conjun o. Has a aho a, cuando p e endíamos usa el mues eado de Gibbs
suponíamos que eníamos que simula dis ibuciones o almen e condicionadas al es o, y que
pa a ello pod iamos ecu i a Me opolis-Has ings, al mues eo po echazo o a cualquie o a
écnica que p opo cione buenos esul ados.
Si bien las an e io es écnicas uncionan bien, es necesa io de e mina en odas y cada una de
ellas una densidad p opues a de la que gene a los alo es que luego se án acep ados o echazados
con una de e minada p obabilidad. El mé odo que a con inuación se p opone es una écnica del
ipo ’caja neg a’ en el sen ido de que no necesi a es a de e minación ealizando una se ie de
hipó esis sob e la densidad de la que se p e ende mues ea .
Empezamos de iniendo algunos concep os. Supongamos que que emos mues ea de una
densidad
g(x)
, la cual podemos conoce has a una cons an e mul iplica i a (cons an e de in e-
g ación), con inua y di e enciable en odo su dominio. Supongamos además que h(x) = ln g(x)
es cónca a en odo D. Sea
Tk= (x1,...,xk)
,
k
pun os de abcisa en los que han sido e aluados
an o
h(x)
como
h0(x)
. A con inuación de inimos un con o no de echazo en
Tk
como
exp uk(x)
,
donde
uk(x)
es una unción lineal de inida a ozos de inida po las angen es de
h(x)
en
Tk
. Pa a
j=1,...,k−1, las angen es de xjyxj+1in e secan en:
zj=h(xj+1)−h(xj)−xj+1h0(xj+1)+xjh0(xj)
h0(xj)−h0(xj+1)(3.8)
Es a unción es una ap oximación ia con o no lineal po encima de
g(x)
.Po an o, pa a
x∈[zj−1,zj], de inimos:
uk(x) = h(xj)+(x−xj)h0(xj)(3.9)
donde
z0
es el ex emo in e io de
D
(o
−∞
si no es á aco ado in e io men e) y
zk
es el ex emo
supe io de
D
(o
∞
si no es á aco ado supe io men e.). Con es os concep os, de inimos la
densidad:
sk(x) = exp uk(x)
RDexp uk(x0)dx0(3.10)
Finalmen e, de inimos la unción es icción in e io en
Tk
como
explk(x)
, donde
lk(x)
es
una unción lineal de inida a ozos po debajo de
g(x)
, o mada en e abcisas adyacen es en
Tk
.
Pa a x∈[xj,xj+1]:
lk(x) = (xj+1−x)h(xj)+(x−xj)h(xj+1)
xj+1−xj
(3.11)
3.5 O as conside aciones 21
La conca idad de
h(x)
asegu a que
lk(X)≤h(x)≤uk(x)
en odo el dominio. T as odas es as
de iniciones, podemos pone en ma cha el algo i mo.
1.
Inicializa la abcisa en
Tk
. Si
D
no es á aco ado po la izquie da, escogemos
x1
al que
h0(x1)>0
. Si
D
es á no aco ado po la de echa, en onces escogemos
xk
al que
h0(xk)<0
.
2. Calculamos las unciones uk(x),sk(x)ylk(x)pa a es os kpun os.
3. Usamos sk(x)pa a mues ea un alo x∗ywde una uni o me unidad.
4. Si w≤exp(lk(x∗)−uk(x∗)) acep a x∗.
5. En o o caso, si w≤exp(h(x∗)−uk(x∗)), acep a x∗. En cualquie o o caso echaza lo.
6.
Si
h(x∗)
ha enido que se e aluado en es e úl imo es , incluímos
x∗
en
Tk
pa a o ma
Tk+1y ol emos al segundo paso has a que engamos npun os mues eados.
Es e algo i mo es el que se implemen a en la mayo ía de so wa e especializado dedicado al
mues eado de Gibbs, como pueden se
BUGS
(o sus múl iples a ian es) o
JAGS
, po lo que
iene g an impo ancia. En capí ulos pos e io es ha emos uso de es as he amien as pa a algunas
a eas. En
R
exis e una implemen ación en el paque e
a s
.
Eje cicio 3.4
Supongamos que que emos mues ea alo es de una Be a(2,3) median e
Adap i e Rejec ion Sampling. U ilizamos el paque e mencionado.
lib a y(a s)
#Mues ea de una Be a(2,3)
n=20
2<- unc ion(x,a,b){(a-1)*log(x)+(b-1)*log(1-x)}
2p ima<- unc ion(x,a,b){(a-1)/x-(b-1)/(1-x)}
mysample2<-a s(20, 2, 2p ima,x=c(0.3,0.6),m=2,
lb=TRUE,xlb=0,ub=TRUE,xub=1,a=2,b=3)
mysample2
3.5 O as conside aciones
En es e apa ado conside a emos algunos aspec os a ene en cuen a cuando se es án eali-
zando simulaciones MCMC en gene al, desde cómo de e mina unos alo es iniciales adecuados
has a eglas pa a de e mina un núme o de i e aciones adecuado.
3.5.1 Valo es iniciales
Si bien como hemos dicho los alo es iniciales
X0
no in luyen en la dis ibución es aciona ia
de la cadena, si és os no se escogen adecuadamen e algunas cadenas pueden mezcla len amen e,
eniendo que desca a muchas i e aciones como calen amien o. Una mane a de de e mina un
buen alo inicial es co e a ias cadenas simul áneamen e con dis in os alo es iniciales (si
se iene la su icien e po encia compu acional en pa alelo). Si se dispone de in o mación p e ia
sob e la dis ibución es aciona ia, un buen pun o de pa ida pod ía se la mediana.
3.5.2 Calen amien o
El núme o de i e aciones de calen amien o, es deci , odas aquellas p e ias que se desca -
an po no pe enece a la dis ibución de la que se p e ende mues ea depende an o de los
alo es iniciales como de la asa de con e gencia de la cadena (hablamos de es o úl imo en el
siguien e apa ado). En mayo medida depende de cuán ce ca es én las dis ibución p opues a
y la es aciona ia. La mane a en la que se suele de e mina el núme o
m
de i e aciones que se
22 Capí ulo 3. Ma ko Chain Mon e Ca lo
debe desca a se suele ealiza median e análisis g á ico de la aza de la cadena. Exis en en la
li e a u a algunos es imado es de
m
analí icos, pe o no suelen e se con demasiada sol u a en las
aplicaciones p obablemen e po que el análisis g á ico p opo ciona la su icien e in o mación.
O a cues ión elacionada es cuándo pa a la cadena pa a ealiza es imaciones Mon e Ca lo.
La medida más usual pa a ello es co e como se ha suge ido an es a ias cadenas en pa alelo, de
al mane a que si con e gen ápidamen e odas hacia una egión de e minada podemos asumi
que ya ha alcanzado su dis ibución es aciona ia y co e es a ez una cadena su icien emen e
la ga pa a ealiza es imaciones.
Gene almen e, las i e aciones pe enecien es al calen amien o se desca an en las es imaciones
de medias y a ianzas pues pueden sesga las.
3.5.3 Análisis de la salida
Como se ha enido explicando median e pinceladas has a aho a, una salida de una simulación
Mon e Ca lo p opo ciona su icien e in o mación pa a:
Realiza un g á ico de líneas de la aza pa a comp oba con e gencia y calen amien o.
Desca a es as úl imas i e aciones en los siguien es análisis.
Es ima medias y a ianzas de la cadena.
Pod íamos incluso es ima in e alos de con ianza, usando las es imaciones de a ianza
aho a mismo omadas o de mane a más sencilla, omando los pe cen iles
α
2
y
1−α
2
de la
simulación.
Rizando un poco el izo, pod iamos es ima las dis ibuciones ma ginales median e alguna
unción núcleo
K(X.i)
. Una elección bas an e co ien e pa a es a unción es la dis ibución
o almen e condicionada de X.i|X.−i, como se explicó en el mues eado de Gibbs.
3.5.4 Tasa de con e gencia
Se dice que una cadena de Ma ko
X
es geomé icamen e e gódica si es ape iódica posi i a
ecu en e y exis e un λ∈(0,1)y una unción V(.) al que se cumple:
∑
j
|Pi j −π(j)| ≤ V(i)λ(3.12)
El
λ
más pequeño pa a el que exis a una unción
V
que sa is aga la condición an e io se
denomina asa de con e gencia. Lo no a emos po
λ∗
. Pa a en ende mejo las implicaciones
de las cadenas geomé icamen e e gódicas ecu imos al análisis espec al. Es o se escapa
ampliamen e del obje i o del abajo, pe o di emos que pa a las cadenas geomé icamen e
e gódigas el p incipal au o alo es
λ0=1
y el es o ( ini os) es án aco ados po el cí culo unidad.
La elocidad a la que con e ge la cadena o
λ∗
es po an o dependien e del segundo au o alo
más g ande.
3.5.5 Es imación de la a ianza
Una de las consecuencias más impo an es de las cadenas e gódicas es que pe mi e la
exis encia de esul ados del ipo eo ema cen al del lími e pa a medias e gódicas de la o ma:
n−E[ (X)] →N(0,σ2)(3.13)
Dicha con e gencia es en dis ibución. Po an o es de i al impo ancia es ima co ec amen e
σ2
. Si bien pueden u iliza se es imado es simples como los p opues os en in eg ación Mon e
Ca lo, se han p opues o algunas al e na i as más so is icadas.
3.5 O as conside aciones 23
Ba ch means
La idea de ás de es e p ocedimien o es co e una cadena de Ma ko en
N=mn
i e aciones,
con nsu icien emen e g ande. No emos po :
Yk=1
n
kn
∑
i=(k−1)n+1
(Xi)(3.14)
con k=1,...,m. De es a mane a, una buena es imación de σ2pod ía eni dada po :
ˆ
σ2≈n
m−1
m
∑
k=1Yk− N(3.15)
Eje cicio 3.5 Una posible implemen ación gene al de ba ch means pod ía se la siguien e:
n=1000
m=5
cadena = exp(n*m, a e = 1/3)
ba ch.means <- unc ion(chain, ,m){
N = mean( (chain))
n = leng h(chain)/m
li = nume ic(m)
o (i in 1:m){
li[i]=mean( (chain[((i-1)*n+1):(i*n)]))
}
sigma = n/(m-1)*sum((li- N)^2)
e u n(sigma)
}
ba ch.means(cadena,iden i y,5)
Es imado es de en ana
O a opción pa a es ima
σ2
es usa la unción de au oco a ianzas mues al. Sin emba go,
es e mé odo p oduce es imado es inconsis en es según aumen a el e a do
i
, po lo que se p opone
una e sión uncada del es imado , que puede exp esa se de la siguien e mane a:
σ2≈ˆ
γ0+2
∞
∑
i
wn(i)ˆ
γi(3.16)
donde los
wn(i)
son pesos e i icando
|wn(i)|<1
y
∑iwn(i) = 1
. Una posible implemen ación
ápida pod ia se la siguien e:
Eje cicio 3.6
cadena = exp(5000, a e = 1/3)
window.es <- unc ion(cadena,weigh s){
au = as.nume ic(ac (cadena,plo =F, ype="co a iance")$ac )
sigma = au [1]+2*sum(pesos*au [2:(leng h(pesos)+1)])
e u n(sigma)
24 Capí ulo 3. Ma ko Chain Mon e Ca lo
}
pesos = ep(1/12,12)
window.es (cadena,pesos)
Dig a os acíclicos
G a o de independencia condicional
Ejemplos de aplicación
4. Modelos g a icos y DAGs
4.1 Dig a os acíclicos
Es e capí ulo se á el más co o dedicado a aplicaciones. Los modelos g á icos se usan cada
ez más ecuencia en modelos bayesianos en los que se aplica MCMC. Las elaciones en e las
a iables del modelo pueden ep esen adas haciendo que los nodos en un g a o ep esen en esas
a iables, y los ejes en e los nodos an e io es ep esen ando la elación (o ausencia) en é minos
de independencia condicional.
Es os g a os, gene almen e se p esen an con una es uc u a je á quica, con aquellos nodos
que eje cen una in luencia más di ec a sob e los da os colocados en la pa e in e io , y según a
disminuyendo dicha in luencia colocando supe io men e. Es os g a os p opo cionan una mane a
ácil de in e p e a la es uc u a condicional del modelo, simpli icando la implemen ación del
algo i mo MCMC en cues ión, indicando qué a iable eje ce in luencia sob e qué o a en su
dis ibución.
Es os g a os se deno an DAGs (Di ec ed Acyclic G aphs). Consis en en una colección de
nodos y ejes di igidos, donde la di ección de es os ejes de e mina el sen ido de la dependencia
en e a iables. Los nodos, a su ez, pueden se de dos ipos: cí culos si se a a de una a iable
desconocida y que po an o hay que es ima simulando, o cuad ados si se a a de una a iable
conocida. Los DAGs se ca ac e izan po se acíclicos, es deci , no podemos ol e al mismo
nodo de pa ida sea cual sea. Un ejemplo de DAG pod ía se el de la igu a 4.1.
Modelos de mix u as ini os
Usando MCMC
Di icul ad en easignación
De e minando el núme o de subpoblacio-
nes
5. Mix u as gaussianas
En es e capí ulo e emos una de las aplicaciones más in e esan es que ienen las écnicas
MCMC. Den o del campo del ap endizaje no supe isado, podemos conside a las mix u as
como pa e de los algo i mos de clus e ing o conglome ados, donde dado un conjun o de da os el
obje i o es desc ibi si el mismo es á compues o de dos o más subpoblaciones bien di e enciadas.
Conc e amen e, el en oque que u iliza emos en es e capí ulo es in en a de e mina si una
mues a p ocede de dos o más dis ibuciones no males. Las mix u as son en ealidad un caso
pa icula de un conjun o de modelos denominados de a iable la en e o ausen e.
5.1 Modelos de mix u as ini os
En la mayo ía de ex os, una densidad de mix u a queda ep esen ada como una suma
ponde ada de unciones de densidad independien es (o una combinación con exa):
k
∑
j=1
pj j(x),;pj≥0;
k
∑
j=1
pj=1 (5.1)
Donde
k
es el núme o de subpoblaciones. En las si uaciones más simples, las dis ibuciones
j
se conocen y lo que in e esa es es ima las can idades
pj
, o en las p obabilidades de que dado un
elemen o de la mues a, pe enezca a una subpoblación u o a. Gene almen e, las dis ibuciones
an e io es pe enecen a una dis ibución pa amé ica, po lo que hay que inclui θ:
k
∑
j=1
pj (x|θj)(5.2)
Dependiendo de la si uación los obje i os de las mix u as pueden se a ios: puede se es ima
la pe enencia a los g upos de las obse aciones
z
(clus e ing), pa a p opo ciona es imado es
de los pa áme os de las dis ibuciones asociadas o incluso es ima el núme o de poblaciones
subyacen es en el modelo.
En cualquie caso, conside emos una mues a alea o ia simple
x= (x1,...,xn)
p oceden e de
un modelo aco de a 5.1. Conside emos la unción de e osimili ud asociada:
34 Capí ulo 5. Mix u as gaussianas
l(θ,p|x) =
n
∏
i=1
k
∑
j=1
pj (xi|θj)(5.3)
A con inuación empezamos el camino pa a in en a en oca es e p oblema a un ma co de
simulación MCMC. Pa a cualquie dis ibución a p io i
π(θ,p)
, la dis ibución a pos e io i de
(θ,p)es á disponible has a cons an e mul iplica i a (como de cos umb e):
π(θ,p|x)∝"n
∏
i=1
k
∑
j=1
pj (xi|θj)#π(θ,p)(5.4)
Pa a en ende mejo cómo amos a mon a el modelo pa a la simulación, conside emos el
en oque de a iable la en e an e io men e mencionado. Pa a cada
xi
, di emos que iene asociada
una a iable la en e
zi
que indica la pe enencia de esa obse ación a una de e minada subpobla-
ción
i
. Una ez en endido es o, comenzamos de iniendo algunas dis ibuciones condicionadas de
mane a adecuada:
zi|p∼Mk(p1,..., pk)(5.5)
En p ime luga , la a iable la en e
zi
condicionada a
p
se dis ibui á según una dis ibución
mul inomial de pa áme os de inidos po p.
xi|zi,θ∼ (.|θzi)(5.6)
Las obse aciones condicionadas a la subpoblación que pe enecen y a los pa áme os de la
subpoblación que pe enecen se dis ibuyen según
. En nues o caso és a unción de densidad
se á la de una no mal, pe o el modelo es ex ensible a cualquie o a dis ibución (con mayo o
meno di icul ad). Con es a es uc u a, podemos de ini de nue o la unción de e osimili ud:
l(θ,p|x,z) =
n
∏
i=1
pzi (xi|θzi)(5.7)
Donde
z= (z1,...,zn)
.Po an o la unción a pos e io i de inida an e io men e pasa ía a ene
la siguien e o ma, sus i uyendo:
π(θ,p|x,z)∝"n
∏
i=1
pzi (xi|θzi)#π(θ,p)(5.8)
Aho a amos a ealiza una expansión con enien e pa a la exp esión del modelo. Deno emos
po
Z=1,...,kn
, el conjun o de los
kn
alo es posibles del ec o
z
.
Z
se puede descompone
i ialmen e usando unión de conjun os
Z=∪τ
j=1Zj
. Siguiendo con es e azonamien o, dado un
ec o de núme o de asignaciones a cada g upo
(n1,...,nk)
, de inimos una se ie de conjun os de
pa ición:
Zj=(z:
n
∑
i=1
Izi=1,...,
n
∑
i=1
Izi=k)(5.9)
Es un poco en e esado segui lo siguien e: lo an e io ep esen a odas las posibles asigna-
ciones dado un núme o de asignaciones
(n1,...,nk)
. Nomb amos a cada uno de es os conjun os
de pa ición
j=j(n1,...,nk)
, usando el ó den lexicog á ico. U ilizando es o, la dis ibución a
pos e io i puede esc ibi se de o ma ce ada:
π(θ,p|x) = ∑
z∈Z
π(θ,p|x,z) =
τ
∑
i=1∑z∈Ziw(z)π(θ,p|x,z)(5.10)
5.2 Usando MCMC 35
Hemos en e esado an o la unción a pos e io i con dos inalidades: la p ime a encon a
una exp esión ce ada pa a la unción a pos e io i, y lo segundo es pa a comp oba que
w(z)
ep esen a la p obabilidad ma ginal pos e io de una asignación a
z
condicionada a
x
. Es po ello
que un es imado Bayes de la dis ibución an e io es inmedia o:
E[θ,p|x] =
τ
∑
i=1∑
z∈Zi
w(z)E[θ,p|x,z](5.11)
Es a descomposición an poco na u al que hemos desa ollado iene bas an e u ilidad en
o os mé odos que p eceden a las simulaciones MCMC que amos a desa olla a pa i de aho a,
conc e amen e al algo i mo EM.
5.2 Usando MCMC
Llegados a es e pun o, nos plan eamos cómo aplica lo que hemos is o a las mix u as. La
dis ibución o almen e condicionada de z|xes á disponible sal o cons an e mul iplica i a:
π(z|x,θ,p)∝
n
∏
i=1
pzi (xi|θzi)(5.12)
A con inuación al a de e mina una se ie de dis ibuciones conjugadas pa a las pos e io i.
Supongamos que
p
y
θ
son independien es a p io i, en onces, dado
z
los ec o es
p
y
x
son
independien es:
π(p|z,x)∝π(p) (z|p) (x|z)∝π(p|z)(5.13)
De mane a simila ,
θ
es ambién independien e a pos e io i de
p|z
y
x
, con densidad
π(θ|z,x)
.
Aquí ya podemos aplica el mues eo de Gibbs, simulando po una pa e
z
condicionado sob e
(p,θ)
y sob e los da os (y al e és). De es a mane a, nues o algo i mo end ía que hace lo
siguien e:
1. Inicio: Escoge p(0)yθ(0)cualesquie a.
2. ( ++)
Pa a cada elemen o de la mues a, gene a
z( )
i
al que
P(zi=j|θ,p)∝p( −1)
j (xi|θ( −1)
j)
.
3. Gene a p( )según π(p|z( ))
4. Gene a θ( )según π(θ|z( ),x)
Pa a simula
p
, lo suyo es escoge una dis ibución conjugada en
z
. La complexidad de la
simulación de los
θj
depende á de la
con la que es emos abajando (en nues o caso, como
son no males y exis en dis ibuciones conjugadas pa a ambos pa áme os no hay p oblema). Pa a
lo p ime o que hemos mencionado, los
zi
se dis ibuyen como hemos dicho an e io men e según
una mul inomial
M(p1,..., pk)
, que pe mi e encon a ácilmen e una dis ibución conjugada en
z. De es a mane a, la dis ibución a p io i de p∼D(γ1,...,γk)se á una Di ichle , con densidad:
Γ(γ1+...+γk)
Γ(γ1)...Γ(γk)pγ1
1...pγk
k(5.14)
De es a mane a, la dis ibución p|zse á Di ichle ambién:
p|z∼D(n1+γ1,...,nk+γk)(5.15)
Eje cicio 5.1
R no ae una unción de base pa a simula a iables alea o ias según una
dis ibución Di ichle . Implemen a una es ácil a a és de alo es alea o ios p oceden es de
dis ibuciones gamma o be a. Implemen emos uno de mane a ápida.
36 Capí ulo 5. Mix u as gaussianas
#!/us /local/en RSc ip
di ichle <- unc ion(n=1,pa ams=c(1,1)){
k = leng h(pa ams)
a ay = ma ix(NA, n ow = n, ncol = k)
o (i in 1:n){
suppo = gamma(k,shape=pa ams)
a ay[i,] = suppo /sum(suppo )
}
e u n(a ay)
}
Pa a las mix u as no males, las µj|σj,z,xson independien es, con dis ibución:
µj|σj,z,x∼N εj(z),σ2
j
nj+lj!(5.16)
donde:
εj(z) = ljεj+njxj(z)
lj+nj
(5.17)
De mane a simila :
σ2
j|x,z∼IG(0,5( j+nj),0,5sj(z)) (5.18)
donde:
sj(z) = s2
j+njs2
j(z)+ ljnj
lj+nj
(εj−xj(z))2(5.19)
Después de an a eo ía, a siendo ho a de que implemen emos nues o mues eo de Gibbs
o ien ado a mix u as gaussianas con
k
componen es. En el código siguien e, asumimos que
lj=1 j=20:
Eje cicio 5.2
#!/us /bin/en RSc ip
# Gibbs Sample (Gaussian mix u es)
gibbsMix u e <- unc ion(da a,k,max_i e =1000){
n = leng h(da a)
mu = mean(da a)
sig = a (da a)
# P ealloca ing da a
z = ep(NA,n)
g oup_size = ep(NA, k)
g oup_x = ep(NA, k)
g oup_sum = ep(NA, k)
mu_ma = sig_ma = pj_ma = ma ix(NA, n ow=max_i e , ncol=k)
mu_ma [1,] = ep(mu,k)
5.2 Usando MCMC 37
pj_ma [1,] = ep(1,k)/k
sig_ma [1,] = ep(sig,k)
## Possible chunk o likelihood
#######
### Main loop
o (i in 2:max_i e ){
# z es ima ion
o (j in 1:n){
p = pj_ma [i-1,]*dno m(da a[j],mu_ma [i-1,],sq (sig_ma [i-1,]))
z[j] = sample(1:k, size = 1, p ob = p)
}
# es o pa ame e s
o (j in 1:k){
g oup_size[j] = leng h(which(z == j))
g oup_x[j] = sum(as.nume ic(z == j)*da a)
g oup_sum[j] = sum(as.nume ic(z==j)*(da a-g oup_x[j]/g oup_size[j])^2)
}
mu_ma [i,] = no m(k, (mean(da a) + g oup_x)/(g oup_size + 1),
sq (sig_ma [i-1,]/g oup_size + 1))
sig_ma [i,] = 1/ gamma(k, .5*(20+g oup_size),
a (da a) + .5* g oup_sum +
.5* g oup_size/(g oup_size + 1)*(mean(da a) -
g oup_x/g oup_size)^2)
pj_ma [i,] = di ichle (pa ams = g oup_size + 0.5)
}
class = se Class("MCMC mix u e",
slo s=c(mu="ma ix", sig="ma ix", p="ma ix", g oup="nume ic"))
es = class(mu = mu_ma , sig = sig_ma , p = pj_ma , g oup = z)
e u n( es)
}
38 Capí ulo 5. Mix u as gaussianas
Eje cicio 5.3
Usamos un ma co de da os clásico como ai h ul, incluído en R pa a pone a
p ueba nues o algo i mo. En p ime luga ep esen amos los da os pa a hace nos una idea
isual de si un modelo de mix u as es adecuado.
da a( ai h ul)
his ( ai h ul$wai ing, b eaks = 20, col="ligh blue",
eq=F, main="Mix u a Gaussiana")
Apa en emen e pa ece que podemos modela dicho conjun o de da os según una mix u a
gaussiana de k=2 componen es. Aplicamos nues o algo i mo y omamos medias:
se .seed(93)
mix u e = gibbsMix u e( ai h ul$wai ing,k = 2,max_i e = 1000)
(mu = apply(mix u e@mu, 2, mean))
[1] 55.11955 80.08556
(sig = apply(mix u e@sig, 2, mean))
[1] 38.28160 33.79188
(p = apply(mix u e@p, 2, mean))
[1] 0.367852 0.632148
Finalmen e, podemos supe pone sob e el his og ama an e io la densidad de nues a
mix u a:
cu e(p[1]*dno m(x,mu[1],sq (sig[1])),add=T, lwd=2)
cu e(p[2]*dno m(x,mu[2],sq (sig[2])),add=T, col=" ed", lwd=2)
5.3 Di icul ad en easignación 39
Los esul ados an e io es pueden lle a a una alsa sob e alo ación del modelo es adís ico
en e manos. Reco demos que el mues eado de Gibbs es á condicionado a los alo es iniciales
que escojamos pa a la cadena. El condicionamien o sob e
z
implica que las cadenas y po an o
an o
θ
como
p
son incapaces de ealiza cambios d ás icos sob e sí mismos y sob e las asigna-
ciones en la siguien e i e ación. O as écnicas an e io es (como el algo i mo EM) ambién son
sensibles a no ealiza cambios d ás icos sob e asignaciones y po an o son igual de ágiles con
espec o a es e sen ido.
Una posible solución es u iliza el mues eo de Gibss en conjun o con alguna o a écni-
ca MCMC, como Me opolis-Has ings. Pod íamos usa en es e sen ido una p obabilidad de
acep ación con una adecuada dis ibución p opues a del ipo:
π(θ0,p0|x)
π(θ,p|x)
q(θ,p|θ0,p0)
q(θ0,p0|θ,p)(5.20)
5.3 Di icul ad en easignación
Como se ha is o en la sección an e io , al mues eo de Gibbs le cues a ealiza cambios
d ás icos sob e easignaciones en
z
. Una ca ac e ís ica cu iosa en los modelos de mix u as es que
ca ece de un o den, es deci , dada una mix u a compues a po dos densidades, es inco ec o deci
que una de las dos es el p ime componen e de la misma o el segundo. De mane a más écnica
decimos que los pa áme os θjno son ma ginalmen e iden i icables.
Po si ue a poco, en una mix u a de
k
componen es se puede p oba que el núme o de modas
en la e osimili ud es del o den de
O(k!)
. Es o es ácilmen e comp obable ya que si
(θ,p)
es
un máximo local en la e osimili ud, en onces una pe mu ación de es os pa áme os ambién
sigue siéndolo. Es más, si usamos una dis ibución a p io i sob e
(θ,p)
la cual es in a iable bajo
pe mu ación de los índices, odas las dis ibuciones a pos e io i ma ginales son idén icas, lo cual
implica ía que los es imado es Bayes son los mismos.
Es e p oblema sob e easignaciones iene dis in as soluciones. Po un lado, una solución
40 Capí ulo 5. Mix u as gaussianas
sencilla pod ía lle a se a cabo imponiendo una es icción de o den en los pa áme os (po
ejemplo, o denando las medias). Es o se educe en su o alidad a unca la dis ibución a p io i:
π(θ,p)Iµ1≤...≤µk(5.21)
Con es a educción del espacio pa amé ico pueden ocu i sucesos inespe ados, ya que las
dis ibuciones a pos e io i no ienen po qué espe a su opología, de al mane a que cuando
ponemos simulaciones a unciona los es imado es no ienen po qué acaba en una de las
k!
mencionadas an e io men e sino en una zona de silla de baja p obabilidad. No obs an e, si bien
es a es icción puede supone un cambio de endimien o en las simulaciones, és as no ienen po
qué ealiza se du an e la misma sino después. En es e sen ido, si la es icción es sob e el o den
en las medias, una ez que la simulación ha acabado, los componen es se pueden easigna de
acue do con es e o den.
Es o se puede ealiza de la siguien e mane a: dada una mues a de amaño
M
de una
simulación MCMC, de inimos el es imado máximo a pos e io i como:
i∗=a gmaxi=1,...,Mπ{(θ,p)(i)|x}(5.22)
En palab as, es el alo simulado que maximiza la unción de densidad a pos e io i. Es a,
como ya sabemos, no es necesa io que incluya la cons an e de in eg ación. Sin emba go, es
bas an e p oblable que es e alo se encuen e en la ecindad de una de las
k!
posibles modas.
Usa emos es e alo óp imo (MAP) como pi o e, en el sen ido de que o dena emos el es o de
las o as i e aciones con espec o a dicha moda.
En luga de escoge es e eo den según una dis ancia euclídea en el espacio pa amé ico, la
de inimos en el espacio p obabilís ico de las asignacioens. Deno emos po
Gk
el conjun o de las
kposibles pe mu aciones y τ∈Gk. Minimizamos en τuna dis ancia de en opía:
h(i,τ) =
n
∑
=1
k
∑
j=1
P(z =j|θ(i∗),p(i∗))×log P(zj=k|θ(i∗),p(i∗))
P(z =j|τ(θ(i),p(i)))!(5.23)
En de ini i a, podemos de ini nues o algo i mo de eo denamien o pi o al como:
En odas las i e aciones, calcula τi=a gmin h(i,τ)
(θ(i),p(i)) = τi{(θ(i),p(i))}
Con es e eo denamien o, en la mayo ía de i e aciones las asignaciones quedan easignadas
a la misma moda, y po an o eliminamos el p oblema de la iden i icación. Después de es e
eo denamien o, los es imado es Mon e-Ca lo siguen siendo los habi uales. Realicemos la
implemen ación de es e eo denamien o en R:
Eje cicio 5.4
Suponemos que enemos una salida de la unción gibbsMix u e, la cual usamos
en es a unción:
pi o alReo <- unc ion(gibbsMix){
mu = gibbsMix@mu
sig = gibbsMix@sig
p = gibbsMix@p
logpos = gibbsMix@logpos
n = gibbsMix@n
k = gibbsMix@k
da a = gibbsMix@da a
5.3 Di icul ad en easignación 41
max_i e = gibbsMix@max_i e
map_indices = o de (logpos , dec easing = T)[1]
map = lis (mu = mu[map_indices,], sig= sig[map_indices,],
p = p[map_indices,])
lili = ma ix(NA, n, k)
alloca ion = ma ix(NA, n, k)
o ( in 1:n){
lili[ ,] = map$p*dno m(da a[ ], mean = map$mu, sd=sq (map$sig))
lili[ ,] = lili[ ,]/sum(lili[ ,])
}
o de ed_mu = ma ix(NA, ncol=k, n ow=max_i e )
o de ed_sig = ma ix(NA, ncol=k, n ow=max_i e )
o de ed_p = ma ix(NA, ncol=k, n ow=max_i e )
equi e(combina )
pe ma = pe mn(k)
o ( in 1:1000){
en opies = ep(0, ac o ial(k))
o (j in 1:n){
alloca ion[j,] = p[ ,]*dno m(da a[j], mean=mu[ ,], sd=sq (sig[ ,]))
alloca ion[j,] = alloca ion[j,]/sum(alloca ion[j,])
o (i in 1: ac o ial(k)){
en opies[i] = en opies[i] + sum(lili[j,]*log(alloca ion[k,pe ma[[i]]]))
}
}
bes _o de ing = o de (en opies, dec easing=T)[1]
o de ed_mu[ ,] = mu[ ,pe ma[[bes _o de ing]]]
o de ed_sig[ ,] = sig[ ,pe ma[[bes _o de ing]]]
o de ed_p[ ,] = p[ ,pe ma[[bes _o de ing]]]
}
es = lis (mu=o de ed_mu, sig=o de ed_sig, p=o de ed_p)
e u n( es)
}
48 Capí ulo 6. Análisis de imagen
#include <RcppA madillo.h>
#include <RcppA madilloEx ensions/sample.h>
#include <ma h.h> /* exp */
// [[Rcpp::depends(RcppA madillo)]]
using namespace Rcpp;
// [[Rcpp::expo ]]
in nei4(Nume icMa ix x, in a, in b, in col){
in n = x.n ow(), m = x.ncol();
in nei = 0;
a = a-1;
b = b-1;
i (a != 0){
i (x(a-1,b) == col){nei++;}
}
i (b != 0){
i (x(a,b-1) == col){nei++;}
}
i (a != (n-1)){
i (x(a+1,b) == col){nei++;}
}
i (b != (m-1)){
i (x(a,b+1) == col){nei++;}
}
e u n nei;}
// [[Rcpp::expo ]]
Nume icMa ix isingSample (Nume icMa ix x, in max_i e , double be a){
in n = x.n ow(), m = x.ncol();
in i;
in k;
in l;
Nume icVec o ows(n);
ows = Rcpp::seq(1,n);
Nume icVec o cols(m);
cols = Rcpp::seq(1,m);
Nume icVec o pos(2);
pos = Rcpp::seq(0,1);
Nume icVec o pe mu _ ows(n);
Nume icVec o pe mu _cols(m);
Nume icVec o p ob(2, 0.5);
in n0;
in n1;
6.1 Recons ucción de imágenes 49
in pos1;
in pos2;
double es;
o (i = 0; i < max_i e ; i++){
pe mu _ ows = RcppA madillo::sample( ows,n,0);
pe mu _cols = RcppA madillo::sample(cols,m,0);
o (k = 1; k <= n; k++){
o (l = 1; l <= m; l++){
n0 = nei4(x, pe mu _ ows[k-1], pe mu _cols[l-1], 0);
n1 = nei4(x, pe mu _ ows[k-1], pe mu _cols[l-1], 1);
p ob[0] = exp(be a*n0);
p ob[1] = exp(be a*n1);
pos1 = pe mu _ ows[k-1]-1;
pos2 = pe mu _cols[l-1]-1;
double suma = p ob[0] + p ob[1];
p ob[0] = p ob[0]/suma;
p ob[1] = p ob[1]/suma;
es = RcppA madillo::sample(pos, 1, 0, p ob)[0];
// es = Rcpp::Func ion(sample(pos,1,0,p ob));
Rcpp::Func ion(gc());
x(pos1, pos2) = es;
}
}
Rcpp::Func ion(gc());
}
e u n x;
}
Un asun o de in e és es comp oba cómo se compo a el algo i mo dependiendo del alo
de βescogido. Podemos ealiza una p ueba empí ica ápida:
be as = seq(0.2, 1.8, by = 0.3)
50 Capí ulo 6. Análisis de imagen
pa (m ow=c(3,2))
se .seed(39)
x = ma ix(sample(c(0,1),50*50, ep=T), n ow=50)
o (be a in be as){
se .seed(93)
y = isingSample (x, 1000, be a)
image(y, main=pas e("be a=",be a))
}
Que p oduce el siguien e esul ado:
10 20 30 40 50
10 20 30 40 50
be a= 0.2
1:50
1:50
10 20 30 40 50
10 20 30 40 50
be a= 0.5
1:50
1:50
10 20 30 40 50
10 20 30 40 50
be a= 0.8
1:50
1:50
10 20 30 40 50
10 20 30 40 50
be a= 1.1
1:50
1:50
10 20 30 40 50
10 20 30 40 50
be a= 1.4
1:50
1:50
10 20 30 40 50
10 20 30 40 50
be a= 1.7
1:50
1:50
Como se puede comp oba , cuan o mayo es
β
, el modelo iende a concen a se bajo
dis ibuciones de un solo colo , y po an o más homogénea la imagen.
6.1 Recons ucción de imágenes 51
6.1.3 Modelo de Po s
El modelo de Po s es una gene alización na u al al modelo de Ising cuando la imagen iene
más de dos colo es, po ejemplo
G
. Es a ez no amos po
ni,g
al núme o de ecinos de
i
con
colo
g
, es deci ,
ni,g=∑j∼iIxj=g
. La dis ibución o almen e condicionada de
xi
se escoge de
mane a idén ica como:
π(xi=g|x−i)∝exp(βni,g)(6.5)
La densidad conjun a del modelo de Po s, como ya espe amos iene la o ma:
π(x)∝exp β∑
j∼i
Ixj=xi!(6.6)
Aquí nos encon amos con el mismo p oblema que en la sección an e io , y es que no
conocemos el alo de
β
, y po an o no podemos calcula la unción de e osimili ud. Po o a
pa e, con un
β
g ande, al y como pasaba an es, es di ícil que un alo se ac ualice y se educe la
elocidad de con e gencia. Se p opone una modi icación basada en pasos Me opolis-Has ings
que ue za ac ualizaciones en cada paso.
1. Inicialización: Pa a cada i∈I, gene a x(0)
i∼UD(1,...,G)
2. I e ación >1. Gene a u, una pe mu ación alea o ia de los elemen os de I.
3. Pa a ldesde 1 has a |I|:
gene a xu,l∼UD1,2,...,x( −1)
u,l−1,x( −1)
u,l+1,...,G
gene a núme o de ecinos n( )
ul,gy
gene a p obabilidad de acep ación pl=maxexp(βnul,xul )/exp(βnul,xul )( ),1
si se acep a, sus i ui xul po x( )
ul
Eje cicio 6.2
De nue o seguimos con la misma mecánica aplicada en el modelo an e io .
Una implemen ación ine icien e del algo i mo en R pod ía se la siguien e:
po sSample <- unc ion(x, num_col=2, max_i e =1000, be a){
n = dim(x)[1]; m = dim(x)[2]
o (i in 1:max_i e ){
pe mu = sample(1:(n*m))
ca ("I e acion", i, " n")
o (k in 1:(n*m)){
xcu = x[pe mu [k]]
a = (pe mu [k]-1)%%n + 1
b = (pe mu [k]-1)%/%n + 1
x ilde = sample((1:num_col)[-xcu ],1)
p ob = be a*(nei4(x,a,b,x ilde)-nei4(x,a,b,xcu ))
i (log( uni (1))<p ob){
x[pe mu [k]] = x ilde
}
}
52 Capí ulo 6. Análisis de imagen
}
e u n(x)
}
Y una a ios ó denes más e icien e, en
C++
.
#include <RcppA madillo.h>
#include <RcppA madilloEx ensions/sample.h>
#include <ma h.h>
// [[Rcpp::depends(RcppA madillo)]]
using namespace Rcpp;
// [[Rcpp::expo ]]
Nume icMa ix po sSample (Nume icMa ix x, in num_col, in max_i e , double be a){
in n = x.n ow(), m = x.ncol();
in i;
in xcu ;
in k, j;
in a, b;
in x ilde;
double p ob;
Nume icVec o pixels(n*m);
pixels = Rcpp::seq(1, n*m);
o (i = 0; i < max_i e ; i++){
Nume icVec o pe mu (n*m);
pe mu = RcppA madillo::sample(pixels,n*m,1);
o (k = 0; k < (n*m); k++){
Nume icVec o colo s(num_col);
colo s = Rcpp::seq(1, num_col);
xcu = x[pe mu [k] - 1];
a = emainde ((pe mu [k]-1),n) + 1;
b = (pe mu [k]-1)/n + 1;
colo s.e ase(xcu - 1);
x ilde = RcppA madillo::sample(colo s, 1, 0)[0];
p ob = be a*(nei4(x, a, b, x ilde) - nei4(x, a, b, xcu ));
double alea = uni (1)[0];
i (log(alea) < p ob){
x[pe mu [k] - 1] = x ilde;
}
}
6.1 Recons ucción de imágenes 53
}
e u n x;
}
De mane a simila a lo que hicimos pa a el modelo de Ising, podemos comp oba de
mane a empí ica el compo amiendo del modelo eniendo en cuen a el alo de β.
be as = seq(0.2, 1.4, by = 0.4)
pa (m ow=c(2,2))
se .seed(39)
x = ma ix(sample(1:4,50*50, ep=T), n ow=50)
o (be a in be as){
se .seed(93)
y = po sSample (x, num_col=4, 1000, be a)
image(x=1:50, y=1:50,z=y, main=pas e("be a=",be a))
}
10 20 30 40 50
10 30 50
be a= 0.2
1:50
1:50
10 20 30 40 50
10 30 50
be a= 0.6
1:50
1:50
10 20 30 40 50
10 30 50
be a= 1
1:50
1:50
10 20 30 40 50
10 30 50
be a= 1.4
1:50
1:50
De mane a simila , cuan o mayo es el pa áme o, más endencia se iene hacia imágenes
más homogéneas.
6.1.4 Sob e la cons an e de in eg ación
En las secciones an e io es hemos supues o que
β
, y que po an o
Z(β)
e an conocidos.
Como cabe espe a , es o no se suele p oduci en la mayo ía de las ocasiones. El manejo de esa
cons an e de in eg ación ha dado luga a mucha li e a u a, pues es un p oblema di ícil. Aquí
e emos una me odología ela i amen e sencilla pa a a e igua la, denominada mues eo po
caminos. Es a écnica es á basada en una ep esen ación de la de i ada de dicha cons an e:
54 Capí ulo 6. Análisis de imagen
dZ(β)
dβ=∑
x
S(x)exp(βS(x)) (6.7)
Es a de i ada puede exp esa se de mane a con enien e como una espe anza:
dZ(β)
dβ=Z(β)∑
x
S(x)exp(βS(x)
Z(β)=Z(β)E[S(x)] (6.8)
y po an o, omando loga i mo an es:
dlogZ(β)
dβ=E[S(x)] (6.9)
Es a ep esen ación pe mi e que podamos exp esa el a io Z(β1
Z(β0)como una in eg al:
log(Z(β1)/(Z(β0)) = Zβ1
β0
E[S(x)]dβ(6.10)
Es a úl ima ecuación es lo que se denomina iden idad del mues eo po caminos. Con es a
exp esión, podemos ecu i a p ocedimien os de simulación es ánda es pa a su ap oximación.
La in eg al puede ap oxima se po mé odos numé icos, y pa a un alo de
β
,
E[S(x)] = (β)
puede simula se usando el modelo de Po s, po ejemplo. Finalmen e ap oximamos
(β)
po
una unción lineal a ozos pa a in eg a ácilmen e.
Eje cicio 6.3
Realicemos una ap oximación de
(β)
usando R. Pa a un alo de e minado
de β:
pa hSampling <- unc ion(x, ncol=2, max_i e = 10^3, be a){
n = dim(x)[1]; m = dim(x)[2]
S = 0
o (i in seq_len(max_i e )){
i (i%%100==0){
ca ("I e ación ",i," n")
}
s = 0
ows = sample(1:n)
cols = sample(1:m)
o (k in seq_len(n)){
o (l in seq_len(m)){
n0 = nei4(x, ows[k], cols[l], x[ ows[k], cols[l]])
col = sample(1:ncol, 1)
n1 = nei4(x, ows[k], cols[l], col)
i (log( uni (1))< (be a*(n1-n0))){
x[ ows[k], cols[l]] = col
n0 = n1
}
s = s+n0
}
}
i (2*i > max_i e ){
6.1 Recons ucción de imágenes 55
S=S+s
}
}
e u n(2*S/max_i e )
}
Y pa a gene a una secuencia, en in e alos an pequeños como se quie a:
gene a e be a <- unc ion(x, max_i e =10^3, ncol=2){
sea ch = seq(0.1, 2, by=0.1)
Z = seq_along(sea ch)
i = 1
o (be a in sea ch){
ca ("Be a: ", be a, " n")
Z[i] = pa hSampling(x, ncol, max_i e , be a)
i = i+1
}
plo (sea ch, Z, main=" (be a) app ox.", ype="l")
e u n(Z)
}
El p oblema en e iciencia es que pa a imágenes de jugue e nues o algo i mo es á
bien, pe o como enimos acos umb ando en es e capí ulo, se ía buena idea implemen a
pa hSampling
en
C++
.
double pa hSampling(Nume icMa ix x, in ncol, in max_i e ,
double be a){
in n = x.n ow(), m = x.ncol();
double S = 0.0, s = 0.0;
in i, k, l, n0, n1, col;
Nume icVec o ows(n);
ows = Rcpp::seq(1, n);
Nume icVec o cols(n);
cols = Rcpp::seq(1, m);
Nume icVec o colo s(ncol);
colo s = Rcpp::seq(1, ncol);
o (i = 0; i < max_i e ; i++){
s = 0.0;
ows = RcppA madillo::sample( ows, n, 0);
cols = RcppA madillo::sample(cols, m, 0);
o (k = 0; k < n; k++){
56 Capí ulo 6. Análisis de imagen
o (l = 0; l < n; l++){
n0 = nei4(x, ows[k], cols[l], x( ows[k]-1, cols[l]-1));
col = RcppA madillo::sample(colo s, 1, 0)[0];
n1 = nei4(x, ows[k], cols[l], col);
i (log( uni (1)[0]) < (be a*(n1-n0))){
x( ows[k]-1, cols[l]-1) = col;
n0 = n1;
}
s = s + n0;
}
}
i (2*i > max_i e ){
S = S + s;
}
}
e u n 2*S/max_i e ;
}
6.1.5 Segmen ación de imágenes
Después de odas es as secciones ya enemos odos los ma e iales disponibles pa a la a ea
que que íamos ealiza . Conside a emos las imágenes como obje os es adís icos, pe o no como
al, sino que conside a emos que enemos una imagen con cie a dis o sión
y
, es deci que el colo
de g is (o cualquie o a escala, po con eniencia) se p esen a con cie a pe u bación. El obje i o
de la segmen ación de imágenes es, dada es a imagen dis o sionada, ealiza un p ocedimien o
de ’clus e ing’ de píxeles a endiendo únicamen e a la es uc u a de dependencia espacial de la
es uc u a.
Es a dependencia espacial ’ eal’ de píxeles se no a po
x
, donde gene almen e
y
puede oma
alo es eales posi i os y
x
solo disc e os, po con eniencia. Es amos pues in e esados en co-
noce la dis ibución a pos e io i de
x
dada
y
, es deci
π(x|y)∝ (y|x)π(x)
. En es a exp esión
(y|x)
ep esen a la e osimili ud de los da os y la unción link en e la imagen o iginal y su
clasi icación, mien as que
π(x)
ep esen a nues a in o mación a p io i (o p opiedades deseadas)
sob e el compo amien o de la imagen ’ eal’.
Gene almen e es a dis ibución a p io i suele se el modelo de Po s ya es udiado con
G
ca ego ías:
π(x|β) = 1
Z(β)exp β∑
i∼j
Ixj=xi!(6.11)
Dado
x
, y po acilidad compu acional, asumimos que las obse aciones en
y
son a iables
independien es no males. Es o se hace po que es más ácil pa ame iza que una dis ibución
mul inomial que ome alo es en 0,...,255:
6.1 Recons ucción de imágenes 57
(y|x,σ2,µ1,...,µG) = ∏
i∈I
1
(2πσ2)(1/2)exp−1
2σ2(yi−µxi)(6.12)
Pa a los pa áme os de es as no males se suelen u iliza dis ibuciones a p io i uni o mes:
β∼U(0,2)(6.13)
µ∼U(0≤µ1≤... ≤µG≤255)(6.14)
π(σ2)∝σ−2I[0,∞)(σ2)(6.15)
Es a úl ima no es más que una dis ibución uni o me en el loga i mo de
σ
. Los
µi
no ienen
po qué o dena se, pe o así e i an el p oblema del eo denamien o pi o al explicado en el capí ulo
de las mix u as. La dis ibución conjun a a pos e io i es la siguien e:
π(x,β,σ2,µ|y)∝π(β,σ2,µ)×1
Z(β)exp β∑
j∼i
Ixj=xi!
×∏
i∈I
1
(2πσ2)1/2exp−1
2σ2(yi−µxi)2(6.16)
Comenzamos a con inuación a cons uí las di e sas dis ibuciones o almen e condicionadas
pa a el abajo del mues eo de Gibbs. La de xi:
P(xi=g|y,β,σ2,µ)∝exp β∑
i∼j
Ixj=g−1
2σ2(yi−µg)2!(6.17)
Aquí puede comp oba se que, una ez que
x
es conocido, los elemen os de dis in a ca ego ía
se sepa an y el es o de pa áme os pueden simula se independien emen e condicionados a
x,y
y
σ2
. Si no amos po
ng=∑i∈IIxi=g
y
sg=∑i∈IIxi=gyi
, la dis ibución o almen e condicionada de
µg
es una dis ibución no mal uncada en
[µg−1,µg+1]
(po azones ob ias,
µ0=0,µg+1=255
),
de media
sg/ng
y a ianza
σ2/ng
. La dis ibución condicionada de
σ2
es una gamma in e sa con
pa áme os |I|2/2 y ∑i∈I(yi−µxi)2/2.
Finalmen e, la dis ibución o almen e condicionada de βes aquella que:
π(β|y)∝1
Z(β)exp β∑
j∼i
Ixj=xi!(6.18)
.
Eje cicio 6.4
Una implemen ación ápida en
R
de la segmen ación de imágenes bayesiana
pod ía se la siguien e:
bay <- unc ion (y, ncol=2, max_i e =10^3)
{numb = dim(y)[1]
x=0*y
mu = ma ix(0, max_i e , 6)
sigma2 = ep(0, max_i e )
mu[1, ] = c(35, 50, 65, 84, 92, 120)
sigma2[1] = 20
be a = ep(1, max_i e )
64 Capí ulo 7. Impu ación múl iple
Es os son solo algunos de los p oblemas que pod ían su gi en nues o análisis. An e io men-
e, exis ía ambién un mé odo de impu ación múl iple basado en MCMC que los da os enían
de inidos po algún modelo de densidad mul i a ian e
P(Y|θ)
. La di icul ad de es e mé odo es
e iden e en el momen o que enemos un ma co de da os con a iables de dis in o ipo. Auna
odas esas a iables bajo un mismo modelo eó ico mul i a ian e (como po ejemplo hace el
algo i mo EM y la no mal mul i a ian e) es po un lado i eal y po o o di ícil.
FCS es un mé odo más na u al en el sen ido de que no pa imos de una densidad mul i a ian-
e. En su luga la de inimos implíci amen e especi icando una densidad uni a ian e pa a cada
a iable del es ilo
P(yj|y−j,θj)
. Es a densidad es la usamos pa a impu a
yaus
j
dado
y−j
median e
algún ipo de eg esión (gene almen e lineal o logís ica) sob e los casos
yobs
j
. Se ealizan an os
paseos sob e odas las densidades como i e aciones del algo i mo sean necesa ias.
Las en ajas de FCS sob e o os mé odos son e iden es: e i a ene que especi ica di ec a-
men e una dis ibución mul i a ian e pe mi e mucha lexibilidad a la ho a de abaja los da os.
Es deci , con e imos un p oblema
k
dimensional en
k
p oblemas unidimensionales. Po o a
pa e, la idea de especi ica un modelo de impu ación di e en e pa a cada a iable es bas an e
na u al. No obs an e, FCS ambién cuen a con algunas di icul ades. El ene que especi ica un
modelo di e en e y adecuado pa a cada a iable puede se po un lado pesado pa a el in es igado ,
y po o o es compu acionalmen e mucho más exigen e. Po o a pa e, e alua la calidad de
las impu aciones puede esul a di ícil ya que la densidad conjun a eó ica no iene po qué exis i .
En cualquie caso, comencemos a de ini el mé odo. Supongamos que enemos
Y= (Y1,...,Yk)
un ec o de
k
a iables alea o ias con una dis ibución
P(Y|θ)
. Asumimos que la an e io dis i-
bución queda o almen e de inida po
θ
. El p ocedimien o gene al comp ende ía los siguien es
pasos:
1. De e mina la dis ibución a pos e io i p(θ|yobs)de θusando los da os obse ados yobs.
2. Mues ea un alo θ∗de p(θ|yobs)
3.
Po úl imo mues ea un alo
y∗
de
p(yaus|yobs,θ=θ∗)
, la dis ibución condicional a
pos e io i de yaus dado θ∗.
Gene almen e es e p ocedimien o es pesado y solamen e se epi en los pasos 2 y 3 un núme o
modes o de eces (gene almen e suelen se su icien es en e 5 y 10 i e aciones). Lógicamen e
desde el ma co eó ico uni a ian e es o pa ece ácil, pe o si
Y
es mul i a ian e (caso
k>1
) lo
lógico es u iliza un ma co de simulación basado en el mues eo de Gibbs. El paso ealmen e
complicado es el p ime o, donde necesi amos de e mina la dis ibución a pos e io i
p(θ|yobs)
.
En nues o caso, aplicamos el mues eo de Gibbs mues eando de dis ibuciones condicionales
de la o ma:
P(Y1|Y−1,θ1)(7.1)
... (7.2)
P(Yk|Y−k,θk)(7.3)
Es deci , u ilizamos los pa áme os
θi
solo en las densidades uni a ian es donde in e ienen
espec i amen e. En o as palab as, es os pa áme os no se ían necesa iamen e el p oduc o de
una dis ibución conjun a
P(Y|θ)
. De mane a más cla a, una i e ación del mues eo de Gibbs
comp ende ía los siguien es pasos:
7.1 Impu ación usando ecuaciones encadenadas 65
θ∗
1∼P(θ1|yobs
1,y −1
2,...,y −1
k)(7.4)
y∗( )
1∼P(yaus
1|yobs
1,y −1
2,...,y −1
k,θ∗( )
1)(7.5)
... (7.6)
θ∗
k∼P(θk|yobs
k,y
2,...,y
k−1)(7.7)
y∗( )
k∼P(yaus
k|yobs
k,y
2,...,y
k−1,θ∗( )
k)(7.8)
(7.9)
Y si nos ijamos, en ealidad ya podemos ecu i a modelos de impu ación donde
y
es uni a-
ian e pa a mues ea de las an e io es dis ibuciones. En la siguien e sección se de allan a ios
modelos dependiendo de la escala de la a iable, con lo que el p oblema queda ía comple amen e
solucionado.
En ningún caso se u iliza in o mación (cuál se iba a usa , en cualquie caso) de
yaus
j
pa a
mues ea de θ∗( )
j. Las simulaciones del modelo pueden se pa alelizadas ácilmen e.
7.1.1 Modelos de impu ación uni a ian e
Pa a los pasos del mues eo de Gibbs es necesa io ecu i a modelos de impu ación uni a-
ian e. Dependiendo de la dis ibución de y|x, los algo i mos son los siguien es:
Da os no malmen e dis ibuídos (escala)
Pa a es e caso, ecu imos a un modelo de eg esión lineal o dina io. Suponemos media
βx
,
a ianza
σ2
como pa áme os. Po o a pa e descomponemos
x= (xobs,xmis)
como la ma iz
n×k
de da os y
n=nobs +naus
pa a la a iable a p edeci . El algo i mo cons a de los siguien es
pasos:
1.
Es ima
β
po una eg esión o dina ia lineal de
y
sob e
x
. Conc e amen e po :
b=
(xobs0xobs)−1xobs0yobs
2. Mues ea g∼χ2(nobs −p)
3. Es ima σ2∗= (yobs −xobsb)0(yobs −xobsb/g)
4. Mues ea w1∼N(0,Ik)
5.
Calcula
b∗=b+σ2∗w1V1/2
, donde
V1/2
es la ma iz iangula supe io de una decom-
posición Cholesky de V= (xobs0xobs)−1
6.
Mues ea
w2∼N(0,Iaus
n)
, donde es e úl imo é mino ep esen a una ma iz iangula de
amaño ndonde los elemen os donde yies é ausen e alen 1 y 0 en caso con a io.
7. Impu a y∗=xausb∗+w2σ∗
Es e algo i mo debe ía se bas an e obus o a al a de no malidad, pe o no obs an e se p opo-
nen un pa de al e na i as pa a ene en cuen a es e enómeno:
P edic i e mean ma ching
. Bas a sus i ui el úl imo paso calculando
yaus =xausb∗
. Pa a
cada alo ausen e, impu a ese alo po el alo de
yobs =xobsb∗
que más ce cano se encuen e
ayaus.
Ho -deck
Reemplaza el penúl imo paso el mues eo po uno de
naus
alo es con eemplaza-
mien o del conjun o de nobs esiduos es anda izados.
66 Capí ulo 7. Impu ación múl iple
Da os de espues a bina ia
Recu imos a un modelo de eg esión logís ica. Asumimos
P(yi|β,xi) = exiβ
(1+exiβ)yi 1−exiβ
(1+exiβ)1−yi!
.
A con inuación, las impu aciones se ob ienen de la siguien e mane a:
1.
Ob ene po un mé odo numé ico habi ual (New on-Raphson po ejemplo) un es imado
b
de βy de Va [β](un hessiano) usando solo los da os comple os.
2. Mues ea b∗∼N(b,Va [b])
3. Pa a cada obse ación ausen e, calcula wi=exib∗
1+exib∗
4.
Pa a cada obse ación ausen e, mues eamos un
ui∼U(0,1)
. Si
ui>wi
, impu amos esa
obse ación como yi=1, en o o caso, yi=0.
Da os ca egó icos
Es e caso es el más sencillo po que no es más que una gene alización na u al del an e io .
No emos las ca ego ías po
0,...,s−1
. La dis ibución de
y
puede ca ac e iza se po
log(P(y=
j|x)/P(y=0|x)) = βjx
. En es e sen ido, el modelo pa a
y
no es más que una se ie de eg esiones
logís icas apiladas conside ando una ca ego ía base (one s es ).
1.
De mane a simila , u iliza un modelo de eg esión logís ica mul inomial, es ima
b=ˆ
β
y
Va (b) = ˆ
Va (β).
2. Mues ea b∗∼N(b,Va (b))
3. Pa a cada obse ación ausen e, calcula πaus
i,j=e−b∗
jxi
1+∑s−1
=1eb∗
xi
4. Pa a cada obse ación ausen e xi, mues ea de 0,...,s−1 con p obabilidad πaus
i,j.
7.2 P obando el algo i mo
En es a co a sección nos dedica emos a p oba con a ios ma cos de da os de dis in a
na u aleza cómo se compo a el algo i mo a ni el de calidad de impu aciones.
Eje cicio 7.1
En p ime luga p oba emos con el ma co de da os
m ca s
del paque e
da ase s
. Es e ma co se ca ac e iza po ene a iables de odo ipo, po lo que se p es a
muy bien a analiza los e o es come idos. Po o a pa e, es un ma co de da os pequeños, po
lo que se i á pa a e alua la calidad de las impu aciones cuando hay pocas obse aciones
in oluc adas.
lib a y(mice)
da a = m ca s
da a$am = ac o (da a$am)
da a$ s = ac o (da a$ s)
da oso iginales = da a
In oducimos de mane a alea o ia da os ausen es:
n = n ow(da a)
m = ncol(da a)
p = 0.15
7.2 P obando el algo i mo 67
## Espe amos ap ox n*m*p da os ausen es
o (i in 1:n){
o (j in 1:m){
u = uni (1)
i (u < p){
da a[i,j] = NA
}
}
}
expec ed = n*m*p
leng h(which(is.na(da a)))
Co emos el algo i mo, usando pa a odas las a iables P edic i e Mean Ma ching:
impu a ion = mice(da a, m=10, maxi =100, me hod = ep(" as pmm", 11))
comple os = comple e(impu a ion, ac ion = 5)
Finalmen e, e aluamos la calidad de las impu aciones median e la aíz del e o cuad á ico
medio (pa a las obse aciones ca egó icas no obs an e ca ece de sen ido es a medida).
a = as.nume ic(comple os[is.na(da a)])
b = as.nume ic(da oso iginales[is.na(da a)])
mse = sq (mean((a-b)^2))
O o ejemplo ambién con a iables de dis in o ipo muy sencillo es el siguien e:
da a(nhanes2)
head(nhanes2)
impu a ion = mice(nhanes2, m=10, maxi =100, me hod = ep(" as pmm", 4))
comple os = comple e(impu a ion, ac ion = 5)
p in (comple os)
Descub imien o de emá icas
Asignación la en e de Di ichle
Fo malización de la gene ación
In e encia y es imación usando mues eo
de Gibbs colapsado
Fo malización del ap endizaje
Ejemplos de aplicación
Análisis de pe il en Twi e
Temá ica según au o es clásicos
8. Mine ía de ex o
8.1 Descub imien o de emá icas
En es e capí ulo abo da emos un p oblema ela i amen e ecien e de la In eligencia A i i-
cial. Supongamos que enemos una g an can idad de documen os que analiza , con el mismo o
di e en es p opósi os. Una de las p oblemá icas que nos puede su gi en e an a in o mación es
la de busca emá icas comunes en e di e sos documen os. Es o puede ene mucha u ilidad en
di e sos campos, pe o especialmen e en mo o es de búsqueda, donde las páginas (especialmen e
de no icias), quedan indexadas según una se ie de palab as cla e. És as se ían pos e io men e
ag upadas au omá icamen e po el mismo mo o pa a p opo ciona in o mación de una mane a
o denada.
En de ini i a, nues o p oblema explicado de mane a in o mal es el siguien e: dados una
se ie de documen os de ex o cualesquie a, de e mina si exis en
k
emá icas ela i as a esos
documen os. Es e núme o
k
de emá icas iene ijado de an emano. En los úl imos años se han
desa ollado mul i ud de écnicas pa a in en a sol en a es e p oblema. Pa a es e abajo amos a
cen a nos en una llamada La en Di ichle Alloca ion (que aduci emos po Asignación la en e
de Di ichle ), cuya es imación de pa áme os puede se solucionada po écnicas MCMC. Nó ese
que es as écnicas no son las únicas pa a esol e el p oblema, pe o como e emos se p es an a
se las más áciles de implemen a en es e caso.
La asignación la en e de Di ichle (la llama emos LDA po comodidad a pa i de aho a)
es á implemen ada en a ios lenguajes de p og amación, siendo la más ácil de u iliza una de
Py hon. És a, además implemen a dicha écnica usando el mues eo de Gibbs (conc e amen e
una e sión llamada mues eo de Gibbs colapsado), po lo que esul a especialmen e adecuada
pa a ejempli ica .
8.2 Asignación la en e de Di ichle
An es de comenza a explica de mane a o mal en qué consis e dicha écnica, amos a
e cómo unciona de mane a cons uc i a: LDA supone que los documen os son mix u as de
emá icas, y que las emá icas a su ez son mix u as de palab as. Es deci , supongamos en p ime
luga que que emos gene a un documen o con espec o a es as suposiciones:
70 Capí ulo 8. Mine ía de ex o
1.
Decidi íamos en p ime luga cuán as palab as
N
a a ene un documen o, po ejemplo de
acue do con una dis ibución de Poisson.
2.
Decidimos una dis ibución pa a las emá icas, po ejemplo de acue do con una dis ibución
Di ichle .
3. Aho a gene amos cada palab a del documen o de la mane a siguien e:
a)
Decidimos la emá ica de la palab a, escogiendo una de acue do a una dis ibución
mul inomial, conjugada con la Di ichle an e io .
b)
Una ez escogida la emá ica, escogemos la palab a de acue do a la dis ibución
mul inomial de la emá ica.
8.2.1 Fo malización de la gene ación
De inamos lo siguien e:
Una palab a es la unidad básica de da o disc e o, de inida como un miemb o de ocabula io
indexado en e
1,...,V
. Rep esen amos las palab as como ec o es uni a ios, que ienen un
solo componen e igual a 1 y el es o igual a 0. La
-ésima palab a del ocabula io queda
ep esen ada po w =1 solo pa a .
Un documen o es una secuencia de Npalab as deno adas po w= (w1,...,wN).
Un cue po es una colección de Mdocumen os deno ados po D=w1,...,wM
De mane a o mal, pa a gene a un documen o lo que ha íamos se ía lo siguien e:
1. Escoge N∼Poisson(ε)
2. Escogemos θ∼Di (α)
3. Pa a cada una de las Npalab as wi:
a) Escogemos una emá ica según zn∼Mul inomial(θ)
b) Escogemos una palab a wicon p obabilidad P(wi|zi,β).
Realiza emos algunas hipó esis de simpli icación pa a con inua . Pa a empeza y como ya
hemos comen ado, la dimensionalidad
k
de la dis ibución Di ichle queda ijada de an emano
(y po an o el núme o de emá icas a encon a ). En segundo luga , las p obabilidades de cada
una de las palab as quedan pa ame izadas po una ma iz
β
de dimensiones
k×V
, donde
βi j =P(wj=1|zj=1)
. En un p incipio a a emos es a ma iz como algo ijado que end emos
que es ima e en ualmen e. Po úl imo, la suposición de que los documen os ienen un núme o
de palab as de acue do con una dis ibución de Poisson puede se igno ada comple amen e.
8.3 In e encia y es imación usando mues eo de Gibbs colapsado
Has a aho a hemos ap endido cómo gene a documen os de acue do con es e modelo. He-
mos empezado así po que LDA es un modelo gene a i o, en el sen ido de que los documen os
ienen gene ados según esa se ie de hipó esis pa a luego ealiza el ap endizaje de mane a más
sencilla. En p ime luga explica emos de mane a simila a como hemos hecho an es el p oce-
dimien o de ap endizaje de mane a in o mal, y luego lo o maliza emos de mane a más de allada.
Nues o caso aho a se cen a en que enemos una se ie de documen os a los cuales que emos
ex ae emá icas:
De mane a alea o ia, pa a cada cada palab a en cada documen o, asignamos una emá ica.
Es a asignación ya nos da una ep esen ación de la dis ibución de las emá icas y de las
palab as (si bien nada buenas). Pa a mejo a las:
•
Pa a cada palab a en cada documen o y pa a cada emá ica calcula 1)
P( em ica |documen o d) =
la p opo ción de palab as en el documen o
d
que es án asignadas a la emá ica
y
2)
P(palab a w| em ica ) =
la p opo ción de asignaciones a la emá ica
debidas a
es a palab a w.
8.3 In e encia y es imación usando mues eo de Gibbs colapsado 71
•
Reasigna a la palab a una nue a emá ica de acue do con
P( em ica |documen o d)×
P(palab a w| em ica )pa a odas las emá icas
Después de epe i es e p ocedimien o un núme o ele ado de eces, alcanza emos un es ado
es able en el que las asignaciones apenas cambian. En es e es ado ya podemos usa las asignacio-
nes pa a es ima las mix u as en cada documen o (simplemen e con ando las palab as asignadas
a cada emá ica en cada documen o) y las palab as asociadas a cada emá ica (con ando el o al
de palab as asignadas a cada emá ica).
8.3.1 Fo malización del ap endizaje
Dados
D
documen os con eniendo
T
emá icas exp esadas bajo
W
palab as únicas, podemos
ep esen a
P(w|z)
como un conjun o de
T
dis ibuciones mul inomiales
φ
sob e las
W
palab as,
de al mane a que
P(w|z=j) = φ(j)
w
y equi alen emen e, ep esen amos
P(z)
como
D
dis ibu-
ciones mul inomiales sob e las
T
emá icas, de al mane a que pa a un documen o cualquie a
P(z=j) = θ(d)
j.
Nues a es a egia pa a descub i las emá icas pasa po explo a la dis ibución a pos e io i de
las emá icas sob e las palab as
P(z|w)
. Una ez e aluada, pod emos ob ene es imaciones an o
de
φ
como de
θ
. (Es as dos úl imas pueden se ep esen adas como ma ices po con eniencia).
U ilizamos el mues eo de Gibbs pa a mues ea de es a dis ibución a pos e io i, pa a ello,
de inimos lo siguien e:
wi|zi,φ(zi)∼Mul inomial(φ(zi))
φ∼Di ichle (β)
zi|θ(di)∼Mul inomial(θ(di))
θ∼Di ichle (α)
Donde
α
y
β
son hipe pa áme os pa a con ola las dis ibuciones a p io i de
φ
y
θ
. Es as
dis ibuciones a p io i son conjugadas con la Mul inomial, pe mi iendo calcula la dis ibución
conjun a
P(w,z) = P(w|z)P(z)
. Si en es a dis ibución in eg amos p ime o sob e
φ
ob end emos:
P(w|z) = Γ(Wβ)
Γ(β)WTT
∏
j=1
∏wΓ(n(w)
j+β)
Γ(n(.)
j+Wβ)(8.1)
donde
n(w)
j
es el núme o de eces que una palab a
w
ha sido asignada a la emá ica
j
en el
ec o de asignaciones z. Si en el segundo é mnino in eg amos sob e θ, ob end íamos:
P(z) = Γ(Tα)
Γ(α)TDD
∏
d=1
∏jΓ(n(d)
j+α)
Γ(n(d)
.+Tα)(8.2)
donde
n(d)
j
es el núme o de eces que una de e minada palab a de un documen o
d
ha sido
asignada al documen o j. Po an o, nues o obje i o se ía e alua la dis ibución a pos e io i:
P(z|w) = P(w,z)
∑zP(w,z)(8.3)
Desg aciadamen e, es a dis ibución no puede se calculada di ec amen e, ya que su cálculo
implica la e aluación de una dis ibución de p obabilidad sob e un espacio de es ados disc e os
muy amplios. No obs an e, podemos usa el mues eo de Gibss, pa a ob ene la exp esión,
72 Capí ulo 8. Mine ía de ex o
median e cancelación de é minos de las dos dis ibuciones an e io es, ob eniendo la siguien e
dis ibución o almen e condicionada:
P(zi=j|z−i,w)∝n(wi)
−i,j+β
n(.)
−i,j+Wβ
n(di)
−i,j+α
n(di)
−i+Tα
(8.4)
donde
n.
−i
es un con eo sin con a la asignación ac ual de
zi
. Es e esul ado es bas an e
in ui i o, ya que de la exp esión an e io el p ime é mino exp esa la p obabilidad de
wi
bajo
la emá ica
j
, y el segundo a io exp esa la p obabilidad de la emá ica
j
sob e el documen o
di
. Una ez ob enida es a dis ibución o almen e condicionada, el algo i mo p ocede como
sigue: Inicializamos odas las palab as de los documen os ( ealmen e sólo aquellas incluídas
en e las
W
del ocabula io) con alo es en e
1
y
T
. A con inuación asignamos emá icas a
las palab as de acue do a las p obabilidades ob enidas con la exp esión an e io . Lo más ade-
cuado, es que en cada i e ación dichas p obabilidades sean halladas únicamen e con aquellas
palab as cuya asignación haya sido ealizado an e io men e. Una ez que la cadena ha co ido un
núme o su icien emen e al o de eces, llega á a un es ado su icien emen e es able (o es aciona io).
Con un conjun o de mues as de la dis ibución a pos e io i ya end íamos esuel o el
p oblema de asignación. No obs an e, si es u ié amos in e esados en conoce la o ma de los
pa áme os dis ibucionales pod íamos es ima los median e dicha mues a de la siguien e mane a:
ˆ
φ(w)
j=n(w)
j+β
n(.)
j+Wβ
(8.5)
ˆ
θ(d)
j=n(d)
j+α
n(d)
.+Tα(8.6)
Escogiendo los pa áme os α,βy el núme o de emá icas T
El algo i mo desc i o en es a sección pod ía ex ende se al caso en el que
α
y
β
ue an a su
ez dis ibuciones dependien es de más hipe pa áme os y mues ea de ellos, pe o ealiza lo
ealmen e con ibuye a un mayo cos e compu acional sin ob ene mejo as de endimien o.
La elección de es os dos pa áme os es impo an e a la ho a de ob ene dis in os esul ados:
en pa icula , un
β
más g ande p oduce un menos núme o de emá icas más amplias. Dados
alo es ijos a es os hipe pa áme os se a a ía de escoge un núme o de emá icas ap opiado,
lo que ealmen e es un p oblema de selección de modelos bayesiano. Lo na u al en es adís ica
bayesiana es e alua la dis ibución a pos e io i de los modelos dados los da os, y escoge aquel
que de mayo densidad de p obabilidad. En nues o caso, nues os da os son las palab as en los
documen os, luego necesi a íamos e alua la e osimili ud P(w|T).
De nue o, es e es un p oblema in a able desde el pun o de is a compu acional, pe o
pod íamos ap oxima dicha dis ibución omando la media ha mónica de un conjun o de alo es
de
P(w|z,T)
cuando
z
se mues ea de la dis ibución a pos e io i
P(z|w,T)
. Es as úl imas pueden
halla se usando 8.1. Una ez ap oximada, podemos u iliza es a can idad pa a escoge el núme o
de emá icas adecuadamen e.
8.4 Ejemplos de aplicación
En es a sección nos dedica emos a aplica es e algo i mo en algunas si uaciones in e esan es
en las que enga sen ido. Desde el pun o de is a p ác ico es a emos u ilizando el paque e
lda
disponible en
PyPI
, an o bajo
Py hon2
como
Py hon3
. En es e caso u iliza emos es a úl ima
implemen ación po e i a p oblemas con la au en icación ia
oau h2
.
8.4 Ejemplos de aplicación 73
8.4.1 Análisis de pe il en Twi e
Con el auge de las nue as edes sociales, un nue o mundo de in o mación ex ual se ab e.
Es a in o mación es ácilmen e accesible pa a cualquie a con una cuen a de po ejemplo Twi e
(con algunas limi aciones) desde cualquie lenguaje de p og amación mode no. En nues o caso,
lo que amos a ealiza con el siguien e código es un análisis de pe il de una de e minada cuen a
de es a ed social. Conc e amen e, es a íamos in e esados en conoce el con enido en cuan o a
emá icas de los wee s de una de e minada pe sona.
Cen émonos en analiza de qué habla Elon Musk, CEO de Tesla Mo o s. Pa a ello necesi a-
mos en p ime luga una se ie de cla es de acceso a la API de Twi e , las cuales son ob enibles
ácilmen e a a és de la página web. Comencemos el análisis.
Eje cicio 8.1
#!/us /bin/en py hon
# -*- coding: u -8 -*-
"""
La en Di ichle Alloca ion
Twi e C awle
@au ho : JJimenez
"""
impo numpy as np
access_ oken = "xxxxx-xxxxxx"
access_ oken_sec e = "xxxxx"
consume _key = "xxxxxxxxxxxxxxxxxxxxxxxxx"
consume _sec e = "xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx"
impo weepy
au h = weepy.OAu hHandle (consume _key, consume _sec e )
au h.se _access_ oken(access_ oken, access_ oken_sec e )
api = weepy.API(au h)
public_ wee s = api.use _ imeline('elonmusk', coun = 1000)
ex s = []
o wee in public_ wee s:
ex s.append( wee . ex )
om sklea n impo ea u e_ex ac ion
pa a ocab = ea u e_ex ac ion. ex .Coun Vec o ize (s op_wo ds = 'english',
max_ ea u es = 50). i ( ex s)
ocab = pa a ocab. ocabula y_
ocab = ocab.keys()
80 BIBLIOGRAFÍA
[MR14]
Jean-Michel Ma in y Ch is ian P Robe . Bayesian Essen ials wi h R. 2014, pági-
na 296. ISBN: 978-1-4614-8686-2. DOI:
10.1007/978-1-4614-8687-9
.URL:
h p://link.sp inge .com/10.1007/978-1-4614-8687-9
.
[Mcm05]
Da a-d i en Mcmc. “MCMC Tu o ial a ICCV Mo i a ion o MCMC”. En: Oc o-
be (2005).
[MJG12] Mo phome ic Mcmc, Lei T Johnson y Cha les J Geye . “Mo phome ic MCMC
(mcmc Package Ve . 0.9)”. En: (2012), páginas 1-20.
[Mcm+05]
T ans-dimensional Mcmc y col. “Re e ences
•
Pe e G een , Re e sible jump
Ma ko chain Mon e Ca lo , in “ Highly S uc u ed S ochas ic Sys ems ”, 2003
•
Paskin & Th un , Robo ic Mapping wi h Polygonal Ou line Model Selec ion
•
Mos common case : in e ence on p ( x | z ), x con inuous p”. En: Oc obe Oc obe
(2005), páginas 1-12.
[Ml]
An Ml. “Mix u e models ( Ch . 16 ) Mix u e models ( con ’ d ) Bayesian es ima ion
( con ’ d ) Mix u e models - missing da a”. En: 1 ().
[Ric+14] Au ho Richa d y col. “Package ‘ BayesFac o ’”. En: (2014).
[Sah00]
Suji Sahu. “Ma ko Chain Mon e Ca lo ( MCMC ) In oduc ion”. En: Augus
(2000).
[VA13]
Sunay Vaishna y P ima y Ad iso . “A Ma ko Chain Mon e Ca lo based app oach
o Image Segmen a ion”. En: (2013).
[WG00]
Michael D Wa d y K is ian Sk ede Gledi sch. “Loca ion , Loca ion , Loca ion : An
MCMC App oach o Modeling Spa ial Con ex wi h Ca ego ical Va iables in he
S udy and P edic ion o Wa 1”. En: (2000).
[Zhu+05]
Song-chun Zhu y col. “Ma ko Chain Mon e Ca lo o Compu e Vision (SLIDES)”.
En: Oc obe Oc obe (2005).
[ZS02]
Zhuowen Tu y Song-Chun Zhu. “Image segmen a ion by da a-d i en ma ko chain
mon e ca lo”. En: IEEE T ansac ions on Pa e n Analysis and Machine In elligence
24.5 (2002), páginas 657-673. ISSN: 01628828. DOI:
10.1109/34.1000239
.