Jennifer Ganeles
2/24/19

Comparing MLE Results: Slide 4.15 vs. Slide 4.19

Using maximum likelihood estimation, both Slide 4.15 and Slide 4.19 show the relationship between education (independent variable) and income (dependent variable) using the “turnout” dataset from the Zelig package. However, Slide 4.15 maximizes the likelihood estimate of mean income (mu), while Slide 4.19 estimates the standard deviation (sigma) of income (so estimated parameters are different). In other words, Slide 4.15 looks at the effect of education on average income (the more common appraoch) while Slide 4.19 looks at the effect of education on the variation of income, or income inequality (less commonly used, but nevertheless useful in this example).

Relationship Between Education and Mean Income (Slide 4.15)

The results on Slide 4.15 use a linear regression model of Y∼N(β0+β1X,σ), in which beta1 represents the y-intercept for mu and beta2 represents the slope for mu. As shown below, a beta1 (y-intercept) of -.65 means that those who have no education (x=0) have an estimated mean income of -.65. A beta2 (slope) of .38 means that an increase in education by one year will have an estimated increase in mean income by .38 units. This result is statistically significant, and suggests a positive correlation between education and income, where those with higher education will also have a higher mean income.

library(Zelig)
data("turnout")
#Looking at the relationship between education and mean income
ols.lf<-function(param) {
  beta<-param[-1]
  sigma<-param[1]
  y<-as.vector(turnout$income)
  x<-cbind(1, turnout$educate)
  mu<-x%*%beta
  sum(dnorm(y, mu, sigma, log=TRUE))
}
library(maxLik)
mle_ols<-maxLik(logLik=ols.lf, start=c(sigma=1, beta1=1, beta2=1))
summary(mle_ols)
--------------------------------------------
Maximum Likelihood estimation
Newton-Raphson maximisation, 12 iterations
Return code 2: successive function values within tolerance limit
Log-Likelihood: -4691.256 
3  free parameters
Estimates:
      Estimate Std. error t value Pr(> t)    
sigma  2.52613    0.03989  63.326 < 2e-16 ***
beta1 -0.65207    0.20827  -3.131 0.00174 ** 
beta2  0.37613    0.01663  22.612 < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
--------------------------------------------

Checking Results:

summary(lm(income~educate,data=turnout))

Call:
lm(formula = income ~ educate, data = turnout)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.2028 -1.7363 -0.4273  1.3150 11.0632 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.65207    0.21016  -3.103  0.00194 ** 
educate      0.37613    0.01677  22.422  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.527 on 1998 degrees of freedom
Multiple R-squared:  0.201, Adjusted R-squared:  0.2006 
F-statistic: 502.8 on 1 and 1998 DF,  p-value: < 2.2e-16

Relationship Between Education and Income Inequality (Slide 4.19)

