6.5 Monte Carlo Markov Chain (MCMC)

Teori Resiko


Kontak : \(\downarrow\)
Email
Instagram https://www.instagram.com/claraevania/
RPubs https://rpubs.com/claradellaevania/

Sebenarnya materi yang dari web untuk 6.5 itu masih sedang dalam penulisan dan belum selesai dalam pengeditan. Jadi apa yang ditulis disini hanya memberikan gambaran besarnya saja.

Ide dari teknik Monte Carlo bergantung pada hukum bilangan besar (yang menjamin konvergensi rata-rata terhadap integral) dan teorema limit pusat (yang digunakan untuk mengukur ketidakpastian dalam perhitungan). Perlu diingat kembali jika (\(X_i\)) adalah urutan ke-i dari variabel acak dengan distribusi F, maka

\[ \frac{1}{\sqrt{n}}\left(\sum_{i=1}^n h(X_i)-\int h(x)dF(x)\right)\overset{\mathcal{L}}{\rightarrow }\mathcal{N}(0,\sigma^2),\text{ as }n\rightarrow\infty , \]

atau beberapa varian \(σ^2>0\) . Namun sebenarnya, teorema ergodik dapat digunakan untuk melemahkan hasil sebelumnya, karena independensi variabel tidak diperlukan. Lebih tepatnya, jika (\(X_i\)) adalah Proses Markov dengan ukuran invarian \(μ\) , di bawah beberapa asumsi teknis tambahan, maka dapat diperoleh

\[ \frac{1}{\sqrt{n}}\left(\sum_{i=1}^n h(X_i)-\int h(x)d\mu(x)\right)\overset{\mathcal{L}}{\rightarrow }\mathcal{N}(0,\sigma_\star^2),\text{ as }n\rightarrow\infty. \]

untuk beberapa varian \(σ^2_⋆>0\) .

Oleh karena itu, dari sifat ini, dapat melihat bahwa tidak selalu mungkin untuk menghasilkan nilai-nilai independen dari F , tetapi untuk menghasilkan proses Markov dengan ukuran invarian F , dan untuk mempertimbangkan rata-rata dari proses (tidak harus independen).

Dengan mempertimbangkan kasus vektor Gaussian terkendala: kami ingin menghasilkan pasangan acak dari vektor acak \(X\) , tetapi kami hanya tertarik pada kasus di mana jumlah komposisinya cukup besar, yang dapat ditulis \(X^T1>m\) untuk nilai nyata \(m\) . Tentu saja, dimungkinkan untuk menggunakan algoritme terima-tolak, tetapi kami telah melihat bahwa ini mungkin sangat tidak efisien. Satu dapat menggunakan Metropolis Hastingsand Gibbs sampler untuk menghasilkan proses Markov dengan ukuran invarian tersebut.

6.5.1 Metropolis Hastings

Algoritma agak sederhana untuk dihasilkan dari \(f\) : dapat dimulai dengan nilai layak \(x_1\) . Kemudian, pada langkah \(t\) , kita perlu menentukan kernel transisi : diberikan \(x_t\) , kita memerlukan distribusi bersyarat untuk \(X_{t+1}\) diberikan \(x_t\) . Algoritme akan bekerja dengan baik jika distribusi bersyarat itu dapat dengan mudah disimulasikan. dengan \(Ï€(â‹…|xt)\) menunjukkan probabilitas itu.

Gambarkan nilai potensial \(x^⋆_{t+1}\) , dan \(u\) , dari distribusi seragam. Selanjutnya Menghitung

\(R= \frac{f(x_{t+1}^\star)}{f(x_t)}\)

jika \(u<r\) , lalu atur \(x_{t+1}=x^⋆_t\) jika \(u≤r\) , maka atur \(x_{t+1}=x_t\)

Di sini r disebut rasio penerimaan selanjutnya menerima nilai baru dengan probabilitas r (atau sebenarnya yang terkecil antara 1 dan r karena r dapat melebihi 1 ).

