A Markov chain
example
Suppose that \(P\) is the transition
matrix of a discrete-time (homogeneous) Markov chain with 5 states:
\[
P_{ij} = Pr(X_t=j|X_{t-1}=i) \qquad \forall t \ \ \mbox{and} \ \forall
i,j
\] Let \(P\) be the following
matrix: \[
P =
\begin{bmatrix}
0.6&0.1&0.1&0.1&0.1\\
0.0&0.4&0.2&0.1&0.3\\
0.1&0.0&0.5&0.3&0.1\\
0.0&0.0&0.1&0.8&0.1\\
0.1&0.1&0.1&0.0&0.7
\end{bmatrix}
\] It is easy to show that the probability of starting at state
\(i\) (time \(t=0\))and reaching state \(j\) after \(k\) (time \(t=k\)) iterations of the Markov chain is
equal to \[
Pr(X_k=j|X_0=i) = (P^k)_{ij}
\]
P = matrix(c(
0.6,0.1,0.1,0.1,0.1,
0.0,0.4,0.2,0.1,0.3,
0.1,0.0,0.5,0.3,0.1,
0.0,0.0,0.1,0.8,0.1,
0.1,0.1,0.1,0.0,0.7),5,5,byrow=TRUE)
P
## [,1] [,2] [,3] [,4] [,5]
## [1,] 0.6 0.1 0.1 0.1 0.1
## [2,] 0.0 0.4 0.2 0.1 0.3
## [3,] 0.1 0.0 0.5 0.3 0.1
## [4,] 0.0 0.0 0.1 0.8 0.1
## [5,] 0.1 0.1 0.1 0.0 0.7
For example \(P^5\) is the following
matrix
P5 = P%*%P%*%P%*%P%*%P
round(P5,4)
## [,1] [,2] [,3] [,4] [,5]
## [1,] 0.1617 0.0811 0.1803 0.3043 0.2726
## [2,] 0.1184 0.0809 0.1869 0.2942 0.3197
## [3,] 0.1098 0.0525 0.1816 0.4106 0.2456
## [4,] 0.0752 0.0415 0.1693 0.4735 0.2405
## [5,] 0.1495 0.0972 0.1828 0.2140 0.3565
such that \[
Pr(X_5=4|X_0=2) = (P^5)_{24} = 0.2942 \ > \ 0.1 = P_{24} =
Pr(X_1=4|X_0=2)
\]
Limiting/equilibrium
distribution
How does \(P^{100}\) looks like?
P1 = P
for (k in 1:100)
P1 = P1%*%P
round(P1,4)
## [,1] [,2] [,3] [,4] [,5]
## [1,] 0.1152 0.0664 0.1777 0.3574 0.2832
## [2,] 0.1152 0.0664 0.1777 0.3574 0.2832
## [3,] 0.1152 0.0664 0.1777 0.3574 0.2832
## [4,] 0.1152 0.0664 0.1777 0.3574 0.2832
## [5,] 0.1152 0.0664 0.1777 0.3574 0.2832
In other words, after a large number of iterations (\(k=100\) here) the rows of the transition
matrix converge to the same probability distribution for the 3 states:
\[
Pr(X_k=j|X_0=i) = Pr(X_k=j) = \pi_j \qquad \forall i,j.
\] The distribution \(\pi=(\pi_1,\pi_2,\pi_3)\) is known as the
limiting distribution of the Markov chain, or its
equilibrium distribution.
Markov chain and
Bayesian computation
Markov chain Monte Carlo (MCMC) schemes are Monte Carlo schemes that
build Markov chains whose equilibrium/limiting distribution is the
posterior distribution of interest. Gibbs sampler, random-walk
Metropolis-Hastings, independent Metropolis-Hastings, Hamiltonian Monte
Carlo, Delayed Metropolis, Multiple-try Metropolis are all schemes based
on well defined Markov chains structures.
Sampling the path of
the Markov chain
Let us know compare the theoretical equilibrium distribution, \[
\pi=(0.1152,0.0664,0.1777,0.3574,0.2832),
\] with its sampling approximation based on the path of the
simulated Markov chain.
xs = rep(0,100000)
x = 1
for (t in 1:100000){
x = sample(1:5,size=1,prob=P[x,])
xs[t] = x
}
xs = xs[95001:100000]
par(mfrow=c(2,1))
ts.plot(xs,ylab="State of the Markov chain",xlab="Iteration (after burn-in of size 95000)")
plot(table(xs)/5000,ylab="Relative frequency",lwd=6,xlab="States")
lines(P1[1,],type="h",col=2,lwd=2)
legend("topleft",legend=c("Theoretical limiting probability","Markov chain simulation"),col=2:1,lwd=3,bty="n")
tab = rbind(P1[1,],table(xs)/5000)
rownames(tab) = c("Theoretical","Simulated")
round(tab,4)
## 1 2 3 4 5
## Theoretical 0.1152 0.0664 0.1777 0.3574 0.2832
## Simulated 0.1110 0.0622 0.1694 0.3580 0.2994
Multiple linear
regression via Gibbs sampler
Suppose the data is \(\mbox{data}=\{(y_1,x_1),\ldots,(y_n,x_n)\}\),
where \(y_i\) is a real number and
\(x_i=(x_{i1},\dots,x_{ip})'\) is a
\(p\)-dimensional vector of explanatory
variables, also known as regressors or features or characteristics of
unit \(i\), for \(i=1,\ldots,n\). One of the most common
models for this kind of arrangement is that of a Gaussian linear
regression, where \[
y_i | x_i,\beta,\sigma^2 \sim N(x_i'\beta,\sigma^2),
\] where \(\beta=(\beta_1,\ldots,\beta_p)'\) and
\(y\)s conditionally independent given
\(x_i,\beta,\sigma^2\).
Likelihood
function
This heavily structured model (linearity, gaussianity, exogeneity,
homokedasticity) has likelihood \[
L(\beta_1,\ldots,\beta_p,\sigma^2) = \prod_{i=1}^n
(2\pi\sigma^2)^{-1/2}\exp\left\{\frac{1}{2\sigma
^2}(y_i-x_i'\beta)^2\right\},
\] for \(\beta_i \in \Re\) and
\(\sigma^2 \in \Re^+\). In matrix
notation, \[
y|X,\beta,\sigma^2 \sim N(X\beta,\sigma^2I_n),
\] where \(y=(y_1,\ldots,y_n)\),
\(X=(x_1,\ldots,x_n)'\) of
dimension \((n \times p)\) and \(I_n\) is the identity matrix of order \(n\). In this case the likelihood becomes
\[
L(\beta,\sigma^2) \propto
(\sigma^2)^{-\frac{n}{2}}\exp\left\{\frac{1}{2\sigma
^2}(y-X\beta)'(y-X\beta)\right\},
\] A further simplication leads to \[
L(\beta,\sigma^2) \propto
(\sigma^2)^{-\frac{n}{2}}\exp\left\{\frac{1}{2\sigma ^2}(\beta'
X'X \beta - 2\beta'X'y+y'y)\right\}.
\] In other words, the sufficient statistics for the Gaussian
linear regression are the cross-products \(S_{xx}=X'X\), \(S_{xy}=X'y\) and \(S_{yy}=y'y\).
Prior
distribution
Let us assume that, for \(i=1,\ldots,p\), \(\beta_i \sim N(0,\tau^2)\), say \(\tau=1\) or \(\tau=10\), and \(p(\sigma^2) \propto 1/\sigma^2\). This way
the prior of \((\beta,\sigma^2)\) is
\[
p(\beta,\sigma^2) \propto
\frac{1}{\sigma^2}\exp\left\{-\frac{1}{2\tau^2}\sum_{j=1}^p
\beta_j^2\right\}=
\frac{1}{\sigma^2}\exp\left\{-\frac{1}{2\tau^2} \beta'\beta\right\}
\]
Posterior
distribution
Combining the likelihood function with the prior density, we arrive
at the kernel of the posterior density \[
p(\beta,\sigma^2|\mbox{data}) \propto
\left[(\sigma^2)^{-1}\exp\left\{-\frac{1}{2\tau^2}
\beta'\beta\right\}\right] \times
\left[(\sigma^2)^{-\frac{n}{2}}\exp\left\{-\frac{1}{2\sigma
^2}(\beta' X'X \beta -
2\beta'X'y+y'y)\right\}\right].
\]
Full conditional
distributions
Suppose we do not recognize this joint posterior, but are able to
derive the full conditionals \[
p(\beta|\sigma^2,\mbox{data}) \ \ \ \mbox{and} \ \ \
p(\sigma^2|\beta,\mbox{data}).
\]
Full conditional of
\(\sigma^2\)
\[
p(\sigma^2|\beta,\mbox{data}) \propto
(\sigma^2)^{-\left(\frac{n}{2}+1\right)}
\exp\left\{-\frac{1}{2\sigma ^2}\left(\beta' X'X \beta -
2\beta'X'y+y'y\right)\right\},
\] which looks like the kernel of an following inverse-gamma
distribution \[
(\sigma^2|\beta,\mbox{data}) \sim IG\left(\frac{n}{2},\frac{\beta'
S_{xx} \beta - 2\beta'S_{xy}+S_{yy})}{2}\right) \equiv
IG\left(\frac{n}{2},\frac{(y-X\beta)'(y-X\beta)}{2}\right).
\]
Full conditional of
\(\beta\)
We can see that \[
p(\beta|\sigma^2,\mbox{data}) \propto \exp\left\{-\frac{1}{2\tau^2}
\beta'\beta\right\}
\exp\left\{-\frac{1}{2\sigma ^2}(\beta' X'X \beta -
2\beta'X'y)\right\} =
\exp\left\{-\frac{1}{2}\left(\beta'(I_p/\tau^2+X'X/\sigma^2)\beta-2\beta'X'y/\sigma^2\right)\right\},
\] which looks like the kernel of a multivariate normal
distribution with mean and variance \[
(I_p/\tau^2+X'X/\sigma^2)^{-1}X'y/\sigma^2 \ \ \ \mbox{and} \ \
\ (I_p/\tau^2+X'X/\sigma^2)^{-1},
\] or \[
(\beta|\sigma^2,\mbox{data}) \sim
N\left[(I_p/\tau^2+S_{xx}/\sigma^2)^{-1}S_{xy}/sigma^2,
(I_p/\tau^2+S_{xx}/\sigma^2)^{-1}\right].
\]
Gibbs sampler
Sample \(\beta\) from \(N\left[(I_p/\tau^2+S_{xx}/\sigma^2)^{-1}S_{xy}/\sigma^2,
(I_p/\tau^2+S_{xx}/\sigma^2)^{-1}\right]\)
Sample \(\sigma^2\) from \(IG\left[n/2,(\beta' S_{xx} \beta -
2\beta'S_{xy}+S_{yy})/2\right]\)
Return on
education
The data is a 1976 Panel Study of Income Dynamics, based on data for
the previous year, 1975. Of the 753 observations, the first 428 are for
women with positive hours worked in 1975, while the remaining 325
observations are for women who did not work for pay in 1975. A more
complete discussion of the data is found in Mroz [1987], Appendix 1.
Thomas A. Mroz (1987) The Sensitivity of an Empirical Model of Married
Women’s Hours of Work to Economic and Statistical Assumptions.
Econometrica, Vol. 55, No. 4 (July 1987), pp. 765-799. Stable URL: http://www.jstor.org/stable/1911029. The variables in
the dataset are as follows:
- LFP “A dummy variable = 1 if woman worked in 1975, else 0”;
- WHRS “Wife’s hours of work in 1975”;
- KL6 “Number of children less than 6 years old in household”;
- K618 “Number of children between ages 6 and 18 in household”;
- WA “Wife’s age”;
- WE “Wife’s educational attainment, in years”;
- WW “Wife’s average hourly earnings, in 1975 dollars”;
- RPWG “Wife’s wage reported at the time of the 1976 interview (not
the same as the 1975 estimated wage). To use the subsample with this
wage, one needs to select 1975 workers with LFP=1, then select only
those women with non-zero RPWG. Only 325 women work in 1975 and have a
non-zero RPWG in 1976.”;
- HHRS “Husband’s hours worked in 1975”;
- HA “Husband’s age”;
- HE “Husband’s educational attainment, in years”;
- HW “Husband’s wage, in 1975 dollars”;
- FAMINC “Family income, in 1975 dollars. This variable is used to
construct the property income variable.”;
- MTR “This is the marginal tax rate facing the wife, and is taken
from published federal tax tables (state and local income taxes are
excluded). The taxable income on which this tax rate is calculated
includes Social Security, if applicable to wife.”;
- WMED “Wife’s mother’s educational attainment, in years”;
- WFED “Wife’s father’s educational attainment, in years”;
- UN “Unemployment rate in county of residence, in percentage points.
This taken from bracketed ranges.”;
- CIT “Dummy variable = 1 if live in large city (SMSA), else 0”;
- AX “Actual years of wife’s previous labor market experience”;
data = read.table("http://hedibert.org/wp-content/uploads/2020/01/mroz-data.txt",header=TRUE)
attach(data)
names= c("working (W)",
"hours work (W)",
"children 6-18",
"hourly earnings (W)",
"hours worked (H)",
"wage (H)",
"marginal tax rate (W)",
"live in large city")
y = scale(log(FAMINC))
X = cbind(LFP,WHRS,K618,WW,HHRS,HW,MTR,CIT)
n = nrow(X)
p = ncol(X)
par(mfrow=c(2,4))
boxplot(y~X[,1],xlab=names[1],main="",outline=FALSE,names=c("No","Yes"),ylab="Family income (log)")
for (i in 2:(p-1))
plot(X[,i],y,xlab=names[i],ylab="Family income (log)")
boxplot(y~X[,p],xlab=names[p],main="",outline=FALSE,names=c("No","Yes"),ylab="Family income (log)")
Ordinary least
squares
fit.ols = lm(y~X-1)
summary(fit.ols)
##
## Call:
## lm(formula = y ~ X - 1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.6820 -0.1883 -0.0045 0.2341 3.0471
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## XLFP -2.267e-01 6.749e-02 -3.360 0.000820 ***
## XWHRS 2.864e-04 3.250e-05 8.812 < 2e-16 ***
## XK618 3.537e-02 1.455e-02 2.432 0.015268 *
## XWW 4.378e-02 7.672e-03 5.706 1.67e-08 ***
## XHHRS 4.553e-04 2.776e-05 16.400 < 2e-16 ***
## XHW 1.271e-01 4.311e-03 29.488 < 2e-16 ***
## XMTR -3.428e+00 1.027e-01 -33.374 < 2e-16 ***
## XCIT 1.528e-01 4.125e-02 3.704 0.000228 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5136 on 745 degrees of freedom
## Multiple R-squared: 0.7387, Adjusted R-squared: 0.7359
## F-statistic: 263.2 on 8 and 745 DF, p-value: < 2.2e-16
beta.ols = fit.ols$coef
sig.ols = sqrt(sum((y-X%*%beta.ols)^2)/(n-p))
Bayesian posterior
inference via Gibbs sampler
library("mvtnorm")
tau2 = 1.0
Sxx = t(X)%*%X
Sxy = t(X)%*%y
a = n/2
Ip = diag(1/tau2,p)
sigma2 = 1
M0 = 10000
M = 10000
niter = M0+M
draws = matrix(0,niter,p+1)
for (iter in 1:niter){
V = solve(Ip+Sxx/sigma2)
m = V%*%Sxy/sigma2
beta = matrix(rmvnorm(1,mean=m,sigma=V),p,1)
sigma2 = 1/rgamma(1,a,sum((y-X%*%beta)^2)/2)
draws[iter,] = c(beta,sqrt(sigma2))
}
draws = draws[(M0+1):niter,]
Graphical
summaries
par(mfrow=c(3,3))
for (i in 1:p){
hist(draws[,i],prob=TRUE,main=names[i],xlab=paste("beta(",i,")",sep=""))
abline(v=0,col=2,lwd=3)
abline(v=beta.ols[i],col=3,lwd=3)
}
hist(draws[,p+1],prob=TRUE,main="sigma",xlab="")
abline(v=sig.ols,col=3,lwd=3)
LS0tCnRpdGxlOiAiQ2xhc3MgbWF0ZXJpYWwiCmF1dGhvcjogIkhlZGliZXJ0IEZyZWl0YXMgTG9wZXMiCmRhdGU6ICI1LzIxLzIwMjQiCm91dHB1dDoKICBodG1sX2RvY3VtZW50OgogICAgdG9jOiB0cnVlCiAgICB0b2NfZGVwdGg6IDMKICAgIHRvY19jb2xsYXBzZWQ6IHRydWUKICAgIGNvZGVfZG93bmxvYWQ6IHllcwogICAgbnVtYmVyX3NlY3Rpb25zOiB0cnVlCi0tLQoKIyBBIE1hcmtvdiBjaGFpbiBleGFtcGxlClN1cHBvc2UgdGhhdCAkUCQgaXMgdGhlIHRyYW5zaXRpb24gbWF0cml4IG9mIGEgZGlzY3JldGUtdGltZSAoaG9tb2dlbmVvdXMpIE1hcmtvdiBjaGFpbiB3aXRoIDUgc3RhdGVzOgokJApQX3tpan0gPSBQcihYX3Q9anxYX3t0LTF9PWkpIFxxcXVhZCBcZm9yYWxsIHQgXCBcIFxtYm94e2FuZH0gXCBcZm9yYWxsIGksagokJApMZXQgJFAkIGJlIHRoZSBmb2xsb3dpbmcgbWF0cml4OgokJApQID0KXGJlZ2lue2JtYXRyaXh9CjAuNiYwLjEmMC4xJjAuMSYwLjFcXAowLjAmMC40JjAuMiYwLjEmMC4zXFwKMC4xJjAuMCYwLjUmMC4zJjAuMVxcCjAuMCYwLjAmMC4xJjAuOCYwLjFcXAowLjEmMC4xJjAuMSYwLjAmMC43ClxlbmR7Ym1hdHJpeH0KJCQKSXQgaXMgZWFzeSB0byBzaG93IHRoYXQgdGhlIHByb2JhYmlsaXR5IG9mIHN0YXJ0aW5nIGF0IHN0YXRlICRpJCAodGltZSAkdD0wJClhbmQgcmVhY2hpbmcgc3RhdGUgJGokIGFmdGVyICRrJCAodGltZSAkdD1rJCkgaXRlcmF0aW9ucyBvZiB0aGUgTWFya292IGNoYWluIGlzIGVxdWFsIHRvCiQkClByKFhfaz1qfFhfMD1pKSA9IChQXmspX3tpan0KJCQKCmBgYHtyfQpQID0gbWF0cml4KGMoCjAuNiwwLjEsMC4xLDAuMSwwLjEsCjAuMCwwLjQsMC4yLDAuMSwwLjMsCjAuMSwwLjAsMC41LDAuMywwLjEsCjAuMCwwLjAsMC4xLDAuOCwwLjEsCjAuMSwwLjEsMC4xLDAuMCwwLjcpLDUsNSxieXJvdz1UUlVFKQpQCmBgYAoKRm9yIGV4YW1wbGUgJFBeNSQgaXMgdGhlIGZvbGxvd2luZyBtYXRyaXgKYGBge3J9ClA1ID0gUCUqJVAlKiVQJSolUCUqJVAKcm91bmQoUDUsNCkKYGBgCnN1Y2ggdGhhdCAKJCQKUHIoWF81PTR8WF8wPTIpID0gKFBeNSlfezI0fSA9IDAuMjk0MiBcID4gXCAwLjEgPSBQX3syNH0gPSBQcihYXzE9NHxYXzA9MikKJCQKCiMjIExpbWl0aW5nL2VxdWlsaWJyaXVtIGRpc3RyaWJ1dGlvbgoKSG93IGRvZXMgJFBeezEwMH0kIGxvb2tzIGxpa2U/CmBgYHtyfQpQMSA9IFAKZm9yIChrIGluIDE6MTAwKQogIFAxID0gUDElKiVQCnJvdW5kKFAxLDQpCmBgYAoKSW4gb3RoZXIgd29yZHMsIGFmdGVyIGEgbGFyZ2UgbnVtYmVyIG9mIGl0ZXJhdGlvbnMgKCRrPTEwMCQgaGVyZSkgdGhlIHJvd3Mgb2YgdGhlIHRyYW5zaXRpb24gbWF0cml4IGNvbnZlcmdlIHRvIHRoZSBzYW1lIHByb2JhYmlsaXR5IGRpc3RyaWJ1dGlvbiBmb3IgdGhlIDMgc3RhdGVzOgokJApQcihYX2s9anxYXzA9aSkgPSBQcihYX2s9aikgPSBccGlfaiBccXF1YWQgXGZvcmFsbCBpLGouCiQkClRoZSBkaXN0cmlidXRpb24gJFxwaT0oXHBpXzEsXHBpXzIsXHBpXzMpJCBpcyBrbm93biBhcyB0aGUgKmxpbWl0aW5nIGRpc3RyaWJ1dGlvbiogb2YgdGhlIE1hcmtvdiBjaGFpbiwgb3IgaXRzICAqZXF1aWxpYnJpdW0gZGlzdHJpYnV0aW9uKi4KCiMjIE1hcmtvdiBjaGFpbiBhbmQgQmF5ZXNpYW4gY29tcHV0YXRpb24KCk1hcmtvdiBjaGFpbiBNb250ZSBDYXJsbyAoTUNNQykgc2NoZW1lcyBhcmUgTW9udGUgQ2FybG8gc2NoZW1lcyB0aGF0IGJ1aWxkIE1hcmtvdiBjaGFpbnMgd2hvc2UgZXF1aWxpYnJpdW0vbGltaXRpbmcgZGlzdHJpYnV0aW9uIGlzIHRoZSBwb3N0ZXJpb3IgZGlzdHJpYnV0aW9uIG9mIGludGVyZXN0LiAgR2liYnMgc2FtcGxlciwgcmFuZG9tLXdhbGsgTWV0cm9wb2xpcy1IYXN0aW5ncywgaW5kZXBlbmRlbnQgTWV0cm9wb2xpcy1IYXN0aW5ncywgSGFtaWx0b25pYW4gTW9udGUgQ2FybG8sIERlbGF5ZWQgTWV0cm9wb2xpcywgTXVsdGlwbGUtdHJ5IE1ldHJvcG9saXMgYXJlIGFsbCBzY2hlbWVzIGJhc2VkIG9uIHdlbGwgZGVmaW5lZCBNYXJrb3YgY2hhaW5zIHN0cnVjdHVyZXMuCgojIyBTYW1wbGluZyB0aGUgcGF0aCBvZiB0aGUgTWFya292IGNoYWluCkxldCB1cyBrbm93IGNvbXBhcmUgdGhlIHRoZW9yZXRpY2FsIGVxdWlsaWJyaXVtIGRpc3RyaWJ1dGlvbiwgCiQkClxwaT0oMC4xMTUyLDAuMDY2NCwwLjE3NzcsMC4zNTc0LDAuMjgzMiksCiQkCndpdGggaXRzIHNhbXBsaW5nIGFwcHJveGltYXRpb24gYmFzZWQgb24gdGhlIHBhdGggb2YgdGhlIHNpbXVsYXRlZCBNYXJrb3YgY2hhaW4uCgpgYGB7ciBmaWcud2lkdGg9NiwgZmlnLmhlaWdodD0xMCwgZmlnLmFsaWduPSdjZW50ZXInfQp4cyA9IHJlcCgwLDEwMDAwMCkKeCAgPSAxCmZvciAodCBpbiAxOjEwMDAwMCl7CiAgeCA9IHNhbXBsZSgxOjUsc2l6ZT0xLHByb2I9UFt4LF0pCiAgeHNbdF0gPSB4Cn0KeHMgPSB4c1s5NTAwMToxMDAwMDBdCgpwYXIobWZyb3c9YygyLDEpKQp0cy5wbG90KHhzLHlsYWI9IlN0YXRlIG9mIHRoZSBNYXJrb3YgY2hhaW4iLHhsYWI9Ikl0ZXJhdGlvbiAoYWZ0ZXIgYnVybi1pbiBvZiBzaXplIDk1MDAwKSIpCnBsb3QodGFibGUoeHMpLzUwMDAseWxhYj0iUmVsYXRpdmUgZnJlcXVlbmN5Iixsd2Q9Nix4bGFiPSJTdGF0ZXMiKQpsaW5lcyhQMVsxLF0sdHlwZT0iaCIsY29sPTIsbHdkPTIpCmxlZ2VuZCgidG9wbGVmdCIsbGVnZW5kPWMoIlRoZW9yZXRpY2FsIGxpbWl0aW5nIHByb2JhYmlsaXR5IiwiTWFya292IGNoYWluIHNpbXVsYXRpb24iKSxjb2w9MjoxLGx3ZD0zLGJ0eT0ibiIpCgp0YWIgPSByYmluZChQMVsxLF0sdGFibGUoeHMpLzUwMDApCnJvd25hbWVzKHRhYikgPSBjKCJUaGVvcmV0aWNhbCIsIlNpbXVsYXRlZCIpCnJvdW5kKHRhYiw0KQpgYGAKCiMgTXVsdGlwbGUgbGluZWFyIHJlZ3Jlc3Npb24gdmlhIEdpYmJzIHNhbXBsZXIKClN1cHBvc2UgdGhlIGRhdGEgaXMgJFxtYm94e2RhdGF9PVx7KHlfMSx4XzEpLFxsZG90cywoeV9uLHhfbilcfSQsIHdoZXJlICR5X2kkIGlzIGEgcmVhbCBudW1iZXIgYW5kICR4X2k9KHhfe2kxfSxcZG90cyx4X3tpcH0pJyQgaXMgYSAkcCQtZGltZW5zaW9uYWwgdmVjdG9yIG9mIGV4cGxhbmF0b3J5IHZhcmlhYmxlcywgYWxzbyBrbm93biBhcyByZWdyZXNzb3JzIG9yIGZlYXR1cmVzIG9yIGNoYXJhY3RlcmlzdGljcyBvZiB1bml0ICRpJCwgZm9yICRpPTEsXGxkb3RzLG4kLiAgT25lIG9mIHRoZSBtb3N0IGNvbW1vbiBtb2RlbHMgZm9yIHRoaXMga2luZCBvZiBhcnJhbmdlbWVudCBpcyB0aGF0IG9mIGEgR2F1c3NpYW4gbGluZWFyIHJlZ3Jlc3Npb24sIHdoZXJlIAokJAp5X2kgfCB4X2ksXGJldGEsXHNpZ21hXjIgXHNpbSBOKHhfaSdcYmV0YSxcc2lnbWFeMiksCiQkCndoZXJlICRcYmV0YT0oXGJldGFfMSxcbGRvdHMsXGJldGFfcCknJCBhbmQgJHkkcyBjb25kaXRpb25hbGx5IGluZGVwZW5kZW50IGdpdmVuICR4X2ksXGJldGEsXHNpZ21hXjIkLiAgCgojIyBMaWtlbGlob29kIGZ1bmN0aW9uClRoaXMgaGVhdmlseSBzdHJ1Y3R1cmVkIG1vZGVsIChsaW5lYXJpdHksIGdhdXNzaWFuaXR5LCBleG9nZW5laXR5LCBob21va2VkYXN0aWNpdHkpIGhhcyBsaWtlbGlob29kIAokJApMKFxiZXRhXzEsXGxkb3RzLFxiZXRhX3AsXHNpZ21hXjIpID0gXHByb2Rfe2k9MX1ebiAoMlxwaVxzaWdtYV4yKV57LTEvMn1cZXhwXGxlZnRce1xmcmFjezF9ezJcc2lnbWEgXjJ9KHlfaS14X2knXGJldGEpXjJccmlnaHRcfSwKJCQKZm9yICRcYmV0YV9pIFxpbiBcUmUkIGFuZCAkXHNpZ21hXjIgXGluIFxSZV4rJC4gIEluIG1hdHJpeCBub3RhdGlvbiwgCiQkCnl8WCxcYmV0YSxcc2lnbWFeMiBcc2ltIE4oWFxiZXRhLFxzaWdtYV4ySV9uKSwKJCQKd2hlcmUgJHk9KHlfMSxcbGRvdHMseV9uKSQsICRYPSh4XzEsXGxkb3RzLHhfbiknJCBvZiBkaW1lbnNpb24gJChuIFx0aW1lcyBwKSQgYW5kICRJX24kIGlzIHRoZSBpZGVudGl0eSBtYXRyaXggb2Ygb3JkZXIgJG4kLiAgSW4gdGhpcyBjYXNlIHRoZSBsaWtlbGlob29kIGJlY29tZXMKJCQKTChcYmV0YSxcc2lnbWFeMikgXHByb3B0byAoXHNpZ21hXjIpXnstXGZyYWN7bn17Mn19XGV4cFxsZWZ0XHtcZnJhY3sxfXsyXHNpZ21hIF4yfSh5LVhcYmV0YSknKHktWFxiZXRhKVxyaWdodFx9LAokJApBIGZ1cnRoZXIgc2ltcGxpY2F0aW9uIGxlYWRzIHRvIAokJApMKFxiZXRhLFxzaWdtYV4yKSBccHJvcHRvIChcc2lnbWFeMileey1cZnJhY3tufXsyfX1cZXhwXGxlZnRce1xmcmFjezF9ezJcc2lnbWEgXjJ9KFxiZXRhJyBYJ1ggXGJldGEgLSAyXGJldGEnWCd5K3kneSlccmlnaHRcfS4KJCQKSW4gb3RoZXIgd29yZHMsIHRoZSBzdWZmaWNpZW50IHN0YXRpc3RpY3MgZm9yIHRoZSBHYXVzc2lhbiBsaW5lYXIgcmVncmVzc2lvbiBhcmUgdGhlIGNyb3NzLXByb2R1Y3RzCiRTX3t4eH09WCdYJCwgJFNfe3h5fT1YJ3kkIGFuZCAkU197eXl9PXkneSQuCgojIyBQcmlvciBkaXN0cmlidXRpb24KTGV0IHVzIGFzc3VtZSB0aGF0LCBmb3IgJGk9MSxcbGRvdHMscCQsICRcYmV0YV9pIFxzaW0gTigwLFx0YXVeMikkLCBzYXkgJFx0YXU9MSQgb3IgJFx0YXU9MTAkLCBhbmQgJHAoXHNpZ21hXjIpIFxwcm9wdG8gMS9cc2lnbWFeMiQuICBUaGlzIHdheSB0aGUgcHJpb3Igb2YgJChcYmV0YSxcc2lnbWFeMikkIGlzIAokJApwKFxiZXRhLFxzaWdtYV4yKSBccHJvcHRvIFxmcmFjezF9e1xzaWdtYV4yfVxleHBcbGVmdFx7LVxmcmFjezF9ezJcdGF1XjJ9XHN1bV97aj0xfV5wIFxiZXRhX2peMlxyaWdodFx9PQpcZnJhY3sxfXtcc2lnbWFeMn1cZXhwXGxlZnRcey1cZnJhY3sxfXsyXHRhdV4yfSBcYmV0YSdcYmV0YVxyaWdodFx9CiQkCgojIyBQb3N0ZXJpb3IgZGlzdHJpYnV0aW9uCgpDb21iaW5pbmcgdGhlIGxpa2VsaWhvb2QgZnVuY3Rpb24gd2l0aCB0aGUgcHJpb3IgZGVuc2l0eSwgd2UgYXJyaXZlIGF0IHRoZSBrZXJuZWwgb2YgdGhlIHBvc3RlcmlvciBkZW5zaXR5CiQkCnAoXGJldGEsXHNpZ21hXjJ8XG1ib3h7ZGF0YX0pIFxwcm9wdG8gClxsZWZ0Wyhcc2lnbWFeMileey0xfVxleHBcbGVmdFx7LVxmcmFjezF9ezJcdGF1XjJ9IFxiZXRhJ1xiZXRhXHJpZ2h0XH1ccmlnaHRdIFx0aW1lcyAKXGxlZnRbKFxzaWdtYV4yKV57LVxmcmFje259ezJ9fVxleHBcbGVmdFx7LVxmcmFjezF9ezJcc2lnbWEgXjJ9KFxiZXRhJyBYJ1ggXGJldGEgLSAyXGJldGEnWCd5K3kneSlccmlnaHRcfVxyaWdodF0uCiQkCgojIyBGdWxsIGNvbmRpdGlvbmFsIGRpc3RyaWJ1dGlvbnMKClN1cHBvc2Ugd2UgZG8gbm90IHJlY29nbml6ZSB0aGlzIGpvaW50IHBvc3RlcmlvciwgYnV0IGFyZSBhYmxlIHRvIGRlcml2ZSB0aGUgZnVsbCBjb25kaXRpb25hbHMgCiQkCnAoXGJldGF8XHNpZ21hXjIsXG1ib3h7ZGF0YX0pIFwgXCBcIFxtYm94e2FuZH0gXCBcIFwgcChcc2lnbWFeMnxcYmV0YSxcbWJveHtkYXRhfSkuCiQkCgojIyMgRnVsbCBjb25kaXRpb25hbCBvZiAkXHNpZ21hXjIkIAoKJCQKcChcc2lnbWFeMnxcYmV0YSxcbWJveHtkYXRhfSkgXHByb3B0byAoXHNpZ21hXjIpXnstXGxlZnQoXGZyYWN7bn17Mn0rMVxyaWdodCl9ClxleHBcbGVmdFx7LVxmcmFjezF9ezJcc2lnbWEgXjJ9XGxlZnQoXGJldGEnIFgnWCBcYmV0YSAtIDJcYmV0YSdYJ3kreSd5XHJpZ2h0KVxyaWdodFx9LAokJAp3aGljaCBsb29rcyBsaWtlIHRoZSBrZXJuZWwgb2YgYW4gZm9sbG93aW5nIGludmVyc2UtZ2FtbWEgZGlzdHJpYnV0aW9uCiQkCihcc2lnbWFeMnxcYmV0YSxcbWJveHtkYXRhfSkgXHNpbSBJR1xsZWZ0KFxmcmFje259ezJ9LFxmcmFje1xiZXRhJyBTX3t4eH0gXGJldGEgLSAyXGJldGEnU197eHl9K1Nfe3l5fSl9ezJ9XHJpZ2h0KSBcZXF1aXYgSUdcbGVmdChcZnJhY3tufXsyfSxcZnJhY3soeS1YXGJldGEpJyh5LVhcYmV0YSl9ezJ9XHJpZ2h0KS4KJCQKCiMjIyBGdWxsIGNvbmRpdGlvbmFsIG9mICRcYmV0YSQKCldlIGNhbiBzZWUgdGhhdAokJApwKFxiZXRhfFxzaWdtYV4yLFxtYm94e2RhdGF9KSBccHJvcHRvIFxleHBcbGVmdFx7LVxmcmFjezF9ezJcdGF1XjJ9IFxiZXRhJ1xiZXRhXHJpZ2h0XH0KXGV4cFxsZWZ0XHstXGZyYWN7MX17MlxzaWdtYSBeMn0oXGJldGEnIFgnWCBcYmV0YSAtIDJcYmV0YSdYJ3kpXHJpZ2h0XH0gPQpcZXhwXGxlZnRcey1cZnJhY3sxfXsyfVxsZWZ0KFxiZXRhJyhJX3AvXHRhdV4yK1gnWC9cc2lnbWFeMilcYmV0YS0yXGJldGEnWCd5L1xzaWdtYV4yXHJpZ2h0KVxyaWdodFx9LAokJAp3aGljaCBsb29rcyBsaWtlIHRoZSBrZXJuZWwgb2YgYSBtdWx0aXZhcmlhdGUgbm9ybWFsIGRpc3RyaWJ1dGlvbiB3aXRoIG1lYW4gYW5kIHZhcmlhbmNlCiQkCihJX3AvXHRhdV4yK1gnWC9cc2lnbWFeMileey0xfVgneS9cc2lnbWFeMiAgXCBcIFwgXG1ib3h7YW5kfSBcIFwgXCAoSV9wL1x0YXVeMitYJ1gvXHNpZ21hXjIpXnstMX0sCiQkCm9yCiQkCihcYmV0YXxcc2lnbWFeMixcbWJveHtkYXRhfSkgXHNpbSBOXGxlZnRbKElfcC9cdGF1XjIrU197eHh9L1xzaWdtYV4yKV57LTF9U197eHl9L3NpZ21hXjIsIChJX3AvXHRhdV4yK1Nfe3h4fS9cc2lnbWFeMileey0xfVxyaWdodF0uCiQkCgojIyBHaWJicyBzYW1wbGVyCgoxKSBTYW1wbGUgJFxiZXRhJCBmcm9tICROXGxlZnRbKElfcC9cdGF1XjIrU197eHh9L1xzaWdtYV4yKV57LTF9U197eHl9L1xzaWdtYV4yLCAoSV9wL1x0YXVeMitTX3t4eH0vXHNpZ21hXjIpXnstMX1ccmlnaHRdJAoKMikgU2FtcGxlICRcc2lnbWFeMiQgZnJvbSAkSUdcbGVmdFtuLzIsKFxiZXRhJyBTX3t4eH0gXGJldGEgLSAyXGJldGEnU197eHl9K1Nfe3l5fSkvMlxyaWdodF0kCgoKCiMgUmV0dXJuIG9uIGVkdWNhdGlvbgoKVGhlIGRhdGEgaXMgYSAxOTc2IFBhbmVsIFN0dWR5IG9mIEluY29tZSBEeW5hbWljcywgYmFzZWQgb24gZGF0YSBmb3IgdGhlIHByZXZpb3VzIHllYXIsIDE5NzUuICBPZiB0aGUgNzUzIG9ic2VydmF0aW9ucywgdGhlIGZpcnN0IDQyOCBhcmUgZm9yIHdvbWVuIHdpdGggcG9zaXRpdmUgaG91cnMgd29ya2VkIGluIDE5NzUsIHdoaWxlIHRoZSByZW1haW5pbmcgMzI1IG9ic2VydmF0aW9ucyBhcmUgZm9yIHdvbWVuIHdobyBkaWQgbm90IHdvcmsgZm9yIHBheSBpbiAxOTc1LiAgQSBtb3JlIGNvbXBsZXRlIGRpc2N1c3Npb24gb2YgdGhlIGRhdGEgaXMgZm91bmQgaW4gTXJveiBbMTk4N10sIEFwcGVuZGl4IDEuICAgVGhvbWFzIEEuIE1yb3ogKDE5ODcpIFRoZSBTZW5zaXRpdml0eSBvZiBhbiBFbXBpcmljYWwgTW9kZWwgb2YgTWFycmllZCBXb21lbidzIEhvdXJzIG9mIFdvcmsgdG8gRWNvbm9taWMgYW5kIFN0YXRpc3RpY2FsIEFzc3VtcHRpb25zLiAgRWNvbm9tZXRyaWNhLCBWb2wuIDU1LCBOby4gNCAoSnVseSAxOTg3KSwgcHAuIDc2NS03OTkuICBTdGFibGUgVVJMOiBodHRwOi8vd3d3LmpzdG9yLm9yZy9zdGFibGUvMTkxMTAyOS4gIFRoZSB2YXJpYWJsZXMgaW4gdGhlIGRhdGFzZXQgYXJlIGFzIGZvbGxvd3M6CgoKKiBMRlAgICJBIGR1bW15IHZhcmlhYmxlID0gMSBpZiB3b21hbiB3b3JrZWQgaW4gMTk3NSwgZWxzZSAwIjsKKiBXSFJTICJXaWZlJ3MgaG91cnMgb2Ygd29yayBpbiAxOTc1IjsKKiBLTDYgICJOdW1iZXIgb2YgY2hpbGRyZW4gbGVzcyB0aGFuIDYgeWVhcnMgb2xkIGluIGhvdXNlaG9sZCI7CiogSzYxOCAiTnVtYmVyIG9mIGNoaWxkcmVuIGJldHdlZW4gYWdlcyA2IGFuZCAxOCBpbiBob3VzZWhvbGQiOwoqIFdBICAgIldpZmUncyBhZ2UiOwoqIFdFICAgIldpZmUncyBlZHVjYXRpb25hbCBhdHRhaW5tZW50LCBpbiB5ZWFycyI7CiogV1cgICAiV2lmZSdzIGF2ZXJhZ2UgaG91cmx5IGVhcm5pbmdzLCBpbiAxOTc1IGRvbGxhcnMiOwoqIFJQV0cgIldpZmUncyB3YWdlIHJlcG9ydGVkIGF0IHRoZSB0aW1lIG9mIHRoZSAxOTc2IGludGVydmlldyAobm90IHRoZSBzYW1lIGFzIHRoZSAxOTc1IGVzdGltYXRlZCB3YWdlKS4gIFRvIHVzZSB0aGUgc3Vic2FtcGxlIHdpdGggdGhpcyB3YWdlLCBvbmUgbmVlZHMgdG8gc2VsZWN0IDE5NzUgd29ya2VycyB3aXRoIExGUD0xLCB0aGVuIHNlbGVjdCBvbmx5IHRob3NlIHdvbWVuIHdpdGggbm9uLXplcm8gUlBXRy4gIE9ubHkgMzI1IHdvbWVuIHdvcmsgaW4gMTk3NSBhbmQgaGF2ZSBhIG5vbi16ZXJvIFJQV0cgaW4gMTk3Ni4iOwoqIEhIUlMgIkh1c2JhbmQncyBob3VycyB3b3JrZWQgaW4gMTk3NSI7CiogSEEgICAiSHVzYmFuZCdzIGFnZSI7CiogSEUgICAiSHVzYmFuZCdzIGVkdWNhdGlvbmFsIGF0dGFpbm1lbnQsIGluIHllYXJzIjsKKiBIVyAgICJIdXNiYW5kJ3Mgd2FnZSwgaW4gMTk3NSBkb2xsYXJzIjsKKiBGQU1JTkMgIkZhbWlseSBpbmNvbWUsIGluIDE5NzUgZG9sbGFycy4gIFRoaXMgdmFyaWFibGUgaXMgdXNlZCB0byBjb25zdHJ1Y3QgdGhlIHByb3BlcnR5IGluY29tZSB2YXJpYWJsZS4iOwoqIE1UUiAgIlRoaXMgaXMgdGhlIG1hcmdpbmFsIHRheCByYXRlIGZhY2luZyB0aGUgd2lmZSwgYW5kIGlzIHRha2VuIGZyb20gcHVibGlzaGVkIGZlZGVyYWwgdGF4IHRhYmxlcyAoc3RhdGUgYW5kIGxvY2FsIGluY29tZSB0YXhlcyBhcmUgZXhjbHVkZWQpLiBUaGUgdGF4YWJsZSBpbmNvbWUgb24gd2hpY2ggdGhpcyB0YXggcmF0ZSBpcyBjYWxjdWxhdGVkIGluY2x1ZGVzIFNvY2lhbCBTZWN1cml0eSwgaWYgYXBwbGljYWJsZSB0byB3aWZlLiI7CiogV01FRCAiV2lmZSdzIG1vdGhlcidzIGVkdWNhdGlvbmFsIGF0dGFpbm1lbnQsIGluIHllYXJzIjsKKiBXRkVEICJXaWZlJ3MgZmF0aGVyJ3MgZWR1Y2F0aW9uYWwgYXR0YWlubWVudCwgaW4geWVhcnMiOwoqIFVOICAgIlVuZW1wbG95bWVudCByYXRlIGluIGNvdW50eSBvZiByZXNpZGVuY2UsIGluIHBlcmNlbnRhZ2UgcG9pbnRzLiAgVGhpcyB0YWtlbiBmcm9tIGJyYWNrZXRlZCByYW5nZXMuIjsKKiBDSVQgICJEdW1teSB2YXJpYWJsZSA9IDEgaWYgbGl2ZSBpbiBsYXJnZSBjaXR5IChTTVNBKSwgZWxzZSAwIjsKKiBBWCAgICJBY3R1YWwgeWVhcnMgb2Ygd2lmZSdzIHByZXZpb3VzIGxhYm9yIG1hcmtldCBleHBlcmllbmNlIjsKCmBgYHtyIGZpZy53aWR0aD0xMCwgZmlnLmhlaWdodD02LCBmaWcuYWxpZ249J2NlbnRlcid9CmRhdGEgPSByZWFkLnRhYmxlKCJodHRwOi8vaGVkaWJlcnQub3JnL3dwLWNvbnRlbnQvdXBsb2Fkcy8yMDIwLzAxL21yb3otZGF0YS50eHQiLGhlYWRlcj1UUlVFKQoKYXR0YWNoKGRhdGEpCgpuYW1lcz0gYygid29ya2luZyAoVykiLAogICAgICAgICAiaG91cnMgd29yayAoVykiLAogICAgICAgICAiY2hpbGRyZW4gNi0xOCIsCiAgICAgICAgICJob3VybHkgZWFybmluZ3MgKFcpIiwKICAgICAgICAgImhvdXJzIHdvcmtlZCAoSCkiLAogICAgICAgICJ3YWdlIChIKSIsCiAgICAgICAgICJtYXJnaW5hbCB0YXggcmF0ZSAoVykiLAogICAgICAgICAibGl2ZSBpbiBsYXJnZSBjaXR5IikKCnkgPSBzY2FsZShsb2coRkFNSU5DKSkKWCA9IGNiaW5kKExGUCxXSFJTLEs2MTgsV1csSEhSUyxIVyxNVFIsQ0lUKQpuID0gbnJvdyhYKQpwID0gbmNvbChYKQoKCnBhcihtZnJvdz1jKDIsNCkpCmJveHBsb3QoeX5YWywxXSx4bGFiPW5hbWVzWzFdLG1haW49IiIsb3V0bGluZT1GQUxTRSxuYW1lcz1jKCJObyIsIlllcyIpLHlsYWI9IkZhbWlseSBpbmNvbWUgKGxvZykiKQpmb3IgKGkgaW4gMjoocC0xKSkKICBwbG90KFhbLGldLHkseGxhYj1uYW1lc1tpXSx5bGFiPSJGYW1pbHkgaW5jb21lIChsb2cpIikKYm94cGxvdCh5flhbLHBdLHhsYWI9bmFtZXNbcF0sbWFpbj0iIixvdXRsaW5lPUZBTFNFLG5hbWVzPWMoIk5vIiwiWWVzIikseWxhYj0iRmFtaWx5IGluY29tZSAobG9nKSIpCmBgYAoKIyMgT3JkaW5hcnkgbGVhc3Qgc3F1YXJlcwpgYGB7cn0KZml0Lm9scyA9IGxtKHl+WC0xKQpzdW1tYXJ5KGZpdC5vbHMpCmJldGEub2xzID0gZml0Lm9scyRjb2VmCnNpZy5vbHMgID0gc3FydChzdW0oKHktWCUqJWJldGEub2xzKV4yKS8obi1wKSkKYGBgCgojIyBCYXllc2lhbiBwb3N0ZXJpb3IgaW5mZXJlbmNlIHZpYSBHaWJicyBzYW1wbGVyCmBgYHtyfQpsaWJyYXJ5KCJtdnRub3JtIikKdGF1MiAgID0gMS4wClN4eCAgICA9IHQoWCklKiVYClN4eSAgICA9IHQoWCklKiV5CmEgICAgICA9IG4vMgpJcCAgICAgPSBkaWFnKDEvdGF1MixwKQpzaWdtYTIgPSAxCk0wICAgICA9IDEwMDAwCk0gICAgICA9IDEwMDAwCm5pdGVyICA9IE0wK00KZHJhd3MgID0gbWF0cml4KDAsbml0ZXIscCsxKQpmb3IgKGl0ZXIgaW4gMTpuaXRlcil7CiAgViAgICAgID0gc29sdmUoSXArU3h4L3NpZ21hMikKICBtICAgICAgPSBWJSolU3h5L3NpZ21hMgogIGJldGEgICA9IG1hdHJpeChybXZub3JtKDEsbWVhbj1tLHNpZ21hPVYpLHAsMSkKICBzaWdtYTIgPSAxL3JnYW1tYSgxLGEsc3VtKCh5LVglKiViZXRhKV4yKS8yKQogIGRyYXdzW2l0ZXIsXSA9IGMoYmV0YSxzcXJ0KHNpZ21hMikpCn0KZHJhd3MgPSBkcmF3c1soTTArMSk6bml0ZXIsXQpgYGAKCiMjIEdyYXBoaWNhbCBzdW1tYXJpZXMKCmBgYHtyIGZpZy53aWR0aD0xMCwgZmlnLmhlaWdodD03LCBmaWcuYWxpZ249J2NlbnRlcid9CnBhcihtZnJvdz1jKDMsMykpCmZvciAoaSBpbiAxOnApewogIGhpc3QoZHJhd3NbLGldLHByb2I9VFJVRSxtYWluPW5hbWVzW2ldLHhsYWI9cGFzdGUoImJldGEoIixpLCIpIixzZXA9IiIpKQogIGFibGluZSh2PTAsY29sPTIsbHdkPTMpCiAgYWJsaW5lKHY9YmV0YS5vbHNbaV0sY29sPTMsbHdkPTMpCn0KaGlzdChkcmF3c1sscCsxXSxwcm9iPVRSVUUsbWFpbj0ic2lnbWEiLHhsYWI9IiIpCmFibGluZSh2PXNpZy5vbHMsY29sPTMsbHdkPTMpCmBgYAo=