The results on Slide 4.19 use a linear regression model of Y∼N(μ,θ0+θ1X, in which theta1 represents the y-intercept for sigma and theta2 represents the slope for sigma. A theta1 (y-intercept) of 1.46 indicates that those with no education (x=0) have an estimated standard deviation of income of 1.46. A theta2 (slope) of .109 means that an increase in education by one year will have an estimated increase in the standard deviation of income by .109 units. These results are also statistically significant and suggest a positive correlation between education and income inequality, in which those with higher education will also see more variation in income.

#Looking at the relationship between education and standard deviation of income (i.e. income inequality)
ols.lf2<-function(param){
  mu<-param[1]
  theta<-param[-1]
  y<-as.vector(turnout$income)
  x<-cbind(1, turnout$educate)
  sigma<-x%*%theta
  sum(dnorm(y,mu,sigma,log=TRUE))
}
library(maxLik)
mle_ols2<-maxLik(logLik=ols.lf2, start=c(mu=1, theta1=1, theta2=1))
summary(mle_ols2)
--------------------------------------------
Maximum Likelihood estimation
Newton-Raphson maximisation, 9 iterations
Return code 2: successive function values within tolerance limit
Log-Likelihood: -4861.964 
3  free parameters
Estimates:
       Estimate Std. error t value Pr(> t)    
mu     3.516764   0.070320   50.01  <2e-16 ***
theta1 1.461011   0.106745   13.69  <2e-16 ***
theta2 0.109081   0.009185   11.88  <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
--------------------------------------------

Introducing Age as a Second Independent Variable for Both Models

If we introduce age as another independent variable, it can be hypothesized that age and mean income will be positively correlated since people become more educated as they get older (and as we just saw, higher education translates to higher mean income). However, given that people eventually retire and earn less as they grow older, it is also possible that age negatively correlates with mean income after a certain age is reached. The following results show an intercept of -.446, which refers to a mean income of -.446 units when education and age equal zero. The results also suggest a negative correlatation between age and mean income (beta3=-.003), implying that an increase in age by one year will decrease mean income by .003 units. However, it is important to note that this relationship is insignificant. Meanwhile, the slope for education is similar to our previous results, showing a statistically significant increase in mean income by .37 (beta2=.37) as education increases by one year. Education, therefore, seems to be a better predictor for mean income than age.

When it comes to age and income inequality, I would hypothesize that age and the standard deviation of income would be positively correlated. Indeed, the results show an intercept of .36 and a positive correlation between age and variation of income (theta3=.018). The slope for age implies that as age increases by one year, the standard deviation of income increases by .018 units. This relationship is statistically significant and suggests that as age increases, so does income inequality. These results also show a statistically significant positive correlation between education and variation of income, with a slightly higher slope for education than our previous result exhibited (theta2=.133). As one can see, income inequality positively correlates with both increases in education and increases in age.

Adding Age to Model 1:

#Looking at the relationship between age and education and mean income
ols.lf3<-function(param) {
  beta<-param[-1]
  sigma<-param[1]
  y<-as.vector(turnout$income)
  x<-cbind(1, turnout$educate, turnout$age)
  mu<-x%*%beta
  sum(dnorm(y, mu, sigma, log=TRUE))
}
mle_ols3<-maxLik(logLik=ols.lf3, start=c(sigma=1, beta1=1, beta2=1, beta3=1))
summary(mle_ols3)
--------------------------------------------
Maximum Likelihood estimation
Newton-Raphson maximisation, 16 iterations
Return code 2: successive function values within tolerance limit
Log-Likelihood: -4690.815 
4  free parameters
Estimates:
       Estimate Std. error t value Pr(> t)    
sigma  2.525576   0.039919  63.268  <2e-16 ***
beta1 -0.446047   0.300583  -1.484   0.138    
beta2  0.371011   0.017493  21.209  <2e-16 ***
beta3 -0.003184   0.003373  -0.944   0.345    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
--------------------------------------------

Checking Results:

#Checking Results
summary(lm(income~educate+age,data=turnout))

Call:
lm(formula = income ~ educate + age, data = turnout)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.2128 -1.7471 -0.4217  1.3042 11.1256 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.446084   0.303955  -1.468    0.142    
educate      0.371013   0.017641  21.031   <2e-16 ***
age         -0.003183   0.003394  -0.938    0.348    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.527 on 1997 degrees of freedom
Multiple R-squared:  0.2014,    Adjusted R-squared:  0.2006 
F-statistic: 251.8 on 2 and 1997 DF,  p-value: < 2.2e-16

Adding Age to Model 2:

#Looking at the relationship between age and educaton and standard deviation of income (income inequality)
ols.lf4<-function(param){
  mu<-param[1]
  theta<-param[-1]
  y<-as.vector(turnout$income)
  x<-cbind(1, turnout$educate, turnout$age)
  sigma<-x%*%theta
  sum(dnorm(y, mu, sigma, log=TRUE))
}
mle_ols4<-maxLik(logLik=ols.lf4, start=c(mu=1, theta1=1, theta2=1, theta3=1), method="bfgs")
summary(mle_ols4)
--------------------------------------------
Maximum Likelihood estimation
BFGS maximization, 150 iterations
Return code 0: successful convergence 
Log-Likelihood: -4843.15 
4  free parameters
Estimates:
       Estimate Std. error t value  Pr(> t)    
mu     3.555011   0.069193  51.378  < 2e-16 ***
theta1 0.362114   0.204550   1.770   0.0767 .  
theta2 0.133349   0.010756  12.398  < 2e-16 ***
theta3 0.017507   0.002852   6.139 8.32e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
--------------------------------------------

Plotting Age and Income Variables

As hypothesized above, conditional mean income appears to increase with age until about age 50, where it then begins to decline.

library(ggplot2)
ggplot(turnout)+
  geom_point(aes(x = age, y = income)) + geom_smooth(aes(x = age, y = income))

LS0tCnRpdGxlOiAiU09DIDcxMjogSG9tZXdvcmsgMyIKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQoqSmVubmlmZXIgR2FuZWxlcyoKPGJyLz4qMi8yNC8xOSoKCiMjQ29tcGFyaW5nIE1MRSBSZXN1bHRzOiBTbGlkZSA0LjE1IHZzLiBTbGlkZSA0LjE5CgpVc2luZyBtYXhpbXVtIGxpa2VsaWhvb2QgZXN0aW1hdGlvbiwgYm90aCBTbGlkZSA0LjE1IGFuZCBTbGlkZSA0LjE5IHNob3cgdGhlIHJlbGF0aW9uc2hpcCBiZXR3ZWVuIGVkdWNhdGlvbiAoaW5kZXBlbmRlbnQgdmFyaWFibGUpIGFuZCBpbmNvbWUgKGRlcGVuZGVudCB2YXJpYWJsZSkgdXNpbmcgdGhlICJ0dXJub3V0IiBkYXRhc2V0IGZyb20gdGhlICpaZWxpZyogcGFja2FnZS4gSG93ZXZlciwgU2xpZGUgNC4xNSBtYXhpbWl6ZXMgdGhlIGxpa2VsaWhvb2QgZXN0aW1hdGUgb2YgKiptZWFuIGluY29tZSAobXUpKiosIHdoaWxlIFNsaWRlIDQuMTkgZXN0aW1hdGVzIHRoZSAqKnN0YW5kYXJkIGRldmlhdGlvbiAoc2lnbWEpIG9mIGluY29tZSoqIChzbyBlc3RpbWF0ZWQgcGFyYW1ldGVycyBhcmUgZGlmZmVyZW50KS4gSW4gb3RoZXIgd29yZHMsIFNsaWRlIDQuMTUgbG9va3MgYXQgdGhlIGVmZmVjdCBvZiBlZHVjYXRpb24gb24gYXZlcmFnZSBpbmNvbWUgKHRoZSBtb3JlIGNvbW1vbiBhcHByYW9jaCkgd2hpbGUgU2xpZGUgNC4xOSBsb29rcyBhdCB0aGUgZWZmZWN0IG9mIGVkdWNhdGlvbiBvbiB0aGUgdmFyaWF0aW9uIG9mIGluY29tZSwgb3IgaW5jb21lIGluZXF1YWxpdHkgKGxlc3MgY29tbW9ubHkgdXNlZCwgYnV0IG5ldmVydGhlbGVzcyB1c2VmdWwgaW4gdGhpcyBleGFtcGxlKS4gCgojIyNSZWxhdGlvbnNoaXAgQmV0d2VlbiBFZHVjYXRpb24gYW5kIE1lYW4gSW5jb21lIChTbGlkZSA0LjE1KQpUaGUgcmVzdWx0cyBvbiBTbGlkZSA0LjE1IHVzZSBhIGxpbmVhciByZWdyZXNzaW9uIG1vZGVsIG9mIFniiLxOKM6yMCvOsjFYLM+DKSwgaW4gd2hpY2ggYmV0YTEgcmVwcmVzZW50cyB0aGUgeS1pbnRlcmNlcHQgZm9yIG11IGFuZCBiZXRhMiByZXByZXNlbnRzIHRoZSBzbG9wZSBmb3IgbXUuIEFzIHNob3duIGJlbG93LCBhIGJldGExICh5LWludGVyY2VwdCkgb2YgLS42NSBtZWFucyB0aGF0IHRob3NlIHdobyBoYXZlIG5vIGVkdWNhdGlvbiAoeD0wKSBoYXZlIGFuIGVzdGltYXRlZCBtZWFuIGluY29tZSBvZiAtLjY1LiBBIGJldGEyIChzbG9wZSkgb2YgLjM4IG1lYW5zIHRoYXQgYW4gaW5jcmVhc2UgaW4gZWR1Y2F0aW9uIGJ5IG9uZSB5ZWFyIHdpbGwgaGF2ZSBhbiBlc3RpbWF0ZWQgaW5jcmVhc2UgaW4gbWVhbiBpbmNvbWUgYnkgLjM4IHVuaXRzLiBUaGlzIHJlc3VsdCBpcyBzdGF0aXN0aWNhbGx5IHNpZ25pZmljYW50LCBhbmQgc3VnZ2VzdHMgYSBwb3NpdGl2ZSBjb3JyZWxhdGlvbiBiZXR3ZWVuIGVkdWNhdGlvbiBhbmQgaW5jb21lLCB3aGVyZSB0aG9zZSB3aXRoIGhpZ2hlciBlZHVjYXRpb24gd2lsbCBhbHNvIGhhdmUgYSBoaWdoZXIgbWVhbiBpbmNvbWUuCmBgYHtyIHdhcm5pbmc9RkFMU0V9CgpsaWJyYXJ5KFplbGlnKQpkYXRhKCJ0dXJub3V0IikKCiNMb29raW5nIGF0IHRoZSByZWxhdGlvbnNoaXAgYmV0d2VlbiBlZHVjYXRpb24gYW5kIG1lYW4gaW5jb21lCgpvbHMubGY8LWZ1bmN0aW9uKHBhcmFtKSB7CiAgYmV0YTwtcGFyYW1bLTFdCiAgc2lnbWE8LXBhcmFtWzFdCiAgeTwtYXMudmVjdG9yKHR1cm5vdXQkaW5jb21lKQogIHg8LWNiaW5kKDEsIHR1cm5vdXQkZWR1Y2F0ZSkKICBtdTwteCUqJWJldGEKICBzdW0oZG5vcm0oeSwgbXUsIHNpZ21hLCBsb2c9VFJVRSkpCn0KCmxpYnJhcnkobWF4TGlrKQptbGVfb2xzPC1tYXhMaWsobG9nTGlrPW9scy5sZiwgc3RhcnQ9YyhzaWdtYT0xLCBiZXRhMT0xLCBiZXRhMj0xKSkKc3VtbWFyeShtbGVfb2xzKQpgYGAKKipDaGVja2luZyBSZXN1bHRzOioqCmBgYHtyfQpzdW1tYXJ5KGxtKGluY29tZX5lZHVjYXRlLGRhdGE9dHVybm91dCkpCmBgYAoKIyMjUmVsYXRpb25zaGlwIEJldHdlZW4gRWR1Y2F0aW9uIGFuZCBJbmNvbWUgSW5lcXVhbGl0eSAoU2xpZGUgNC4xOSkKVGhlIHJlc3VsdHMgb24gU2xpZGUgNC4xOSB1c2UgYSBsaW5lYXIgcmVncmVzc2lvbiBtb2RlbCBvZiBZ4oi8TijOvCzOuDArzrgxWCwgaW4gd2hpY2ggdGhldGExIHJlcHJlc2VudHMgdGhlIHktaW50ZXJjZXB0IGZvciBzaWdtYSBhbmQgdGhldGEyIHJlcHJlc2VudHMgdGhlIHNsb3BlIGZvciBzaWdtYS4gQSB0aGV0YTEgKHktaW50ZXJjZXB0KSBvZiAxLjQ2IGluZGljYXRlcyB0aGF0IHRob3NlIHdpdGggbm8gZWR1Y2F0aW9uICh4PTApIGhhdmUgYW4gZXN0aW1hdGVkIHN0YW5kYXJkIGRldmlhdGlvbiBvZiBpbmNvbWUgb2YgMS40Ni4gQSB0aGV0YTIgKHNsb3BlKSBvZiAuMTA5IG1lYW5zIHRoYXQgYW4gaW5jcmVhc2UgaW4gZWR1Y2F0aW9uIGJ5IG9uZSB5ZWFyIHdpbGwgaGF2ZSBhbiBlc3RpbWF0ZWQgaW5jcmVhc2UgaW4gdGhlIHN0YW5kYXJkIGRldmlhdGlvbiBvZiBpbmNvbWUgYnkgLjEwOSB1bml0cy4gVGhlc2UgcmVzdWx0cyBhcmUgYWxzbyBzdGF0aXN0aWNhbGx5IHNpZ25pZmljYW50IGFuZCBzdWdnZXN0IGEgcG9zaXRpdmUgY29ycmVsYXRpb24gYmV0d2VlbiBlZHVjYXRpb24gYW5kIGluY29tZSBpbmVxdWFsaXR5LCBpbiB3aGljaCB0aG9zZSB3aXRoIGhpZ2hlciBlZHVjYXRpb24gd2lsbCBhbHNvIHNlZSBtb3JlIHZhcmlhdGlvbiBpbiBpbmNvbWUuCmBgYHtyIHdhcm5pbmc9RkFMU0V9CiNMb29raW5nIGF0IHRoZSByZWxhdGlvbnNoaXAgYmV0d2VlbiBlZHVjYXRpb24gYW5kIHN0YW5kYXJkIGRldmlhdGlvbiBvZiBpbmNvbWUgKGkuZS4gaW5jb21lIGluZXF1YWxpdHkpCgpvbHMubGYyPC1mdW5jdGlvbihwYXJhbSl7CiAgbXU8LXBhcmFtWzFdCiAgdGhldGE8LXBhcmFtWy0xXQogIHk8LWFzLnZlY3Rvcih0dXJub3V0JGluY29tZSkKICB4PC1jYmluZCgxLCB0dXJub3V0JGVkdWNhdGUpCiAgc2lnbWE8LXglKiV0aGV0YQogIHN1bShkbm9ybSh5LG11LHNpZ21hLGxvZz1UUlVFKSkKfQoKbGlicmFyeShtYXhMaWspCm1sZV9vbHMyPC1tYXhMaWsobG9nTGlrPW9scy5sZjIsIHN0YXJ0PWMobXU9MSwgdGhldGExPTEsIHRoZXRhMj0xKSkKc3VtbWFyeShtbGVfb2xzMikKYGBgCgojI0ludHJvZHVjaW5nIEFnZSBhcyBhIFNlY29uZCBJbmRlcGVuZGVudCBWYXJpYWJsZSBmb3IgQm90aCBNb2RlbHMKCklmIHdlIGludHJvZHVjZSBhZ2UgYXMgYW5vdGhlciBpbmRlcGVuZGVudCB2YXJpYWJsZSwgaXQgY2FuIGJlIGh5cG90aGVzaXplZCB0aGF0IGFnZSBhbmQgbWVhbiBpbmNvbWUgd2lsbCBiZSBwb3NpdGl2ZWx5IGNvcnJlbGF0ZWQgc2luY2UgcGVvcGxlIGJlY29tZSBtb3JlIGVkdWNhdGVkIGFzIHRoZXkgZ2V0IG9sZGVyIChhbmQgYXMgd2UganVzdCBzYXcsIGhpZ2hlciBlZHVjYXRpb24gdHJhbnNsYXRlcyB0byBoaWdoZXIgbWVhbiBpbmNvbWUpLiBIb3dldmVyLCBnaXZlbiB0aGF0IHBlb3BsZSBldmVudHVhbGx5IHJldGlyZSBhbmQgZWFybiBsZXNzIGFzIHRoZXkgZ3JvdyBvbGRlciwgaXQgaXMgYWxzbyBwb3NzaWJsZSB0aGF0IGFnZSBuZWdhdGl2ZWx5IGNvcnJlbGF0ZXMgd2l0aCBtZWFuIGluY29tZSBhZnRlciBhIGNlcnRhaW4gYWdlIGlzIHJlYWNoZWQuIFRoZSBmb2xsb3dpbmcgcmVzdWx0cyBzaG93IGFuIGludGVyY2VwdCBvZiAtLjQ0Niwgd2hpY2ggcmVmZXJzIHRvIGEgbWVhbiBpbmNvbWUgb2YgLS40NDYgdW5pdHMgd2hlbiBlZHVjYXRpb24gYW5kIGFnZSBlcXVhbCB6ZXJvLiBUaGUgcmVzdWx0cyBhbHNvIHN1Z2dlc3QgYSBuZWdhdGl2ZSBjb3JyZWxhdGF0aW9uIGJldHdlZW4gYWdlIGFuZCBtZWFuIGluY29tZSAoYmV0YTM9LS4wMDMpLCBpbXBseWluZyB0aGF0IGFuIGluY3JlYXNlIGluIGFnZSBieSBvbmUgeWVhciB3aWxsIGRlY3JlYXNlIG1lYW4gaW5jb21lIGJ5IC4wMDMgdW5pdHMuIEhvd2V2ZXIsIGl0IGlzIGltcG9ydGFudCB0byBub3RlIHRoYXQgdGhpcyByZWxhdGlvbnNoaXAgaXMgaW5zaWduaWZpY2FudC4gTWVhbndoaWxlLCB0aGUgc2xvcGUgZm9yIGVkdWNhdGlvbiBpcyBzaW1pbGFyIHRvIG91ciBwcmV2aW91cyByZXN1bHRzLCBzaG93aW5nIGEgc3RhdGlzdGljYWxseSBzaWduaWZpY2FudCBpbmNyZWFzZSBpbiBtZWFuIGluY29tZSBieSAuMzcgKGJldGEyPS4zNykgYXMgZWR1Y2F0aW9uIGluY3JlYXNlcyBieSBvbmUgeWVhci4gRWR1Y2F0aW9uLCB0aGVyZWZvcmUsIHNlZW1zIHRvIGJlIGEgYmV0dGVyIHByZWRpY3RvciBmb3IgbWVhbiBpbmNvbWUgdGhhbiBhZ2UuCgpXaGVuIGl0IGNvbWVzIHRvIGFnZSBhbmQgaW5jb21lIGluZXF1YWxpdHksIEkgd291bGQgaHlwb3RoZXNpemUgdGhhdCBhZ2UgYW5kIHRoZSBzdGFuZGFyZCBkZXZpYXRpb24gb2YgaW5jb21lIHdvdWxkIGJlIHBvc2l0aXZlbHkgY29ycmVsYXRlZC4gSW5kZWVkLCB0aGUgcmVzdWx0cyBzaG93IGFuIGludGVyY2VwdCBvZiAuMzYgYW5kIGEgcG9zaXRpdmUgY29ycmVsYXRpb24gYmV0d2VlbiBhZ2UgYW5kIHZhcmlhdGlvbiBvZiBpbmNvbWUgKHRoZXRhMz0uMDE4KS4gVGhlIHNsb3BlIGZvciBhZ2UgaW1wbGllcyB0aGF0IGFzIGFnZSBpbmNyZWFzZXMgYnkgb25lIHllYXIsIHRoZSBzdGFuZGFyZCBkZXZpYXRpb24gb2YgaW5jb21lIGluY3JlYXNlcyBieSAuMDE4IHVuaXRzLiBUaGlzIHJlbGF0aW9uc2hpcCBpcyBzdGF0aXN0aWNhbGx5IHNpZ25pZmljYW50IGFuZCBzdWdnZXN0cyB0aGF0IGFzIGFnZSBpbmNyZWFzZXMsIHNvIGRvZXMgaW5jb21lIGluZXF1YWxpdHkuIFRoZXNlIHJlc3VsdHMgYWxzbyBzaG93IGEgc3RhdGlzdGljYWxseSBzaWduaWZpY2FudCBwb3NpdGl2ZSBjb3JyZWxhdGlvbiBiZXR3ZWVuIGVkdWNhdGlvbiBhbmQgdmFyaWF0aW9uIG9mIGluY29tZSwgd2l0aCBhIHNsaWdodGx5IGhpZ2hlciBzbG9wZSBmb3IgZWR1Y2F0aW9uIHRoYW4gb3VyIHByZXZpb3VzIHJlc3VsdCBleGhpYml0ZWQgKHRoZXRhMj0uMTMzKS4gQXMgb25lIGNhbiBzZWUsIGluY29tZSBpbmVxdWFsaXR5IHBvc2l0aXZlbHkgY29ycmVsYXRlcyB3aXRoIGJvdGggaW5jcmVhc2VzIGluIGVkdWNhdGlvbiBhbmQgaW5jcmVhc2VzIGluIGFnZS4KCiMjI0FkZGluZyBBZ2UgdG8gTW9kZWwgMToKYGBge3Igd2FybmluZz1GQUxTRX0KI0xvb2tpbmcgYXQgdGhlIHJlbGF0aW9uc2hpcCBiZXR3ZWVuIGFnZSBhbmQgZWR1Y2F0aW9uIGFuZCBtZWFuIGluY29tZQoKb2xzLmxmMzwtZnVuY3Rpb24ocGFyYW0pIHsKICBiZXRhPC1wYXJhbVstMV0KICBzaWdtYTwtcGFyYW1bMV0KICB5PC1hcy52ZWN0b3IodHVybm91dCRpbmNvbWUpCiAgeDwtY2JpbmQoMSwgdHVybm91dCRlZHVjYXRlLCB0dXJub3V0JGFnZSkKICBtdTwteCUqJWJldGEKICBzdW0oZG5vcm0oeSwgbXUsIHNpZ21hLCBsb2c9VFJVRSkpCn0KCm1sZV9vbHMzPC1tYXhMaWsobG9nTGlrPW9scy5sZjMsIHN0YXJ0PWMoc2lnbWE9MSwgYmV0YTE9MSwgYmV0YTI9MSwgYmV0YTM9MSkpCnN1bW1hcnkobWxlX29sczMpCgpgYGAKKipDaGVja2luZyBSZXN1bHRzOioqCmBgYHtyfQpzdW1tYXJ5KGxtKGluY29tZX5lZHVjYXRlK2FnZSxkYXRhPXR1cm5vdXQpKQpgYGAKCiMjI0FkZGluZyBBZ2UgdG8gTW9kZWwgMjoKYGBge3Igd2FybmluZz1GQUxTRX0KI0xvb2tpbmcgYXQgdGhlIHJlbGF0aW9uc2hpcCBiZXR3ZWVuIGFnZSBhbmQgZWR1Y2F0b24gYW5kIHN0YW5kYXJkIGRldmlhdGlvbiBvZiBpbmNvbWUgKGluY29tZSBpbmVxdWFsaXR5KQoKb2xzLmxmNDwtZnVuY3Rpb24ocGFyYW0pewogIG11PC1wYXJhbVsxXQogIHRoZXRhPC1wYXJhbVstMV0KICB5PC1hcy52ZWN0b3IodHVybm91dCRpbmNvbWUpCiAgeDwtY2JpbmQoMSwgdHVybm91dCRlZHVjYXRlLCB0dXJub3V0JGFnZSkKICBzaWdtYTwteCUqJXRoZXRhCiAgc3VtKGRub3JtKHksIG11LCBzaWdtYSwgbG9nPVRSVUUpKQp9CgptbGVfb2xzNDwtbWF4TGlrKGxvZ0xpaz1vbHMubGY0LCBzdGFydD1jKG11PTEsIHRoZXRhMT0xLCB0aGV0YTI9MSwgdGhldGEzPTEpLCBtZXRob2Q9ImJmZ3MiKQpzdW1tYXJ5KG1sZV9vbHM0KQpgYGAKCgojIyNQbG90dGluZyBBZ2UgYW5kIEluY29tZSBWYXJpYWJsZXMKQXMgaHlwb3RoZXNpemVkIGFib3ZlLCBjb25kaXRpb25hbCBtZWFuIGluY29tZSBhcHBlYXJzIHRvIGluY3JlYXNlIHdpdGggYWdlIHVudGlsIGFib3V0IGFnZSA1MCwgd2hlcmUgaXQgdGhlbiBiZWdpbnMgdG8gZGVjbGluZS4KYGBge3IgbWVzc2FnZT1GQUxTRX0KbGlicmFyeShnZ3Bsb3QyKQpnZ3Bsb3QodHVybm91dCkrCiAgZ2VvbV9wb2ludChhZXMoeCA9IGFnZSwgeSA9IGluY29tZSkpICsgZ2VvbV9zbW9vdGgoYWVzKHggPSBhZ2UsIHkgPSBpbmNvbWUpKQpgYGAK