Misalnya, asumsikan bahwa \(f(⋅|xt)\) seragam pada \([x_t−ε,x_t+ε]\) untuk beberapa \(ε>0\) , dan di mana$ $f (distribusi target kita) adalah \(N(0,1)\) . Kami tidak akan pernah menarik dari \(f\) , tetapi kami akan menggunakannya untuk menghitung rasio penerimaan kami di setiap langkah.

metrop1 <- function(n=1000,eps=0.5){
 vec <- matrix(NA, n, 3)
 x=0
 vec[1] <- x
 for (i in 2:n) {
 innov <- runif(1,-eps,eps)
 mov <- x+innov
 R <- min(1,dnorm(mov)/dnorm(x))
 u <- runif(1)
 if (u < R) x <- mov
 vec[i,] <- c(x,mov,R)
 }
return(vec)}
#install.packages('gifski')
#if (packageVersion('knitr') < '1.20.14') {
#  remotes::install_github('yihui/knitr')
#}
vec <- metrop1(25)
u=seq(-3,3,by=.01)
pic_ani = function(k){
  plot(1:k,vec[1:k,1],pch=19,xlim=c(0,25),ylim=c(-2,2),xlab="",ylab="")
    if(vec[k+1,1]==vec[k+1,2]) points(k+1,vec[k+1,1],col="blue",pch=19)
    if(vec[k+1,1]!=vec[k+1,2]) points(k+1,vec[k+1,1],col="red",pch=19)
  points(k+1,vec[k+1,2],cex=1.5)
  arrows(k+1,vec[k,1]-.5,k+1,vec[k,1]+.5,col="green",angle=90,code = 3,length=.1)
  polygon(c(k+dnorm(u)*10,rep(k,length(u))),c(u,rev(u)),col=rgb(0,1,0,.3),
          border=NA)  
  segments(k,vec[k,1],k+dnorm(vec[k,1])*10,vec[k,1])
  segments(k,vec[k+1,2],k+dnorm(vec[k+1,2])*10,vec[k+1,2])
  text(k,2,round(vec[k+1,3],digits=3))
}
for (k in 2:23) {pic_ani(k)}

Selanjutnya dapat menggunakan simulasi, maka didapat

vec <- metrop1(10000)
simx <- vec[1000:10000,1]
par(mfrow=c(1,4))
plot(simx,type="l")
hist(simx,probability = TRUE,col="light blue",border="white")
lines(u,dnorm(u),col="red")
qqnorm(simx)
acf(simx,lag=100,lwd=2,col="light blue")

6.5.2 Gibbs Sampler

Dapat mempertimbangkan beberapa vektor \(X=(X_1,⋯,X_d)\) dengan komponen independen, \(X_i∼E(λ_i)\) . Selanjutnya mengambil sampel untuk sampel dari \(X\) yang diberikan \(X^T1>s\) untuk beberapa ambang batas \(s>0\) .

  • beberapa titik awal x0 ,
  • Mengambil secara acak \(i∈{1,⋯,d}\)
  • \(X_i\) mengingat \(X_i>s−x^T_{(−i)}1\) berdistribusi Eksponensial \(E(λ_i)\)
  • Menggambar \(Y∼E(λ_i)\) dan atur \(x_i=y+(s−x^T_{(−i)}1)_+\) hingga \(x^T_{(−i)}1+x_i>s\)
sim <- NULL
 lambda <- c(1,2)
 X <- c(3,3)
 s <- 5
 for(k in 1:1000){
 i <- sample(1:2,1)
 X[i] <- rexp(1,lambda[i])+max(0,s-sum(X[-i]))
 while(sum(X)<s){
 X[i] <- rexp(1,lambda[i])+max(0,s-sum(X[-i])) }
 sim <- rbind(sim,X) }
 
plot(sim,xlim=c(1,11),ylim=c(0,4.3))
polygon(c(-1,-1,6),c(-1,6,-1),col="red",density=15,border=NA)
abline(5,-1,col="red")

Konstruksi urutan (algoritma MCMC bersifat iteratif) dapat divisualisasikan di bawah ini

