Sampling importance
resampling (SIR)
Em geral, o problema computacional da inferência Bayesiana é
conseguir amostrar da posteriori \(p(\theta|dados)\), que é obtida pela
combinação de uma verossimilhanca \(L(\theta|dados)\) com uma priori \(p(\theta)\), i.e. \[
p(\theta|dados) = \frac1kp(\theta)L(\theta|dados).
\] Como sabemos avaliar \(p(\theta)\) e \(L(\theta|dados)\) pontualmente, sabemos
também avaliar \(p(\theta|dados)\)
pontualmente, a menos de uma constante (\(1/k\) aqui!).
O que o algoritmo SIR faz é amostrar de uma distribuição auxiliar (ou
proposta ou candidata) \(q(\theta)\)
que i) seja fácil de amsotrar; ii) possa ser avaliada pontualmente e
iii) se “aproxime” da posterior \(p(\theta|dados)\).
Daí o SIR funciona assim:
Uma amostra \(\{\theta^{(1)},\ldots,\theta^{(M)}\}\) é
obtida da densidade proposta \(q(\theta)\). O tamanho da amostra Monte
Carlo pode ser \(M=1.000\) ou \(M=10.000\) ou ainda \(M=100.000\), ou qualquer outra quantidade
dependendo da complexidade contextual.
Cada componente \(\theta^{(i)}\)
da amostra recebe um peso \[
w^{(i)} = \frac{p(\theta^{(i)}|dados)}{q(\theta^{(i)})},
\] ou seja, \[
w^{(i)} = \frac1k
\frac{p(\theta^{(i)}L(\theta^{(i)}|dados)}{q(\theta^{(i)})}.
\] A constante \(k\) é
irrelevante pois os pesos podem ser recalculados dividindo-os pela sua
soma; daí o \(k\) é cancelado no
numerador e no denominador.
Reamostra-se \(N \leq M\)
elementos do conjunto \(\{\theta^{(1)},\ldots,\theta^{(M)}\}\) com
pesos \(\{w^{(1)},\ldots,w^{(M)}\}\) e
com reposição. Assim, obtem-se uma amostra da posterior \(p(\theta|dados)\):
\[
\left\{\theta_*^{(1)},\ldots,\theta_*^{(N)}\right\}.
\]
Caso particular:
Utilização da priori como proposta
Nesse caso \(q(\theta) = p(\theta)\)
e os pesos ficam \[
w^{(i)} = \frac1k
\frac{p(\theta^{(i)})L(\theta^{(i)}|dados)}{p(\theta^{(i)})} =
\frac1k L(\theta^{(i)}|dados).
\] Portanto, os pesos do SIR só serão a verossimilhança se a
densidade proposta for a própria priori.
Observação: Apesar da convenienca de se amostrar da
priori e calcular pesos só com a verossimilhanca, o grande problema
reside no fato que em geral priori e verossimilhança podem estar em
regiões distintas e portanto amostrar da priori pode não ser uma boa
estratégia. Lembre-se que eu disse acima que desejamos uma proposta que
se “aproxime” da posterior \(p(\theta|dados)\).
Exemplos
Proposta envelopa a
distribuição de interesse
Distribuição de interesse: \(h(\theta) =
\frac1k exp(-0.5\theta^2)\)
Aqui \(k\) é a constante
normalizadora.
Distribuição proposta: \(N(0,5^2)\)
M = 10000
N = 10000
theta = rnorm(M,0,5)
w = exp(-0.5*theta^2)/dnorm(theta,0,5)
w = w/sum(w)
thetastar = sample(theta,size=N,replace=TRUE,prob=w)
par(mfrow=c(1,2))
thetas = seq(-5,5,length=1000)
hist(theta,prob=TRUE,main="Amostras de q(theta)",ylim=c(0,0.4))
lines(thetas,dnorm(thetas),lwd=2,col=2)
hist(thetastar,xlab="theta",prob=TRUE,main="Amostras de h(theta)")
lines(thetas,dnorm(thetas),lwd=2,col=2)
Proposta NÃO envelopa
a distribuição de interesse
Distribuição de interesse: \(h(\theta) =
\frac1k exp(-0.5\theta^2)\)
Novamente, \(k\) é a constante
normalizadora.
Distribuição proposta: \(N(-15,5^2)\)
M = 10000
N = 10000
theta = rnorm(M,-15,5)
w = exp(-0.5*theta^2)/dnorm(theta,-15,5)
w = w/sum(w)
thetastar = sample(theta,size=N,replace=TRUE,prob=w)
par(mfrow=c(1,2))
thetas = seq(-5,5,length=1000)
hist(theta,prob=TRUE,main="Amostras de q(theta)",ylim=c(0,0.4))
lines(thetas,dnorm(thetas),lwd=2,col=2)
hist(thetastar,xlab="theta",prob=TRUE,main="Amostras de h(theta)")
lines(thetas,dnorm(thetas),lwd=2,col=2)
SIR na inferência
Bayesiana
Dados = \(\{1,0,1,1,0,1,0,0,1,1\}\)
Modelo: \(x_i|\theta \ iid \
Bernoulli(\theta)\) for \(i=1,\ldots,n\) and \(\theta \in (0,1)\).
Priori: \(\theta \sim Beta(6,2)\),
such that \(E(\theta)=0.75\) and \(V(\theta)=0.144^2\).
Posteriori (distribuição de interesse) \[
p(\theta|dados) = \frac1k \theta^{6-1}(1-\theta)^{2-1}
\theta^6(1-\theta)^4
\] Novamente, \(k\) é a
constante de normalização.
Distribuição proposta: \(Beta(2,2)\)
Você pode trocar essa proposta por qualquer distribuição no intervalo
\((0,1)\).
M = 10000
N = 10000
theta = rbeta(M,2,2)
w = dbeta(theta,6,2)*theta^6*(1-theta)^4/dbeta(theta,2,2)
w = w/sum(w)
thetastar = sample(theta,size=N,replace=TRUE,prob=w)
par(mfrow=c(1,2))
thetas = seq(0,1,length=1000)
hist(theta,prob=TRUE,main="Amostras da proposta",ylim=c(0,4),xlim=c(0,1))
lines(thetas,dbeta(thetas,12,6),lwd=2,col=2)
hist(thetastar,xlab="theta",prob=TRUE,main="Amostras da posteriori",ylim=c(0,4),xlim=c(0,1))
lines(thetas,dbeta(thetas,12,6),lwd=2,col=2)
c(mean(thetastar), sqrt(var(thetastar)))
## [1] 0.6659571 0.1076641
quantile(thetastar,c(0.025,0.5,0.975))
## 2.5% 50% 97.5%
## 0.4426397 0.6704314 0.8562359
Vamos repetir 2.3)
mas usando a priori como proposta
Distribuição proposta = Priori: \(\theta
\sim Beta(6,2)\).
M = 10000
N = 10000
theta = rbeta(M,6,2)
w = theta^6*(1-theta)^4
w = w/sum(w)
thetastar = sample(theta,size=N,replace=TRUE,prob=w)
par(mfrow=c(1,2))
thetas = seq(0,1,length=1000)
hist(theta,prob=TRUE,main="Amostras da proposta",ylim=c(0,4),xlim=c(0,1))
lines(thetas,dbeta(thetas,12,6),lwd=2,col=2)
hist(thetastar,xlab="theta",prob=TRUE,main="Amostras da posteriori",,ylim=c(0,4),xlim=c(0,1))
lines(thetas,dbeta(thetas,12,6),lwd=2,col=2)
c(mean(thetastar), sqrt(var(thetastar)))
## [1] 0.667744 0.107167
quantile(thetastar,c(0.025,0.5,0.975))
## 2.5% 50% 97.5%
## 0.4460907 0.6722582 0.8584376
Vamos repetir 2.4)
mas agora com MUITO mais dados e priori MUITO longe da
verossimilhança
Dados = \(\{1,0,1,1,0,1,0,0,1,1,1,0,1,1,0,1,0,0,1,1,1,0,1,1,0,1,0,0,1,1\}\)
Priori: \(\theta \sim
Beta(3,12)\)
Verossimilhança: muito mais concentrada! \[
L(\theta|dados) = \theta^{18}(1-\theta)^{12}
\]
M = 10000
N = 10000
theta = rbeta(M,3,12)
w = theta^(18)*(1-theta)^(12)
w = w/sum(w)
thetastar = sample(theta,size=N,replace=TRUE,prob=w)
par(mfrow=c(1,2))
thetas = seq(0,1,length=1000)
hist(theta,prob=TRUE,main="Amostras da proposta",ylim=c(0,6),xlim=c(0,1))
lines(thetas,dbeta(thetas,21,24),lwd=2,col=2)
hist(thetastar,xlab="theta",prob=TRUE,main="Amostras da posteriori",ylim=c(0,6),xlim=c(0,1))
lines(thetas,dbeta(thetas,21,24),lwd=2,col=2)
c(mean(thetastar), sqrt(var(thetastar)))
## [1] 0.45943438 0.07241472
quantile(thetastar,c(0.025,0.5,0.975))
## 2.5% 50% 97.5%
## 0.3217701 0.4584958 0.6163900
LS0tCnRpdGxlOiAiU2FtcGxpbmcgSW1wb3J0YW5jZSBSZXNhbXBsaW5nIgpzdWJ0aXRsZTogIk1haXMgdW5zIGV4ZW1wbG9zIgphdXRob3I6ICJIZWRpYmVydCBGcmVpdGFzIExvcGVzIgpkYXRlOiAiNS8xMy8yMDI0IgpvdXRwdXQ6CiAgaHRtbF9kb2N1bWVudDoKICAgIHRvYzogdHJ1ZQogICAgdG9jX2RlcHRoOiAzCiAgICB0b2NfZmxvYXQ6IHRydWUKICAgIHRvY19jb2xsYXBzZWQ6IHRydWUKICAgIGNvZGVfZG93bmxvYWQ6IHllcwogICAgbnVtYmVyX3NlY3Rpb25zOiB0cnVlCi0tLQoKIyBTYW1wbGluZyBpbXBvcnRhbmNlIHJlc2FtcGxpbmcgKFNJUikKCkVtIGdlcmFsLCBvIHByb2JsZW1hIGNvbXB1dGFjaW9uYWwgZGEgaW5mZXLDqm5jaWEgQmF5ZXNpYW5hIArDqSBjb25zZWd1aXIgKmFtb3N0cmFyKiBkYSBwb3N0ZXJpb3JpICRwKFx0aGV0YXxkYWRvcykkLCBxdWUgw6kgCm9idGlkYSBwZWxhIGNvbWJpbmHDp8OjbyBkZSB1bWEgdmVyb3NzaW1pbGhhbmNhICRMKFx0aGV0YXxkYWRvcykkIApjb20gdW1hIHByaW9yaSAkcChcdGhldGEpJCwgaS5lLgokJApwKFx0aGV0YXxkYWRvcykgPSBcZnJhYzFrcChcdGhldGEpTChcdGhldGF8ZGFkb3MpLgokJApDb21vIHNhYmVtb3MgYXZhbGlhciAkcChcdGhldGEpJCBlICRMKFx0aGV0YXxkYWRvcykkIHBvbnR1YWxtZW50ZSwgCnNhYmVtb3MgdGFtYsOpbSBhdmFsaWFyICRwKFx0aGV0YXxkYWRvcykkIHBvbnR1YWxtZW50ZSwgYSBtZW5vcwpkZSB1bWEgY29uc3RhbnRlICgkMS9rJCBhcXVpISkuCgpPIHF1ZSBvIGFsZ29yaXRtbyBTSVIgZmF6IMOpIGFtb3N0cmFyIGRlIHVtYSBkaXN0cmlidWnDp8OjbyAKYXV4aWxpYXIgKG91IHByb3Bvc3RhIG91IGNhbmRpZGF0YSkgJHEoXHRoZXRhKSQgcXVlIGkpIHNlamEgCmbDoWNpbCBkZSBhbXNvdHJhcjsgaWkpIHBvc3NhIHNlciBhdmFsaWFkYSBwb250dWFsbWVudGUgZSBpaWkpIApzZSAiYXByb3hpbWUiIGRhIHBvc3RlcmlvciAkcChcdGhldGF8ZGFkb3MpJC4gIAoKRGHDrSBvIFNJUiBmdW5jaW9uYSBhc3NpbToKCjEpIFVtYSBhbW9zdHJhICRce1x0aGV0YV57KDEpfSxcbGRvdHMsXHRoZXRhXnsoTSl9XH0kIMOpIG9idGlkYSBkYSBkZW5zaWRhZGUgcHJvcG9zdGEgJHEoXHRoZXRhKSQuCk8gdGFtYW5obyBkYSBhbW9zdHJhIE1vbnRlIENhcmxvIHBvZGUgc2VyICRNPTEuMDAwJCBvdSAkTT0xMC4wMDAkIG91IGFpbmRhICRNPTEwMC4wMDAkLCBvdSAKcXVhbHF1ZXIgb3V0cmEgcXVhbnRpZGFkZSBkZXBlbmRlbmRvIGRhIGNvbXBsZXhpZGFkZSBjb250ZXh0dWFsLgoKMikgQ2FkYSBjb21wb25lbnRlICRcdGhldGFeeyhpKX0kIGRhIGFtb3N0cmEgcmVjZWJlIHVtIHBlc28KJCQKd157KGkpfSA9ICBcZnJhY3twKFx0aGV0YV57KGkpfXxkYWRvcyl9e3EoXHRoZXRhXnsoaSl9KX0sCiQkCm91IHNlamEsCiQkCndeeyhpKX0gPSBcZnJhYzFrIFxmcmFje3AoXHRoZXRhXnsoaSl9TChcdGhldGFeeyhpKX18ZGFkb3MpfXtxKFx0aGV0YV57KGkpfSl9LgokJApBIGNvbnN0YW50ZSAkayQgw6kgaXJyZWxldmFudGUgcG9pcyBvcyBwZXNvcyBwb2RlbSBzZXIgcmVjYWxjdWxhZG9zIGRpdmlkaW5kby1vcyBwZWxhIHN1YSBzb21hOyAKZGHDrSBvICRrJCDDqSBjYW5jZWxhZG8gbm8gbnVtZXJhZG9yIGUgbm8gZGVub21pbmFkb3IuCgozKSBSZWFtb3N0cmEtc2UgJE4gXGxlcSBNJCBlbGVtZW50b3MgZG8gY29uanVudG8gJFx7XHRoZXRhXnsoMSl9LFxsZG90cyxcdGhldGFeeyhNKX1cfSQgY29tIApwZXNvcyAkXHt3XnsoMSl9LFxsZG90cyx3XnsoTSl9XH0kIGUgY29tIHJlcG9zacOnw6NvLiAgQXNzaW0sIG9idGVtLXNlIHVtYSBhbW9zdHJhIGRhIHBvc3RlcmlvciAKJHAoXHRoZXRhfGRhZG9zKSQ6CgokJApcbGVmdFx7XHRoZXRhXypeeygxKX0sXGxkb3RzLFx0aGV0YV8qXnsoTil9XHJpZ2h0XH0uCiQkICAgCgojIyBDYXNvIHBhcnRpY3VsYXI6ICBVdGlsaXphw6fDo28gZGEgcHJpb3JpIGNvbW8gcHJvcG9zdGEKTmVzc2UgY2FzbyAkcShcdGhldGEpID0gcChcdGhldGEpJCBlIG9zIHBlc29zIGZpY2FtCiQkCndeeyhpKX0gPSBcZnJhYzFrIFxmcmFje3AoXHRoZXRhXnsoaSl9KUwoXHRoZXRhXnsoaSl9fGRhZG9zKX17cChcdGhldGFeeyhpKX0pfSA9IApcZnJhYzFrIEwoXHRoZXRhXnsoaSl9fGRhZG9zKS4KJCQKUG9ydGFudG8sIG9zIHBlc29zIGRvIFNJUiBzw7Mgc2Vyw6NvIGEgdmVyb3NzaW1pbGhhbsOnYSBzZSBhIGRlbnNpZGFkZSBwcm9wb3N0YSBmb3IgYSBwcsOzcHJpYSBwcmlvcmkuCgoqKk9ic2VydmHDp8OjbzoqKiAgQXBlc2FyIGRhIGNvbnZlbmllbmNhIGRlIHNlIGFtb3N0cmFyIGRhIHByaW9yaSAKZSBjYWxjdWxhciBwZXNvcyBzw7MgY29tIGEgdmVyb3NzaW1pbGhhbmNhLCBvIGdyYW5kZSBwcm9ibGVtYSAKcmVzaWRlIG5vIGZhdG8gcXVlIGVtIGdlcmFsIHByaW9yaSBlIHZlcm9zc2ltaWxoYW7Dp2EgcG9kZW0gCmVzdGFyIGVtIHJlZ2nDtWVzIGRpc3RpbnRhcyBlIHBvcnRhbnRvIGFtb3N0cmFyIGRhIHByaW9yaSBwb2RlIApuw6NvIHNlciB1bWEgYm9hIGVzdHJhdMOpZ2lhLiAgTGVtYnJlLXNlIHF1ZSBldSBkaXNzZSBhY2ltYSBxdWUgCmRlc2VqYW1vcyB1bWEgcHJvcG9zdGEgcXVlIHNlICJhcHJveGltZSIgZGEgcG9zdGVyaW9yIAokcChcdGhldGF8ZGFkb3MpJC4KCgojIEV4ZW1wbG9zIAoKCiMjIFByb3Bvc3RhIGVudmVsb3BhIGEgZGlzdHJpYnVpw6fDo28gZGUgaW50ZXJlc3NlCgpEaXN0cmlidWnDp8OjbyBkZSBpbnRlcmVzc2U6ICAkaChcdGhldGEpID0gXGZyYWMxayBleHAoLTAuNVx0aGV0YV4yKSQKCkFxdWkgJGskIMOpIGEgY29uc3RhbnRlIG5vcm1hbGl6YWRvcmEuCgpEaXN0cmlidWnDp8OjbyBwcm9wb3N0YTogJE4oMCw1XjIpJAoKYGBge3J9Ck0gICAgID0gMTAwMDAKTiAgICAgPSAxMDAwMAp0aGV0YSA9IHJub3JtKE0sMCw1KQp3ICAgICA9IGV4cCgtMC41KnRoZXRhXjIpL2Rub3JtKHRoZXRhLDAsNSkKdyAgICAgPSB3L3N1bSh3KSAKdGhldGFzdGFyID0gc2FtcGxlKHRoZXRhLHNpemU9TixyZXBsYWNlPVRSVUUscHJvYj13KQoKcGFyKG1mcm93PWMoMSwyKSkKdGhldGFzID0gc2VxKC01LDUsbGVuZ3RoPTEwMDApCmhpc3QodGhldGEscHJvYj1UUlVFLG1haW49IkFtb3N0cmFzIGRlIHEodGhldGEpIix5bGltPWMoMCwwLjQpKQpsaW5lcyh0aGV0YXMsZG5vcm0odGhldGFzKSxsd2Q9Mixjb2w9MikKaGlzdCh0aGV0YXN0YXIseGxhYj0idGhldGEiLHByb2I9VFJVRSxtYWluPSJBbW9zdHJhcyBkZSBoKHRoZXRhKSIpCmxpbmVzKHRoZXRhcyxkbm9ybSh0aGV0YXMpLGx3ZD0yLGNvbD0yKQpgYGAKCgojIyBQcm9wb3N0YSBOw4NPIGVudmVsb3BhIGEgZGlzdHJpYnVpw6fDo28gZGUgaW50ZXJlc3NlCgpEaXN0cmlidWnDp8OjbyBkZSBpbnRlcmVzc2U6ICAkaChcdGhldGEpID0gXGZyYWMxayBleHAoLTAuNVx0aGV0YV4yKSQKCk5vdmFtZW50ZSwgJGskIMOpIGEgY29uc3RhbnRlIG5vcm1hbGl6YWRvcmEuCgpEaXN0cmlidWnDp8OjbyBwcm9wb3N0YTogJE4oLTE1LDVeMikkCgpgYGB7cn0KTSAgICAgPSAxMDAwMApOICAgICA9IDEwMDAwCnRoZXRhID0gcm5vcm0oTSwtMTUsNSkKdyA9IGV4cCgtMC41KnRoZXRhXjIpL2Rub3JtKHRoZXRhLC0xNSw1KQp3ID0gdy9zdW0odykgCnRoZXRhc3RhciA9IHNhbXBsZSh0aGV0YSxzaXplPU4scmVwbGFjZT1UUlVFLHByb2I9dykKCnBhcihtZnJvdz1jKDEsMikpCnRoZXRhcyA9IHNlcSgtNSw1LGxlbmd0aD0xMDAwKQpoaXN0KHRoZXRhLHByb2I9VFJVRSxtYWluPSJBbW9zdHJhcyBkZSBxKHRoZXRhKSIseWxpbT1jKDAsMC40KSkKbGluZXModGhldGFzLGRub3JtKHRoZXRhcyksbHdkPTIsY29sPTIpCmhpc3QodGhldGFzdGFyLHhsYWI9InRoZXRhIixwcm9iPVRSVUUsbWFpbj0iQW1vc3RyYXMgZGUgaCh0aGV0YSkiKQpsaW5lcyh0aGV0YXMsZG5vcm0odGhldGFzKSxsd2Q9Mixjb2w9MikKYGBgCgoKCiMjIFNJUiBuYSBpbmZlcsOqbmNpYSBCYXllc2lhbmEKCkRhZG9zID0gJFx7MSwwLDEsMSwwLDEsMCwwLDEsMVx9JAoKTW9kZWxvOiAkeF9pfFx0aGV0YSBcIGlpZCBcIEJlcm5vdWxsaShcdGhldGEpJCBmb3IgJGk9MSxcbGRvdHMsbiQgYW5kICRcdGhldGEgXGluICgwLDEpJC4KClByaW9yaTogJFx0aGV0YSBcc2ltIEJldGEoNiwyKSQsIHN1Y2ggdGhhdCAkRShcdGhldGEpPTAuNzUkIGFuZCAkVihcdGhldGEpPTAuMTQ0XjIkLgoKUG9zdGVyaW9yaSAoZGlzdHJpYnVpw6fDo28gZGUgaW50ZXJlc3NlKQokJCAgICAgCnAoXHRoZXRhfGRhZG9zKSA9IFxmcmFjMWsgXHRoZXRhXns2LTF9KDEtXHRoZXRhKV57Mi0xfSBcdGhldGFeNigxLVx0aGV0YSleNAokJApOb3ZhbWVudGUsICRrJCDDqSBhIGNvbnN0YW50ZSBkZSBub3JtYWxpemHDp8Ojby4KCkRpc3RyaWJ1acOnw6NvIHByb3Bvc3RhOiAkQmV0YSgyLDIpJAoKVm9jw6ogcG9kZSB0cm9jYXIgZXNzYSBwcm9wb3N0YSBwb3IgcXVhbHF1ZXIgZGlzdHJpYnVpw6fDo28gbm8gaW50ZXJ2YWxvICQoMCwxKSQuCgpgYGB7cn0KTSAgICAgPSAxMDAwMApOICAgICA9IDEwMDAwCnRoZXRhID0gcmJldGEoTSwyLDIpCncgPSBkYmV0YSh0aGV0YSw2LDIpKnRoZXRhXjYqKDEtdGhldGEpXjQvZGJldGEodGhldGEsMiwyKQp3ID0gdy9zdW0odykgCnRoZXRhc3RhciA9IHNhbXBsZSh0aGV0YSxzaXplPU4scmVwbGFjZT1UUlVFLHByb2I9dykKCnBhcihtZnJvdz1jKDEsMikpCnRoZXRhcyA9IHNlcSgwLDEsbGVuZ3RoPTEwMDApCmhpc3QodGhldGEscHJvYj1UUlVFLG1haW49IkFtb3N0cmFzIGRhIHByb3Bvc3RhIix5bGltPWMoMCw0KSx4bGltPWMoMCwxKSkKbGluZXModGhldGFzLGRiZXRhKHRoZXRhcywxMiw2KSxsd2Q9Mixjb2w9MikKaGlzdCh0aGV0YXN0YXIseGxhYj0idGhldGEiLHByb2I9VFJVRSxtYWluPSJBbW9zdHJhcyBkYSBwb3N0ZXJpb3JpIix5bGltPWMoMCw0KSx4bGltPWMoMCwxKSkKbGluZXModGhldGFzLGRiZXRhKHRoZXRhcywxMiw2KSxsd2Q9Mixjb2w9MikKICAgICAgCmMobWVhbih0aGV0YXN0YXIpLCBzcXJ0KHZhcih0aGV0YXN0YXIpKSkKcXVhbnRpbGUodGhldGFzdGFyLGMoMC4wMjUsMC41LDAuOTc1KSkKYGBgCiAgIAogICAKIyMgVmFtb3MgcmVwZXRpciAyLjMpIG1hcyB1c2FuZG8gYSBwcmlvcmkgY29tbyBwcm9wb3N0YQoKRGlzdHJpYnVpw6fDo28gcHJvcG9zdGEgPSBQcmlvcmk6ICRcdGhldGEgXHNpbSBCZXRhKDYsMikkLgoKYGBge3J9Ck0gICAgID0gMTAwMDAKTiAgICAgPSAxMDAwMAp0aGV0YSA9IHJiZXRhKE0sNiwyKQp3ID0gdGhldGFeNiooMS10aGV0YSleNAp3ID0gdy9zdW0odykgCnRoZXRhc3RhciA9IHNhbXBsZSh0aGV0YSxzaXplPU4scmVwbGFjZT1UUlVFLHByb2I9dykKCnBhcihtZnJvdz1jKDEsMikpCnRoZXRhcyA9IHNlcSgwLDEsbGVuZ3RoPTEwMDApCmhpc3QodGhldGEscHJvYj1UUlVFLG1haW49IkFtb3N0cmFzIGRhIHByb3Bvc3RhIix5bGltPWMoMCw0KSx4bGltPWMoMCwxKSkKbGluZXModGhldGFzLGRiZXRhKHRoZXRhcywxMiw2KSxsd2Q9Mixjb2w9MikKaGlzdCh0aGV0YXN0YXIseGxhYj0idGhldGEiLHByb2I9VFJVRSxtYWluPSJBbW9zdHJhcyBkYSBwb3N0ZXJpb3JpIiwseWxpbT1jKDAsNCkseGxpbT1jKDAsMSkpCmxpbmVzKHRoZXRhcyxkYmV0YSh0aGV0YXMsMTIsNiksbHdkPTIsY29sPTIpCiAgICAgIApjKG1lYW4odGhldGFzdGFyKSwgc3FydCh2YXIodGhldGFzdGFyKSkpCnF1YW50aWxlKHRoZXRhc3RhcixjKDAuMDI1LDAuNSwwLjk3NSkpCmBgYAoKIyMgVmFtb3MgcmVwZXRpciAyLjQpIG1hcyBhZ29yYSBjb20gTVVJVE8gbWFpcyBkYWRvcyBlIHByaW9yaSBNVUlUTyBsb25nZSBkYSB2ZXJvc3NpbWlsaGFuw6dhCiAgIApEYWRvcyA9ICRcezEsMCwxLDEsMCwxLDAsMCwxLDEsMSwwLDEsMSwwLDEsMCwwLDEsMSwxLDAsMSwxLDAsMSwwLDAsMSwxXH0kCiAgIAkgICAgICAKUHJpb3JpOiAkXHRoZXRhIFxzaW0gQmV0YSgzLDEyKSQKClZlcm9zc2ltaWxoYW7Dp2E6ICBtdWl0byBtYWlzIGNvbmNlbnRyYWRhIQokJApMKFx0aGV0YXxkYWRvcykgPSBcdGhldGFeezE4fSgxLVx0aGV0YSleezEyfQkgICAgICAKJCQKYGBge3J9IApNICAgICA9IDEwMDAwCk4gICAgID0gMTAwMDAKdGhldGEgPSByYmV0YShNLDMsMTIpCncgPSB0aGV0YV4oMTgpKigxLXRoZXRhKV4oMTIpCncgPSB3L3N1bSh3KSAKdGhldGFzdGFyID0gc2FtcGxlKHRoZXRhLHNpemU9TixyZXBsYWNlPVRSVUUscHJvYj13KQoKcGFyKG1mcm93PWMoMSwyKSkKdGhldGFzID0gc2VxKDAsMSxsZW5ndGg9MTAwMCkKaGlzdCh0aGV0YSxwcm9iPVRSVUUsbWFpbj0iQW1vc3RyYXMgZGEgcHJvcG9zdGEiLHlsaW09YygwLDYpLHhsaW09YygwLDEpKQpsaW5lcyh0aGV0YXMsZGJldGEodGhldGFzLDIxLDI0KSxsd2Q9Mixjb2w9MikKaGlzdCh0aGV0YXN0YXIseGxhYj0idGhldGEiLHByb2I9VFJVRSxtYWluPSJBbW9zdHJhcyBkYSBwb3N0ZXJpb3JpIix5bGltPWMoMCw2KSx4bGltPWMoMCwxKSkKbGluZXModGhldGFzLGRiZXRhKHRoZXRhcywyMSwyNCksbHdkPTIsY29sPTIpCiAgICAgIApjKG1lYW4odGhldGFzdGFyKSwgc3FydCh2YXIodGhldGFzdGFyKSkpCnF1YW50aWxlKHRoZXRhc3RhcixjKDAuMDI1LDAuNSwwLjk3NSkpCmBgYAo=