lambda <- c(1,2)
X <- c(3,3)
sim <- X
s <- 5
for(k in 1:100){
 set.seed(k)
 i <- sample(1:2,1)
 X[i] <- rexp(1,lambda[i])+max(0,s-sum(X[-i]))
 while(sum(X)<s){
 X[i] <- rexp(1,lambda[i])+max(0,s-sum(X[-i])) }
 sim <- rbind(sim,X) }
pic_ani = function(n){
plot(sim[1:n,],xlim=c(1,11),ylim=c(0,5),xlab="",ylab="")
i=which(apply(sim[(n-1):n,],2,diff)==0)
if(i==1) abline(v=sim[n,1],col="grey")
if(i==2) abline(h=sim[n,2],col="grey")
if(n>=1) points(sim[n,1],sim[n,2],pch=19,col="blue",cex=1.4)
if(n>=2) points(sim[n-1,1],sim[n-1,2],pch=19,col="red",cex=1.4)
polygon(c(-1,-1,6),c(-1,6,-1),col="red",density=15,border=NA)
abline(5,-1,col="red")
}
for (i in 2:100) {pic_ani(i)}

LS0tDQp0aXRsZTogIjYuNSBNb250ZSBDYXJsbyBNYXJrb3YgQ2hhaW4gKE1DTUMpIg0Kc3VidGl0bGU6ICJUZW9yaSBSZXNpa28iDQphdXRob3I6ICJDbGFyYSBEZWxsYSBFdmFuaWEgKDIwMjA0OTIwMDE4KSINCmRhdGU6ICAiYHIgZm9ybWF0KFN5cy5EYXRlKCksICclQiAlZCwgJVknKWAiDQpvdXRwdXQ6DQogIHJtZGZvcm1hdHM6OnJvYm9ib29rOiAgICMgaHR0cHM6Ly9naXRodWIuY29tL2p1YmEvcm1kZm9ybWF0cw0KICAgIHNlbGZfY29udGFpbmVkOiB0cnVlDQogICAgdGh1bWJuYWlsczogdHJ1ZQ0KICAgIGxpZ2h0Ym94OiB0cnVlDQogICAgZ2FsbGVyeTogdHJ1ZQ0KICAgIGxpYl9kaXI6IGxpYnMNCiAgICBkZl9wcmludDogInBhZ2VkIg0KICAgIGNvZGVfZm9sZGluZzogInNob3ciDQogICAgY29kZV9kb3dubG9hZDogeWVzDQogICAgY3NzOiAic3R5bGUuY3NzIg0KDQotLS0NCg0KDQpgYGB7ciBpbmNsdWRlPUZBTFNFfQ0Ka25pdHI6Om9wdHNfY2h1bmskc2V0KGNsYXNzLnNvdXJjZSA9ICJub2NvcHkiLA0KICAgICAgICAgICAgICAgICAgICAgIGNsYXNzLm91dHB1dCA9ICJub2NvcHkiLA0KICAgICAgICAgICAgICAgICAgICAgIG1lc3NhZ2UgPSBGLA0KICAgICAgICAgICAgICAgICAgICAgIHdhcm5pbmcgPSBGKQ0KYGBgDQoNCg0KPGJyPg0KDQoNCjxpbWcgc3R5bGU9ImZsb2F0OiByaWdodDsgbWFyZ2luOiAtNTBweCA1MHB4IDBweCA1MHB4OyB3aWR0aDoyNSUiIHNyYz0ibWUuanBlZyIvPiANCg0KfA0KOi0tLS0gfDotLS0tDQpLb250YWt8IDogJFxkb3duYXJyb3ckDQpFbWFpbHwgY2xhcmEuZXZhbmlhQHN0dWRlbnQubWF0YW5hdW5pdmVyc2l0eS5hYy5pZA0KSW5zdGFncmFtIHwgaHR0cHM6Ly93d3cuaW5zdGFncmFtLmNvbS9jbGFyYWV2YW5pYS8gDQpSUHVicyAgfCBodHRwczovL3JwdWJzLmNvbS9jbGFyYWRlbGxhZXZhbmlhLyANCg0KU2ViZW5hcm55YSBtYXRlcmkgeWFuZyBkYXJpIHdlYiB1bnR1ayA2LjUgaXR1IG1hc2loIHNlZGFuZyBkYWxhbSBwZW51bGlzYW4gZGFuIGJlbHVtIHNlbGVzYWkgZGFsYW0gcGVuZ2VkaXRhbi4gSmFkaSBhcGEgeWFuZyBkaXR1bGlzIGRpc2luaSBoYW55YSBtZW1iZXJpa2FuIGdhbWJhcmFuIGJlc2FybnlhIHNhamEuDQoNCklkZSBkYXJpIHRla25payBNb250ZSBDYXJsbyBiZXJnYW50dW5nIHBhZGEgaHVrdW0gYmlsYW5nYW4gYmVzYXIgKHlhbmcgbWVuamFtaW4ga29udmVyZ2Vuc2kgcmF0YS1yYXRhIHRlcmhhZGFwIGludGVncmFsKSBkYW4gdGVvcmVtYSBsaW1pdCBwdXNhdCAoeWFuZyBkaWd1bmFrYW4gdW50dWsgbWVuZ3VrdXIga2V0aWRha3Bhc3RpYW4gZGFsYW0gcGVyaGl0dW5nYW4pLiBQZXJsdSBkaWluZ2F0IGtlbWJhbGkgamlrYSAoJFhfaSQpIGFkYWxhaCB1cnV0YW4ga2UtaSBkYXJpIHZhcmlhYmVsIGFjYWsgZGVuZ2FuIGRpc3RyaWJ1c2kgRiwgbWFrYQ0KDQokJA0KXGZyYWN7MX17XHNxcnR7bn19XGxlZnQoXHN1bV97aT0xfV5uIGgoWF9pKS1caW50IGgoeClkRih4KVxyaWdodClcb3ZlcnNldHtcbWF0aGNhbHtMfX17XHJpZ2h0YXJyb3cgfVxtYXRoY2Fse059KDAsXHNpZ21hXjIpLFx0ZXh0eyBhcyB9blxyaWdodGFycm93XGluZnR5ICwNCiQkDQoNCmF0YXUgYmViZXJhcGEgdmFyaWFuICTPg14yPjAkIC4gTmFtdW4gc2ViZW5hcm55YSwgdGVvcmVtYSBlcmdvZGlrIGRhcGF0IGRpZ3VuYWthbiB1bnR1ayBtZWxlbWFoa2FuIGhhc2lsIHNlYmVsdW1ueWEsIGthcmVuYSBpbmRlcGVuZGVuc2kgdmFyaWFiZWwgdGlkYWsgZGlwZXJsdWthbi4gTGViaWggdGVwYXRueWEsIGppa2EgKCRYX2kkKSBhZGFsYWggUHJvc2VzIE1hcmtvdiBkZW5nYW4gdWt1cmFuIGludmFyaWFuICTOvCQgLCBkaSBiYXdhaCBiZWJlcmFwYSBhc3Vtc2kgdGVrbmlzIHRhbWJhaGFuLCBtYWthIGRhcGF0IGRpcGVyb2xlaCANCg0KDQokJA0KXGZyYWN7MX17XHNxcnR7bn19XGxlZnQoXHN1bV97aT0xfV5uIGgoWF9pKS1caW50IGgoeClkXG11KHgpXHJpZ2h0KVxvdmVyc2V0e1xtYXRoY2Fse0x9fXtccmlnaHRhcnJvdyB9XG1hdGhjYWx7Tn0oMCxcc2lnbWFfXHN0YXJeMiksXHRleHR7IGFzIH1uXHJpZ2h0YXJyb3dcaW5mdHkuDQokJA0KDQp1bnR1ayBiZWJlcmFwYSB2YXJpYW4gJM+DXjJf4ouGPjAkIC4NCg0KT2xlaCBrYXJlbmEgaXR1LCBkYXJpIHNpZmF0IGluaSwgZGFwYXQgbWVsaWhhdCBiYWh3YSB0aWRhayBzZWxhbHUgbXVuZ2tpbiB1bnR1ayBtZW5naGFzaWxrYW4gbmlsYWktbmlsYWkgaW5kZXBlbmRlbiBkYXJpIEYgLCB0ZXRhcGkgdW50dWsgbWVuZ2hhc2lsa2FuIHByb3NlcyBNYXJrb3YgZGVuZ2FuIHVrdXJhbiBpbnZhcmlhbiBGICwgZGFuIHVudHVrIG1lbXBlcnRpbWJhbmdrYW4gcmF0YS1yYXRhIGRhcmkgcHJvc2VzICh0aWRhayBoYXJ1cyBpbmRlcGVuZGVuKS4NCg0KRGVuZ2FuIG1lbXBlcnRpbWJhbmdrYW4ga2FzdXMgdmVrdG9yIEdhdXNzaWFuIHRlcmtlbmRhbGE6IGthbWkgaW5naW4gbWVuZ2hhc2lsa2FuIHBhc2FuZ2FuIGFjYWsgZGFyaSB2ZWt0b3IgYWNhayAkWCQgLCB0ZXRhcGkga2FtaSBoYW55YSB0ZXJ0YXJpayBwYWRhIGthc3VzIGRpIG1hbmEganVtbGFoIGtvbXBvc2lzaW55YSBjdWt1cCBiZXNhciwgeWFuZyBkYXBhdCBkaXR1bGlzICRYXlQxPm0kIHVudHVrIG5pbGFpIG55YXRhICRtJCAuIFRlbnR1IHNhamEsIGRpbXVuZ2tpbmthbiB1bnR1ayBtZW5nZ3VuYWthbiBhbGdvcml0bWUgdGVyaW1hLXRvbGFrLCB0ZXRhcGkga2FtaSB0ZWxhaCBtZWxpaGF0IGJhaHdhIGluaSBtdW5na2luIHNhbmdhdCB0aWRhayBlZmlzaWVuLiBTYXR1IGRhcGF0IG1lbmdndW5ha2FuIE1ldHJvcG9saXMgSGFzdGluZ3NhbmQgR2liYnMgc2FtcGxlciB1bnR1ayBtZW5naGFzaWxrYW4gcHJvc2VzIE1hcmtvdiBkZW5nYW4gdWt1cmFuIGludmFyaWFuIHRlcnNlYnV0Lg0KDQojIyA2LjUuMSBNZXRyb3BvbGlzIEhhc3RpbmdzDQoNCkFsZ29yaXRtYSBhZ2FrIHNlZGVyaGFuYSB1bnR1ayBkaWhhc2lsa2FuIGRhcmkgJGYkIDogZGFwYXQgZGltdWxhaSBkZW5nYW4gbmlsYWkgbGF5YWsgJHhfMSQgLiBLZW11ZGlhbiwgcGFkYSBsYW5na2FoICR0JCAsIGtpdGEgcGVybHUgbWVuZW50dWthbiBrZXJuZWwgdHJhbnNpc2kgOiBkaWJlcmlrYW4gJHhfdCQgLCBraXRhIG1lbWVybHVrYW4gZGlzdHJpYnVzaSBiZXJzeWFyYXQgdW50dWsgJFhfe3QrMX0kIGRpYmVyaWthbiAkeF90JCAuIEFsZ29yaXRtZSBha2FuIGJla2VyamEgZGVuZ2FuIGJhaWsgamlrYSBkaXN0cmlidXNpIGJlcnN5YXJhdCBpdHUgZGFwYXQgZGVuZ2FuIG11ZGFoIGRpc2ltdWxhc2lrYW4uIGRlbmdhbiAkz4Ao4ouFfHh0KSQgbWVudW5qdWtrYW4gcHJvYmFiaWxpdGFzIGl0dS4NCg0KR2FtYmFya2FuIG5pbGFpIHBvdGVuc2lhbCAkeF7ii4Zfe3QrMX0kICwgZGFuICR1JCAsIGRhcmkNCmRpc3RyaWJ1c2kgc2VyYWdhbS4gU2VsYW5qdXRueWEgTWVuZ2hpdHVuZw0KIA0KJFI9ICBcZnJhY3tmKHhfe3QrMX1eXHN0YXIpfXtmKHhfdCl9JA0KDQoNCmppa2EgJHU8ciQgLCBsYWx1IGF0dXIgJHhfe3QrMX09eF7ii4ZfdCQNCmppa2EgJHXiiaRyJCAsIG1ha2EgYXR1ciAkeF97dCsxfT14X3QkDQoNCg0KRGkgc2luaSByIGRpc2VidXQgcmFzaW8gcGVuZXJpbWFhbiBzZWxhbmp1dG55YSBtZW5lcmltYSBuaWxhaSBiYXJ1IGRlbmdhbiBwcm9iYWJpbGl0YXMgciAoYXRhdSBzZWJlbmFybnlhIHlhbmcgdGVya2VjaWwgYW50YXJhIDEgZGFuIHIga2FyZW5hIHIgZGFwYXQgbWVsZWJpaGkgMSApLg0KDQpNaXNhbG55YSwgYXN1bXNpa2FuIGJhaHdhICRmKOKLhXx4dCkkIHNlcmFnYW0gcGFkYSAkW3hfdOKIks61LHhfdCvOtV0kIHVudHVrIGJlYmVyYXBhICTOtT4wJCAsIGRhbiBkaSBtYW5hJCAkZiAoZGlzdHJpYnVzaSB0YXJnZXQga2l0YSkgYWRhbGFoICROKDAsMSkkIC4gS2FtaSB0aWRhayBha2FuIHBlcm5haCBtZW5hcmlrIGRhcmkgJGYkICwgdGV0YXBpIGthbWkgYWthbiBtZW5nZ3VuYWthbm55YSB1bnR1ayBtZW5naGl0dW5nIHJhc2lvIHBlbmVyaW1hYW4ga2FtaSBkaSBzZXRpYXAgbGFuZ2thaC4NCg0KYGBge3J9DQptZXRyb3AxIDwtIGZ1bmN0aW9uKG49MTAwMCxlcHM9MC41KXsNCiB2ZWMgPC0gbWF0cml4KE5BLCBuLCAzKQ0KIHg9MA0KIHZlY1sxXSA8LSB4DQogZm9yIChpIGluIDI6bikgew0KIGlubm92IDwtIHJ1bmlmKDEsLWVwcyxlcHMpDQogbW92IDwtIHgraW5ub3YNCiBSIDwtIG1pbigxLGRub3JtKG1vdikvZG5vcm0oeCkpDQogdSA8LSBydW5pZigxKQ0KIGlmICh1IDwgUikgeCA8LSBtb3YNCiB2ZWNbaSxdIDwtIGMoeCxtb3YsUikNCiB9DQpyZXR1cm4odmVjKX0NCmBgYA0KDQpgYGB7cn0NCiNpbnN0YWxsLnBhY2thZ2VzKCdnaWZza2knKQ0KI2lmIChwYWNrYWdlVmVyc2lvbigna25pdHInKSA8ICcxLjIwLjE0Jykgew0KIyAgcmVtb3Rlczo6aW5zdGFsbF9naXRodWIoJ3lpaHVpL2tuaXRyJykNCiN9DQp2ZWMgPC0gbWV0cm9wMSgyNSkNCnU9c2VxKC0zLDMsYnk9LjAxKQ0KcGljX2FuaSA9IGZ1bmN0aW9uKGspew0KICBwbG90KDE6ayx2ZWNbMTprLDFdLHBjaD0xOSx4bGltPWMoMCwyNSkseWxpbT1jKC0yLDIpLHhsYWI9IiIseWxhYj0iIikNCiAgICBpZih2ZWNbaysxLDFdPT12ZWNbaysxLDJdKSBwb2ludHMoaysxLHZlY1trKzEsMV0sY29sPSJibHVlIixwY2g9MTkpDQogICAgaWYodmVjW2srMSwxXSE9dmVjW2srMSwyXSkgcG9pbnRzKGsrMSx2ZWNbaysxLDFdLGNvbD0icmVkIixwY2g9MTkpDQogIHBvaW50cyhrKzEsdmVjW2srMSwyXSxjZXg9MS41KQ0KICBhcnJvd3MoaysxLHZlY1trLDFdLS41LGsrMSx2ZWNbaywxXSsuNSxjb2w9ImdyZWVuIixhbmdsZT05MCxjb2RlID0gMyxsZW5ndGg9LjEpDQogIHBvbHlnb24oYyhrK2Rub3JtKHUpKjEwLHJlcChrLGxlbmd0aCh1KSkpLGModSxyZXYodSkpLGNvbD1yZ2IoMCwxLDAsLjMpLA0KICAgICAgICAgIGJvcmRlcj1OQSkgIA0KICBzZWdtZW50cyhrLHZlY1trLDFdLGsrZG5vcm0odmVjW2ssMV0pKjEwLHZlY1trLDFdKQ0KICBzZWdtZW50cyhrLHZlY1trKzEsMl0saytkbm9ybSh2ZWNbaysxLDJdKSoxMCx2ZWNbaysxLDJdKQ0KICB0ZXh0KGssMixyb3VuZCh2ZWNbaysxLDNdLGRpZ2l0cz0zKSkNCn0NCmZvciAoayBpbiAyOjIzKSB7cGljX2FuaShrKX0NCmBgYA0KDQpTZWxhbmp1dG55YSBkYXBhdCBtZW5nZ3VuYWthbiBzaW11bGFzaSwgbWFrYSBkaWRhcGF0DQoNCg0KYGBge3J9DQp2ZWMgPC0gbWV0cm9wMSgxMDAwMCkNCnNpbXggPC0gdmVjWzEwMDA6MTAwMDAsMV0NCnBhcihtZnJvdz1jKDEsNCkpDQpwbG90KHNpbXgsdHlwZT0ibCIpDQpoaXN0KHNpbXgscHJvYmFiaWxpdHkgPSBUUlVFLGNvbD0ibGlnaHQgYmx1ZSIsYm9yZGVyPSJ3aGl0ZSIpDQpsaW5lcyh1LGRub3JtKHUpLGNvbD0icmVkIikNCnFxbm9ybShzaW14KQ0KYWNmKHNpbXgsbGFnPTEwMCxsd2Q9Mixjb2w9ImxpZ2h0IGJsdWUiKQ0KYGBgDQoNCiMjIDYuNS4yIEdpYmJzIFNhbXBsZXINCg0KRGFwYXQgbWVtcGVydGltYmFuZ2thbiBiZWJlcmFwYSB2ZWt0b3IgJFg9KFhfMSzii68sWF9kKSQgZGVuZ2FuIGtvbXBvbmVuIGluZGVwZW5kZW4sICRYX2niiLxFKM67X2kpJCAuIFNlbGFuanV0bnlhIG1lbmdhbWJpbCBzYW1wZWwgdW50dWsgc2FtcGVsIGRhcmkgJFgkIHlhbmcgZGliZXJpa2FuICRYXlQxPnMkIHVudHVrIGJlYmVyYXBhIGFtYmFuZyBiYXRhcyAkcz4wJCAuDQoNCi0gICBiZWJlcmFwYSB0aXRpayBhd2FsIHgwICwNCi0gICBNZW5nYW1iaWwgc2VjYXJhIGFjYWsgJGniiIh7MSzii68sZH0kDQotICAgJFhfaSQgbWVuZ2luZ2F0ICRYX2k+c+KIknheVF97KOKIkmkpfTEkIGJlcmRpc3RyaWJ1c2kgRWtzcG9uZW5zaWFsICRFKM67X2kpJA0KLSAgIE1lbmdnYW1iYXIgJFniiLxFKM67X2kpJCBkYW4gYXR1ciAkeF9pPXkrKHPiiJJ4XlRfeyjiiJJpKX0xKV8rJCBoaW5nZ2EgJHheVF97KOKIkmkpfTEreF9pPnMkDQoNCg0KYGBge3J9DQpzaW0gPC0gTlVMTA0KIGxhbWJkYSA8LSBjKDEsMikNCiBYIDwtIGMoMywzKQ0KIHMgPC0gNQ0KIGZvcihrIGluIDE6MTAwMCl7DQogaSA8LSBzYW1wbGUoMToyLDEpDQogWFtpXSA8LSByZXhwKDEsbGFtYmRhW2ldKSttYXgoMCxzLXN1bShYWy1pXSkpDQogd2hpbGUoc3VtKFgpPHMpew0KIFhbaV0gPC0gcmV4cCgxLGxhbWJkYVtpXSkrbWF4KDAscy1zdW0oWFstaV0pKSB9DQogc2ltIDwtIHJiaW5kKHNpbSxYKSB9DQogDQpwbG90KHNpbSx4bGltPWMoMSwxMSkseWxpbT1jKDAsNC4zKSkNCnBvbHlnb24oYygtMSwtMSw2KSxjKC0xLDYsLTEpLGNvbD0icmVkIixkZW5zaXR5PTE1LGJvcmRlcj1OQSkNCmFibGluZSg1LC0xLGNvbD0icmVkIikNCg0KYGBgDQoNCktvbnN0cnVrc2kgdXJ1dGFuIChhbGdvcml0bWEgTUNNQyBiZXJzaWZhdCBpdGVyYXRpZikgZGFwYXQgZGl2aXN1YWxpc2FzaWthbiBkaSBiYXdhaCBpbmkNCg0KYGBge3J9DQpsYW1iZGEgPC0gYygxLDIpDQpYIDwtIGMoMywzKQ0Kc2ltIDwtIFgNCnMgPC0gNQ0KZm9yKGsgaW4gMToxMDApew0KIHNldC5zZWVkKGspDQogaSA8LSBzYW1wbGUoMToyLDEpDQogWFtpXSA8LSByZXhwKDEsbGFtYmRhW2ldKSttYXgoMCxzLXN1bShYWy1pXSkpDQogd2hpbGUoc3VtKFgpPHMpew0KIFhbaV0gPC0gcmV4cCgxLGxhbWJkYVtpXSkrbWF4KDAscy1zdW0oWFstaV0pKSB9DQogc2ltIDwtIHJiaW5kKHNpbSxYKSB9DQpwaWNfYW5pID0gZnVuY3Rpb24obil7DQpwbG90KHNpbVsxOm4sXSx4bGltPWMoMSwxMSkseWxpbT1jKDAsNSkseGxhYj0iIix5bGFiPSIiKQ0KaT13aGljaChhcHBseShzaW1bKG4tMSk6bixdLDIsZGlmZik9PTApDQppZihpPT0xKSBhYmxpbmUodj1zaW1bbiwxXSxjb2w9ImdyZXkiKQ0KaWYoaT09MikgYWJsaW5lKGg9c2ltW24sMl0sY29sPSJncmV5IikNCmlmKG4+PTEpIHBvaW50cyhzaW1bbiwxXSxzaW1bbiwyXSxwY2g9MTksY29sPSJibHVlIixjZXg9MS40KQ0KaWYobj49MikgcG9pbnRzKHNpbVtuLTEsMV0sc2ltW24tMSwyXSxwY2g9MTksY29sPSJyZWQiLGNleD0xLjQpDQpwb2x5Z29uKGMoLTEsLTEsNiksYygtMSw2LC0xKSxjb2w9InJlZCIsZGVuc2l0eT0xNSxib3JkZXI9TkEpDQphYmxpbmUoNSwtMSxjb2w9InJlZCIpDQp9DQpmb3IgKGkgaW4gMjoxMDApIHtwaWNfYW5pKGkpfQ0KYGBgDQoNCiMgUmVmZXJlbnNpIA0KLSBodHRwczovL29wZW5hY3R0ZXh0cy5naXRodWIuaW8vTG9zcy1EYXRhLUFuYWx5dGljcy9DaGFwRnJlcXVlbmN5LU1vZGVsaW5nLmh0bWwjUzpnb29kbmVzcy1vZi1maXQ=