Bayesian data analysis project

Data and the problem

The chosen data is Government debt as percent of GDP (available at: https://www.google.com/publicdata/explore?ds=ds22a34krhq5p_#! ). The data contains the yearly debt progression of European countries out of which a subset of countries is analyzed here. The range of data was limited between years 2002-2016 for consistent data availability.

The problem addressed is the prediction of the debt progress based on history data using gaussian processes.

library(ggplot2)
library(reshape2)
library(ggthemes)
data <- read.csv("data.tsv",header=T,sep="\t",na.strings = ":",row.names = 1)
data2 <- data[,-c(1:7)]
l <- rownames(data2)
n <- data2$name
data2 <- as.data.frame(t(data2[,-1]))
colnames(data2) <- n
data2$myfactor <- factor(row.names(data2))
colnames(data2) <- l
data2 <- data2[,-c(33)]
data2$year <- c(2003:2016)
cnames <- colnames(data2)
df_melt <- melt(data2[, cnames], id="year")
gg <- ggplot(df_melt, aes(x=year))
gg + geom_line(aes(y=value, color=variable), size=1) + scale_color_discrete(name="Legend") +
  ylab("% of GDP") +
  ggtitle("General government gross debt % of GDP") +
  theme_igray()

Model description (separate model for each country)

We have developed two models here. We start with separate models for each countries. In the second part of this report we develop the model further and add hierarchical components as well.

The data is a time-series data by nature which makes Gaussian Processes a natural choice. The growing trend of the data is incorporated into the model as a linear component:

μ_{i} = μ_{0} + μ_{b} \cdot x_{i}

$\begin{equation} \mu_i = \mu_0 + \mu_b \cdot x_i \end{equation}$

Then we have the exponentiated quadratic kernel as a covariance:

K_{i, j} = α^{2} \exp (- \frac{1}{2 l^{2}} (x_{i} - x_{j})^{2})

$\begin{equation} K_{i,j} = \alpha^2 \exp \Big( - \frac{1}{2 l^2} (x_i - x_j)^2 \Big) \end{equation}$

And together these compose the full separate model (in addition to the random noise sigma):

\begin{aligned} f & \sim M u l t i N o r m a l (μ, K) \\ y_{i} & \sim N o r m a l (f_{i}, σ) \end{aligned}

$\begin{align} f & \sim MultiNormal(\mu, K) \\ y_i & \sim Normal(f_i, \sigma) \end{align}$

Prior choices

There are in total four initial prior choices in the model. We have listed our inital choices for the priors here, however, these will adjusted later due to some convergence issues. The initial prior choices are:

$\mu_0 \sim Half\_Normal(0, 100)$
$\mu_b \sim Normal(0, 5)$
$\alpha \sim Normal(0, 1)$
$length\_param \sim Weibull(10, 2)$

In addition, there are two parameters $eta$ and $sigma$ for which standard prior choices were used. $Eta$ is needed for the Cholesky decomposition to generate the covariance matrix, $sigma$ is used to incorporate the random, unexplainable variation in the model. Both use standard prior choices $Normal(0, 1)$ .

$\mu_0 , \mu_b$

These two parameters form the linear trend for the model. The intercept $mu\_0$ is restricted to positive and left vague, $mu\_0 \sim Half\_Normal(0, 100)$ , which covers the range of government debts in year 2002. The slope parameter is also fairly vague but possibly non-positive as well, $mu\_b \sim Normal(0, 5)$ .

$\alpha$

Alpha is the strength parameter for the Gaussian Process. The vague prior $alpha \sim Normal(0,1)$ was the initial choice and seems to provide sufficient strength for the variation.

length_param

Length_param controls the speed at which the correlation effect decays as distance between data points grows (lower value = faster decay). Since the amount of data points is low and the data behaves extremely irregularly, using a vague prior for the length parameter resulted in low values ( < 1.0) value that overfit the data. Using prior choice $length\_param \sim Weibull(10, 2)$ , the model puts less emphasis on the high frequency fluctuation and fit becomes smoother. We think this is reasonable, because the financial crisis at 2008 is impossible to explain based only on this dataset, so we want to prevent overfitting as much as possible.

Stan (separate model)

data {
  // Samples used for fitting
    int<lower=1> N_is;
    real x_is[N_is];
    vector[N_is] y_is;
    
    // Samples used for prediction
    int<lower=1> N_pred;
    real x_pred[N_pred];
}

transformed data {
  int<lower=1> N = N_is + N_pred;
  real x[N];
  vector[N] xv;
  real jitter = 1e-9;

  for (n in 1:N_is) x[n] = x_is[n];
  for (n in 1:N_pred) x[N_is + n] = x_pred[n]; // Join the in sample and prediction features.
  xv = to_vector(x);
}

parameters {
  real<lower=0> alpha;   // Strength parameter for GP
  real<lower=1> length_param;  // Lag parameter
  real<lower=0> sigma;  // Random, unexplained noise

  vector[N] eta;   // Gaussian 0,1 for generating the covariance from Cholesky decompose
  
  real<lower=0> mu_0;
  real mu_b;
}

transformed parameters {
  vector[N] f;
  {
    matrix[N, N] K = cov_exp_quad(x, alpha, length_param) + diag_matrix(rep_vector(jitter, N));
    matrix[N, N] L_K = cholesky_decompose(K);
    f = mu_0 + mu_b * xv + L_K * eta;
  }
}

model {
  alpha ~ normal(0, 1);
  //length_param ~ gamma(4,4); // Vague prior
  length_param ~ weibull(10, 3); // More informative prior
  sigma ~ normal(0, 1);
  eta ~ normal(0, 1);
  
  mu_0 ~ normal(0, 100);
  mu_b ~ normal(0, 5);
  
  y_is ~ normal(f[1:N_is], sigma);
}

generated quantities {
  vector[N] y_pred;
  vector[N_is] log_lik;
  for (n in 1:N) y_pred[n] = normal_rng(f[n], sigma);
  for (n in 1:N_is) log_lik[n] = normal_lpdf(y_is[n] | f[n], sigma);
}

R code (separate model)

Training of the initial separate model and plots for different countries.

library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
data = read.csv('data.csv', row.names=1)
data = data.frame(t(data))
data[,'year'] = seq(2002, 2016, by=1)
gp <- function(country_id,country_name) {
  cc = country_id
  y = data[,cc]
  x = 1:length(y)
  
  ss = 3
  y_is = y[1:(length(y) - ss)]
  x_is = x[1:(length(x) - ss)]
  N_is = length(x_is)
  
  y_os = tail(y, (ss + 1))
  x_os = tail(x, (ss + 1))
  
  x_pred = tail(x, ss)
  N_pred = length(x_pred)
  
  gp_model = 'debt_gp_linear.stan'
  
  fit_gp <- stan(file=gp_model,
                 data=list(x_is = x_is,
                           y_is = y_is,
                           N_is = N_is,
                           x_pred = x_pred,
                           N_pred),
                 iter=200,
                 chains=3)
  
  sample = extract(fit_gp)
  sample_preds = data.frame(sample$y_pred)
    
  sample_quantiles <- apply(sample_preds,2,function(x){quantile(x,c(0.05,0.5,0.95))})
  
  sample_data <- data.frame(low = sample_quantiles[1,],
                            med = sample_quantiles[2,],
                            high = sample_quantiles[3,],
                            x = x)
  training_data <- data.frame(x=x_is, obs=y_is)
  validation_data <- data.frame(x=x_os, obs=y_os)
  
  cols <- c("GP median"="black","GP 90% percentile"="lightblue","True debt"="red")
  ggplot(sample_data, aes(x=x, y=med)) +
    geom_ribbon(aes(x=x,ymin=low,ymax=high,colour="GP 90% percentile"),alpha=0.6,fill="lightblue") +
    geom_line(aes(x=x,y=med,colour="GP median")) +
    geom_point() +
    geom_line(data=training_data, aes(x=x, y=obs,colour="True debt")) +
    geom_line(data=validation_data, aes(x=x, y=obs,colour="True debt"), linetype="dashed") +
    xlab('year') + 
    ylab('Debt (% of GDP)') + 
    ggtitle(sprintf('Government debt, country = \'%s\'', country_name)) +
    scale_colour_manual(name="Colour",values=cols) +
    theme(plot.title = element_text(hjust = 0.5))
}
gp('FI','Finland')

gp('SE','Sweden')

gp('EL','Greece')

Convergence diagnostics

Analysis is done using a separate model generated for Finland.

cc = 'FI'
y = data[,cc]
x = 1:length(y)
ss = 3
y_is = y[1:(length(y) - ss)]
x_is = x[1:(length(x) - ss)]
N_is = length(x_is)
y_os = tail(y, (ss + 1))
x_os = tail(x, (ss + 1))
x_pred = tail(x, ss)
N_pred = length(x_pred)
gp_model = 'debt_gp_linear_w.stan'
fit_gp <- stan(file=gp_model,
               data=list(x_is = x_is,
                         y_is = y_is,
                         N_is = N_is,
                         x_pred = x_pred,
                         N_pred),
               iter=1000,
               chains=10,control=list(adapt_delta=0.99,max_treedepth=15))

$\hat{R}$ diagnostics

fit_gp

Inference for Stan model: debt_gp_linear_w.
10 chains, each with iter=1000; warmup=500; thin=1; 
post-warmup draws per chain=500, total post-warmup draws=5000.

               mean se_mean   sd   2.5%    25%    50%    75%  97.5% n_eff Rhat
alpha          3.19    0.01 0.54   2.18   2.82   3.17   3.55   4.33  2463 1.00
length_param   1.80    0.01 0.27   1.21   1.63   1.83   2.00   2.23   660 1.01
sigma          1.62    0.02 0.54   0.48   1.28   1.59   1.93   2.77   659 1.01
eta[1]         1.17    0.02 0.88  -0.55   0.58   1.16   1.76   2.91  2583 1.00
eta[2]         0.63    0.01 0.70  -0.81   0.20   0.66   1.10   1.97  3879 1.00
eta[3]        -0.55    0.01 0.83  -2.21  -1.11  -0.52   0.01   1.03  3463 1.00
eta[4]        -1.53    0.02 0.93  -3.36  -2.17  -1.51  -0.86   0.20  2826 1.00
eta[5]        -2.01    0.02 0.88  -3.71  -2.60  -2.02  -1.41  -0.32  2377 1.00
eta[6]        -1.20    0.04 1.14  -3.51  -1.97  -1.14  -0.41   0.93   961 1.01
eta[7]         0.27    0.01 0.93  -1.57  -0.36   0.29   0.89   2.08  5000 1.00
eta[8]         0.50    0.02 0.86  -1.25  -0.06   0.50   1.07   2.15  2587 1.00
eta[9]         0.21    0.02 0.98  -1.63  -0.49   0.22   0.88   2.13  2093 1.01
eta[10]        0.76    0.01 0.91  -1.01   0.16   0.75   1.36   2.56  5000 1.00
eta[11]        0.73    0.01 0.97  -1.20   0.07   0.73   1.40   2.59  5000 1.00
eta[12]        0.17    0.01 0.97  -1.72  -0.48   0.19   0.84   2.03  5000 1.00
eta[13]        0.00    0.01 1.00  -1.94  -0.67  -0.02   0.66   1.97  5000 1.00
eta[14]       -0.02    0.01 0.98  -1.91  -0.68  -0.01   0.64   1.84  5000 1.00
eta[15]        0.00    0.01 1.00  -1.96  -0.69  -0.02   0.65   2.02  5000 1.00
mu_0          34.99    0.07 3.38  28.29  32.75  35.04  37.21  41.56  2148 1.00
mu_b           1.40    0.01 0.43   0.55   1.11   1.40   1.68   2.28  2223 1.00
f[1]          40.07    0.02 1.41  37.13  39.20  40.11  40.94  42.93  5000 1.00
f[2]          41.99    0.02 1.16  39.44  41.30  42.08  42.78  44.06  3577 1.00
f[3]          42.11    0.02 1.11  39.72  41.44  42.19  42.85  44.13  5000 1.00
f[4]          40.35    0.01 1.05  38.13  39.70  40.37  41.04  42.37  5000 1.00
f[5]          37.49    0.02 1.08  35.38  36.80  37.49  38.18  39.68  5000 1.00
f[6]          35.11    0.03 1.39  32.83  34.11  34.94  35.94  38.26  2122 1.00
f[7]          35.99    0.05 1.57  32.99  34.96  35.94  36.95  39.25   920 1.01
f[8]          40.70    0.02 1.10  38.56  40.00  40.69  41.41  42.87  4150 1.00
f[9]          45.94    0.02 1.10  43.69  45.23  45.96  46.70  47.99  2672 1.00
f[10]         49.87    0.02 1.08  47.73  49.15  49.86  50.59  52.00  2352 1.00
f[11]         53.17    0.02 1.10  50.87  52.51  53.22  53.90  55.20  5000 1.00
f[12]         55.50    0.03 1.49  52.22  54.62  55.65  56.52  58.10  2674 1.00
f[13]         56.37    0.04 2.59  51.02  54.71  56.42  58.15  61.26  5000 1.00
f[14]         56.53    0.06 3.92  48.75  53.93  56.57  59.15  64.11  3893 1.00
f[15]         56.90    0.08 4.80  47.55  53.66  56.89  60.03  66.23  3491 1.00
y_pred[1]     40.09    0.03 2.21  35.55  38.78  40.12  41.47  44.46  5000 1.00
y_pred[2]     41.97    0.03 2.08  37.41  40.78  42.11  43.24  45.87  5000 1.00
y_pred[3]     42.07    0.03 2.06  37.62  40.90  42.17  43.39  46.01  4622 1.00
y_pred[4]     40.37    0.03 2.01  36.12  39.22  40.37  41.61  44.32  5000 1.00
y_pred[5]     37.46    0.03 1.99  33.42  36.20  37.48  38.66  41.55  4954 1.00
y_pred[6]     35.11    0.04 2.19  31.30  33.68  34.89  36.37  40.12  3319 1.00
y_pred[7]     36.01    0.06 2.32  32.19  34.38  35.82  37.38  41.16  1598 1.00
y_pred[8]     40.70    0.03 2.05  36.54  39.44  40.72  41.90  44.90  5000 1.00
y_pred[9]     45.96    0.03 2.04  41.60  44.75  46.06  47.23  49.88  5000 1.00
y_pred[10]    49.87    0.03 2.05  45.76  48.62  49.82  51.17  54.05  4480 1.00
y_pred[11]    53.22    0.03 2.06  48.98  51.98  53.30  54.43  57.22  5000 1.00
y_pred[12]    55.53    0.04 2.25  50.50  54.22  55.69  56.96  59.55  4060 1.00
y_pred[13]    56.35    0.04 3.08  49.90  54.44  56.43  58.41  62.17  5000 1.00
y_pred[14]    56.53    0.07 4.28  47.90  53.64  56.66  59.38  64.76  4141 1.00
y_pred[15]    56.90    0.08 5.09  47.04  53.49  56.81  60.32  66.59  3727 1.00
log_lik[1]    -1.70    0.02 0.64  -3.22  -1.96  -1.63  -1.35  -0.53   830 1.01
log_lik[2]    -1.68    0.02 0.59  -3.05  -1.96  -1.62  -1.35  -0.54   820 1.01
log_lik[3]    -1.62    0.02 0.53  -2.79  -1.88  -1.59  -1.33  -0.51   741 1.01
log_lik[4]    -1.60    0.02 0.49  -2.61  -1.84  -1.59  -1.35  -0.48   767 1.01
log_lik[5]    -1.70    0.02 0.58  -3.04  -1.94  -1.66  -1.41  -0.53   848 1.01
log_lik[6]    -1.80    0.02 0.73  -3.54  -2.16  -1.69  -1.34  -0.55   989 1.01
log_lik[7]    -3.65    0.05 1.61  -7.31  -4.57  -3.52  -2.56  -0.76  1089 1.00
log_lik[8]    -1.82    0.02 0.66  -3.45  -2.08  -1.74  -1.48  -0.54  1023 1.01
log_lik[9]    -1.82    0.02 0.68  -3.40  -2.14  -1.77  -1.46  -0.48   864 1.01
log_lik[10]   -2.01    0.02 0.79  -3.87  -2.32  -1.88  -1.58  -0.56  1116 1.01
log_lik[11]   -1.67    0.02 0.57  -2.99  -1.91  -1.62  -1.36  -0.50   863 1.01
log_lik[12]   -1.87    0.03 0.84  -3.97  -2.18  -1.71  -1.38  -0.51  1119 1.00
lp__         -20.82    0.20 4.32 -29.40 -23.34 -20.79 -18.36 -11.40   482 1.02

Samples were drawn using NUTS(diag_e) at Sat Dec  9 11:11:30 2017.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

All $\hat{R}$ values are below 1.1 so our MCMC simulations convergences well.

Tree Depth

source("stan_utility.R")
check_treedepth(fit_gp)

[1] "0 of 5000 iterations saturated the maximum tree depth of 10 (0%)"

Monte Carlo simulations did not hit the maximum trajectory length, so we don’t need to increase the tree depth.

E-BFMI

check_energy(fit_gp)

[1] "E-BFMI indicated no pathological behavior"

Divergence

check_div(fit_gp)

[1] "21 of 5000 iterations ended with a divergence (0.42%)"
[1] "  Try running with larger adapt_delta to remove the divergences"

Divergences does not vanish completely even though adapt_delta value is 0.999. Let’s check the posterior of length_param (length-scale) and sigma, which is the unexpleinable variance in the data, to see if divergences are false positives or caused by an actual problem.

c_dark <- c("#8F272780")
green <- c("#00FF0080")
partition <- partition_div(fit_gp)
div_params <- partition[[1]]
nondiv_params <- partition[[2]]
par(mar = c(4, 4, 0.5, 0.5))
plot(nondiv_params$length_param, nondiv_params$sigma,
     col=c_dark, pch=16, cex=0.8, xlab="length_param", ylab="sigma",
     xlim=c(0, 4), ylim=c(0,7))
points(div_params$length_param, div_params$sigma,
       col=green, pch=16, cex=0.8)

As can be seen from the plot, all the divergences are clustered to the bottom left where the posterior abrupts quickly. So there is indeed problem with our model. Monte Carlo is not able to explore that area properly so our posterior is slightly biased.

Fixing the model

One issue we see in the posterior plot above is that some simulations of the length parameter in our exponential gaussian kernel fall below 1. This is a possible cause of overfit that might result in divergences, since our data contains only yearly observations.

We try to fix our model by adjusting our prior for length_param (length-scale). At the moment we have $Weibull(10,2)$ , which might give us values between 0 and 1. We can fix this by moving our prior further away from 1 and clipping the prior so that it can only get values higher than 1. Our prior for length-scale is $Weibull(15,3)$ with minimum of 1. After adjusting the prior, the convergence analytics are performed again:

check_treedepth(fit_gp)

[1] "0 of 5000 iterations saturated the maximum tree depth of 10 (0%)"

check_energy(fit_gp)

[1] "E-BFMI indicated no pathological behavior"

check_div(fit_gp)

[1] "0 of 5000 iterations ended with a divergence (0%)"

Changing the prior of the length parameter seems to fix our divergence issues.

Visualization

sample = extract(fit_gp)
sample_preds = data.frame(sample$y_pred)
  
sample_quantiles <- apply(sample_preds,2,function(x){quantile(x,c(0.05,0.5,0.95))})
sample_data <- data.frame(low = sample_quantiles[1,],
                          med = sample_quantiles[2,],
                          high = sample_quantiles[3,],
                          x = x)
training_data <- data.frame(x=x_is, obs=y_is)
validation_data <- data.frame(x=x_os, obs=y_os)
cols <- c("GP median"="black","GP 90% percentile"="lightblue","True debt"="red")
ggplot(sample_data, aes(x=x, y=med)) +
  geom_ribbon(aes(x=x,ymin=low,ymax=high,colour="GP 90% percentile"),alpha=0.6,fill="lightblue") +
  geom_line(aes(x=x,y=med,colour="GP median")) +
  geom_point() +
  geom_line(data=training_data, aes(x=x, y=obs,colour="True debt")) +
  geom_line(data=validation_data, aes(x=x, y=obs,colour="True debt"), linetype="dashed") +
  xlab('year') + 
  ylab('Debt (% of GDP)') + 
  ggtitle(sprintf('Government debt, country = Finland')) +
  scale_colour_manual(name="Colour",values=cols) +
  theme(plot.title = element_text(hjust = 0.5))

The empirical data stays withing the 90% confidence interval in all cases apart from the Year 7 (corresponding to year 2008). This is the year of bank crisis and is therefore an anomaly to the process as does not follow the natural progression explained by the model.

Changing the prior for the length param smoothens the median compared to the original prior choice. This is natural, as increasing the length parameter increases time scale at which the correlation decays. In other words, the effect of the time steps further away contribute more to the value, thus reducing the frequency of the signal.

The out of sample predictions (dashed red line) fall also within the cofidence interval. However, the dataset is small and therefore making accurate predictions is not a trivial task without additional knowledge.

Posterior predictive checking

library(bayesplot)

This is bayesplot version 1.4.0
- Plotting theme set to bayesplot::theme_default()
- Online documentation at mc-stan.org/bayesplot

sample = extract(fit_gp)
ppc_dens_overlay(y=y_is, yrep=sample$y_pred[,1:(length(y) - ss)])+ggtitle("Posterior predictive for Finland") + xlab("Debt (% of GDP)")

Simulated draws from the posterior distribution follow the empirical data distribution roughly. Huge variance in draws at 40-45 is due the fact that our informative prior $Weibull(15,3)$ smoothens the dip (2008 financial crisis) quite heavily so our posterior gives more probablity mass around 40-45 compared to the empirical data.

Sensitivity analysis

We are focusing our sensitivity analysis on different priors for length_param which is the length-scale parameter for gaussian process. Ideally the sensitivity analysis should be performed for the other priors as well, however, due to time constraints we have chose the the length parameter as changing it is expected to have the most notable effect.

Length_param

Vague prior: $length\_param \sim \Gamma(4,4)$
Non-informative prior: $length\_param \sim Uniform(0,10)$
Informative prior: $length\_param \sim Weibull(15, 3)$

cc = 'FI'
y = data[,cc]
x = 1:length(y)
ss = 3
y_is = y[1:(length(y) - ss)]
x_is = x[1:(length(x) - ss)]
N_is = length(x_is)
y_os = tail(y, (ss + 1))
x_os = tail(x, (ss + 1))
x_pred = tail(x, ss)
N_pred = length(x_pred)
prior_check <- function(model_name) {
  fit_gp_tmp <- stan(file=model_name,
               data=list(x_is = x_is,
                         y_is = y_is,
                         N_is = N_is,
                         x_pred = x_pred,
                         N_pred),
               iter=1000,
               chains=10,control=list(adapt_delta=0.99,max_treedepth=15))
  sample <- extract(fit_gp_tmp)
  sample_preds <- data.frame(sample$y_pred)
  
  sample_quantiles <- apply(sample_preds,2,function(x){quantile(x,c(0.05,0.5,0.95))})
  
  sample_data <- data.frame(low = sample_quantiles[1,],
                          med = sample_quantiles[2,],
                          high = sample_quantiles[3,],
                          x = x)
  training_data <- data.frame(x=x_is, obs=y_is)
  validation_data <- data.frame(x=x_os, obs=y_os)
  return(list(s=sample_data, t=training_data, v=validation_data,m=fit_gp_tmp))
}
sample1 <- prior_check('debt_gp_linear.stan')
sample2 <- prior_check('debt_gp_linear2.stan') 
sample3 <- prior_check('debt_gp_linear3.stan')
#cols <- c("Weibull(10,2)"="black","Gamma(4,4)"="green","Uniform prior"="blue","True debt"="red")
ggplot(sample1$s, aes(x=x, y=med)) +
  geom_line(aes(x=x,y=med,colour="Weibull(15,3)")) +
  geom_line(data=sample2$s, aes(x=x,y=med,colour="Gamma(4,4)")) + 
  geom_line(data=sample3$s, aes(x=x,y=med,colour="Uniform prior")) + 
  geom_line(data=sample1$t, aes(x=x, y=obs,colour="True debt")) +
  geom_line(data=sample1$v, aes(x=x, y=obs,colour="True debt"), linetype="dashed") +
  xlab('year') + 
  ylab('Debt (% of GDP)') + 
  ggtitle(sprintf('Government debt, country = Finland')) +
  #scale_colour_manual(name="Colour",values=cols) +
  theme(plot.title = element_text(hjust = 0.5))

As can be seen from the plot, prior information about length-scale has an noticable effect on the outcome. Our choice for informative prior smoothens the curve the most. Using non-informative and vague prior overfits model quite much. Based only on visual analysis our informative prior seems the most reasonable, because the data contain datapoints that are basically impossible to explain based only on this data (financial crisis). However, our model is sensitive to the choice of prior for length scale.

Model comparison for different priors (LOO)

Informative prior

library(loo)

This is loo version 1.1.0

ll <- extract_log_lik(sample1$m)
l <- loo(ll)

Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.

Computed from 5000 by 12 log-likelihood matrix

         Estimate  SE
elpd_loo    -30.8 4.2
p_loo         5.8 2.1
looic        61.6 8.4

Pareto k diagnostic values:
                         Count  Pct 
(-Inf, 0.5]   (good)     6     50.0%
 (0.5, 0.7]   (ok)       5     41.7%
   (0.7, 1]   (bad)      1      8.3%
   (1, Inf)   (very bad) 0      0.0%
See help('pareto-k-diagnostic') for details.

vague prior

ll <- extract_log_lik(sample2$m)
l <- loo(ll)

Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.

Computed from 5000 by 12 log-likelihood matrix

         Estimate  SE
elpd_loo    -21.4 2.3
p_loo        16.1 2.1
looic        42.9 4.6

Pareto k diagnostic values:
                         Count  Pct 
(-Inf, 0.5]   (good)     0      0.0%
 (0.5, 0.7]   (ok)       6     50.0%
   (0.7, 1]   (bad)      6     50.0%
   (1, Inf)   (very bad) 0      0.0%
See help('pareto-k-diagnostic') for details.

non-informative prior

ll <- extract_log_lik(sample3$m)
l <- loo(ll)

Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.

Computed from 5000 by 12 log-likelihood matrix

         Estimate  SE
elpd_loo    -25.4 3.4
p_loo        14.8 3.1
looic        50.9 6.8

Pareto k diagnostic values:
                         Count  Pct 
(-Inf, 0.5]   (good)     0      0.0%
 (0.5, 0.7]   (ok)       6     50.0%
   (0.7, 1]   (bad)      6     50.0%
   (1, Inf)   (very bad) 0      0.0%
See help('pareto-k-diagnostic') for details.

Comparing the models with different prior using psis-LOO is not possible, because models for non informative prior and vague prior get Pareto k values >1.0. Therefore we should use more computationally expenisive methods, such as standard cross-validation, to evalute the performance of our models. However pareto-k values for our informative prior model are quite reasonble: only 2 values are above 0.7 so it is possible to compare that model to our hierarchical model.

Second model (hierarchical)

Our second model tries to incorporate more information to each model by adding hierarchical structure. The idea is to add global trend and gaussian process to the separate model in the previous section. The global components will be affected by all the EU countries in our dataset.

Similarly to the separate model, the components are formed as:

μ_{i} = μ_{0} + μ_{b} \cdot x_{i}

$\begin{equation} \mu_i = \mu_0 + \mu_b \cdot x_i \end{equation}$

And:

K_{i, j} = α^{2} \exp (- \frac{1}{2 l^{2}} (x_{i} - x_{j})^{2})

$\begin{equation} K_{i,j} = \alpha^2 \exp \Big( - \frac{1}{2 l^2} (x_i - x_j)^2 \Big) \end{equation}$ Together these compose the final

\begin{aligned} f & \sim M u l t i N o r m a l (μ_{s} + μ_{g}, K_{c} + K_{g}) \\ y_{i} & \sim N o r m a l (f_{i}, σ) \end{aligned}

$\begin{align} f & \sim MultiNormal(\mu_s + \mu_g, K_c + K_g) \\ y_i & \sim Normal(f_i, \sigma) \end{align}$

Stan (hierarchical)

data {
  // Samples used for fitting
    int<lower=1> N_is;
    int country_ind_is[N_is];
    int year_ind_is[N_is];
    
    vector[N_is] y_is;
    
    int<lower=1> N_years_is;
    int<lower=1> N_years_pred;

    int<lower=1> N_countries;
    // The amount of predictions
}

transformed data {
  // All data points
  int<lower=1> N = N_is + N_years_pred * N_countries;
  int<lower=1> year_ind[N];
  int<lower=1> country_ind[N];
  
  int<lower=1> N_years = N_years_is + N_years_pred;
  real<lower=1> years[N_years];

  real jitter = 1e-9;

  for (year in 1:N_years) years[year] = year;

  // Add in sample year indices
  for (n in 1:N_is){
    year_ind[n] = year_ind_is[n];
    country_ind[n] = country_ind_is[n];
  }
    
  // Create out sample year indices
  for (year in 1:N_years_pred)
    for (country in 1:N_countries) {
      country_ind[N_is + (year - 1) * N_countries + country] = country;
      year_ind[N_is + (year - 1) * N_countries + country] = N_years_is + year;
    }

}

parameters {
  real<lower=0> alpha;   // Strength parameter for GP
  real<lower=0> length_param;  // Lag parameter
  
  real<lower=0> global_alpha;
  real<lower=0> global_length;
  
  real<lower=0> sigma;  // Random, unexplained noise

  matrix[N_years, N_countries] eta;   // Gaussian 0,1 for generating the covariance from Cholesky decompose
  vector[N_years] global_eta;
  
  real mu_0[N_countries];
  real mu_b[N_countries];
  
  real<lower=0> global_mu_0;
  real global_mu_b;
}

transformed parameters {
  
  matrix[N_years, N_countries] GP_country;
  vector[N_years] GP_global;
  
  vector[N] f;
  matrix[N_years, N_countries] linear_country;
  
  real linear_global[N_years];
  
  {
    matrix[N_years, N_years] K = cov_exp_quad(years, alpha, length_param) + diag_matrix(rep_vector(jitter, N_years));
    matrix[N_years, N_years] K_global = cov_exp_quad(years, global_alpha, global_length) + diag_matrix(rep_vector(jitter, N_years));
  
    matrix[N_years, N_years] L_K = cholesky_decompose(K);
    matrix[N_years, N_years] L_K_global = cholesky_decompose(K_global);
    
    GP_country = L_K * eta;
    GP_global = L_K_global * global_eta;
    
    for (year in 1:N_years) {
      linear_global[year] = global_mu_0 + global_mu_b * year;
      
      for (country in 1:N_countries){
        linear_country[year, country] = mu_0[country] + mu_b[country] * year;
      }
    }
      
    for (n in 1:N){
      f[n] = linear_country[year_ind[n], country_ind[n]]
        + linear_global[year_ind[n]]
        + GP_country[year_ind[n], country_ind[n]]
        + GP_global[year_ind[n]];
    }
  }
}

model {
  alpha ~ normal(0, 1);
  global_alpha ~ normal(0, 1);
  
  //length_param ~ gamma(4,4); // Vague prior
  length_param ~ weibull(15, 3 ); // More informative prior
  global_length ~ weibull(15, 3);
  
  sigma ~ normal(0, 1);
  
  to_vector(eta) ~ normal(0, 1);
  to_vector(global_eta) ~ normal(0, 1);
  
  mu_0 ~ normal(0, 100);
  mu_b ~ normal(0, 5);
  
  global_mu_0 ~ normal(0, 100);
  global_mu_b ~ normal(0, 5);
  
  y_is ~ normal(f[1:N_is], sigma);
}

generated quantities {
  vector[N] y_pred;
  vector[N_is] log_lik;
  vector[N] year_indicators;
  vector[N] country_indicators;

  year_indicators = to_vector(year_ind);
  country_indicators = to_vector(country_ind);

  for (n in 1:N) y_pred[n] = normal_rng(f[n], sigma);
  for (n in 1:N_is) log_lik[n] = normal_lpdf(y_is[n] | f[n], sigma);
}

Model

data = read.csv('data.csv', row.names=1)
data = data.frame(t(data))
y_is = unlist(data,use.names = F)
country_ind_is = rep(1:ncol(data),each=15)
year_ind_is = rep(1:15,ncol(data))
N_is = length(y_is)
N_countries = ncol(data)
N_years_pred = 3
N_years_is = 15
# Simple GP + linear model
gp_model_hier = 'debt_gp_hierarchical.stan'
fit_gp_h <- stan(file=gp_model_hier,
               data=list(y_is = y_is,
                         country_ind_is = country_ind_is,
                         year_ind_is = year_ind_is,
                         N_years_is = 15,
                         N_years_pred = N_years_pred,
                         N_is = N_is,
                         N_countries = N_countries),
               iter=1000,
               chains=10)
sample <- extract(fit_gp_h)

sample_preds <- data.frame(sample$y_pred[,sample$country_indicators[1,] == match("FI",colnames(data))])
  
sample_quantiles <- apply(sample_preds,2,function(x){quantile(x,c(0.05,0.5,0.95))})
x = 1:18
sample_data <- data.frame(low = sample_quantiles[1,],
                          med = sample_quantiles[2,],
                          high = sample_quantiles[3,],
                          x = x)
training_data <- data.frame(x=1:15, obs=data[,"FI"])
#validation_data <- data.frame(x=15:18, obs=y_os)
cols <- c("GP median"="black","GP 90% percentile"="lightblue","True debt"="red")
ggplot(sample_data, aes(x=x, y=med)) +
  geom_ribbon(aes(x=x,ymin=low,ymax=high,colour="GP 90% percentile"),alpha=0.6,fill="lightblue") +
  geom_line(aes(x=x,y=med,colour="GP median")) +
  geom_point() +
  geom_line(data=training_data, aes(x=x, y=obs,colour="True debt")) +
  #geom_line(data=validation_data, aes(x=x, y=obs,colour="True debt"), linetype="dashed") +
  xlab('year') + 
  ylab('Debt (% of GDP)') + 
  ggtitle(sprintf('Government debt, country = Finland (hierarchical model)')) +
  scale_colour_manual(name="Colour",values=cols) +
  theme(plot.title = element_text(hjust = 0.5))

By adding more data using hierarcical gaussian processes, our model gives tighter error bounds than the separate model for Finland. As expected the uncertainty increases after there are no observed data points.

Convergence diagnostics

$\hat{R}$ values

fit_gp_h

Inference for Stan model: debt_gp_hierarchical.
10 chains, each with iter=1000; warmup=500; thin=1; 
post-warmup draws per chain=500, total post-warmup draws=5000.

                           mean se_mean    sd    2.5%     25%     50%     75%   97.5% n_eff Rhat
alpha                      6.85    0.01  0.39    6.14    6.58    6.84    7.12    7.66  1422 1.00
length_param               2.37    0.00  0.12    2.14    2.29    2.36    2.44    2.60   549 1.01
global_alpha               3.44    0.01  0.52    2.52    3.07    3.41    3.77    4.54  2938 1.00
global_length              1.69    0.00  0.15    1.38    1.59    1.69    1.79    1.98  1508 1.01
sigma                      1.72    0.00  0.09    1.57    1.66    1.72    1.78    1.90  1482 1.00
eta[1,1]                  -0.37    0.01  0.85   -2.05   -0.92   -0.37    0.18    1.32  3960 1.00
eta[1,2]                   0.38    0.01  0.86   -1.36   -0.20    0.41    0.96    2.03  3763 1.00
eta[1,3]                   1.29    0.01  0.87   -0.40    0.70    1.27    1.88    3.04  3786 1.00
eta[1,4]                   0.95    0.01  0.87   -0.78    0.37    0.96    1.54    2.64  4085 1.00
eta[1,5]                  -0.29    0.01  0.87   -2.02   -0.86   -0.29    0.30    1.39  4065 1.00
eta[1,6]                  -0.68    0.01  0.86   -2.31   -1.26   -0.70   -0.13    1.10  3756 1.00
eta[1,7]                   0.37    0.01  0.86   -1.31   -0.25    0.37    0.97    2.06  3705 1.00
eta[1,8]                  -0.13    0.01  0.89   -1.86   -0.73   -0.15    0.49    1.66  3796 1.00
eta[1,9]                  -0.15    0.01  0.88   -1.85   -0.75   -0.15    0.44    1.57  4239 1.00
eta[1,10]                 -0.16    0.01  0.89   -1.91   -0.77   -0.16    0.46    1.55  3606 1.00
eta[1,11]                  0.27    0.01  0.86   -1.41   -0.31    0.26    0.82    1.95  3848 1.00
eta[1,12]                  0.97    0.02  0.87   -0.79    0.40    0.98    1.55    2.72  3184 1.00
eta[1,13]                 -0.15    0.01  0.88   -1.86   -0.74   -0.15    0.46    1.56  4222 1.00
eta[1,14]                 -0.17    0.01  0.85   -1.85   -0.74   -0.18    0.40    1.51  3684 1.00
eta[1,15]                  0.35    0.01  0.87   -1.29   -0.24    0.35    0.93    2.06  3951 1.00
eta[1,16]                 -0.01    0.01  0.86   -1.69   -0.58   -0.01    0.57    1.66  3747 1.00
eta[1,17]                  0.11    0.02  0.86   -1.59   -0.47    0.10    0.68    1.81  3069 1.00
eta[1,18]                 -0.72    0.01  0.88   -2.41   -1.32   -0.72   -0.12    1.00  3579 1.00
eta[1,19]                 -0.34    0.01  0.86   -2.02   -0.92   -0.34    0.22    1.34  3418 1.00
eta[1,20]                  0.00    0.01  0.86   -1.74   -0.58    0.00    0.57    1.68  3897 1.00
eta[1,21]                  0.01    0.01  0.88   -1.73   -0.59    0.01    0.59    1.77  3495 1.00
eta[1,22]                 -0.31    0.01  0.86   -2.01   -0.88   -0.31    0.26    1.36  4004 1.00
eta[1,23]                 -0.07    0.01  0.86   -1.74   -0.66   -0.07    0.50    1.63  3575 1.00
eta[1,24]                 -0.31    0.01  0.88   -2.06   -0.90   -0.33    0.28    1.47  3550 1.00
eta[1,25]                 -0.08    0.02  0.87   -1.79   -0.66   -0.07    0.51    1.63  3085 1.00
eta[1,26]                 -0.14    0.01  0.87   -1.85   -0.73   -0.13    0.45    1.56  3883 1.00
eta[1,27]                 -0.14    0.01  0.87   -1.84   -0.72   -0.14    0.43    1.57  3952 1.00
eta[1,28]                  0.32    0.01  0.87   -1.43   -0.26    0.32    0.90    2.05  3353 1.00
eta[1,29]                 -0.01    0.01  0.86   -1.68   -0.61   -0.03    0.57    1.67  3491 1.00
eta[1,30]                  0.49    0.02  0.89   -1.27   -0.10    0.48    1.07    2.26  3148 1.00
eta[1,31]                  0.55    0.01  0.87   -1.12   -0.04    0.55    1.14    2.29  3616 1.00
eta[1,32]                 -0.25    0.02  0.88   -1.97   -0.83   -0.25    0.35    1.44  3403 1.00
eta[2,1]                  -0.18    0.01  0.50   -1.17   -0.51   -0.18    0.16    0.80  5000 1.00
eta[2,2]                  -0.67    0.01  0.51   -1.68   -1.02   -0.67   -0.32    0.33  3493 1.00
eta[2,3]                  -1.15    0.01  0.52   -2.21   -1.49   -1.14   -0.80   -0.12  5000 1.00
eta[2,4]                   0.56    0.01  0.51   -0.44    0.22    0.56    0.90    1.57  5000 1.00
eta[2,5]                   0.37    0.01  0.51   -0.63    0.02    0.38    0.71    1.34  5000 1.00
eta[2,6]                   0.88    0.01  0.52   -0.13    0.53    0.89    1.23    1.89  5000 1.00
eta[2,7]                  -0.04    0.01  0.52   -1.08   -0.39   -0.03    0.31    0.97  3539 1.00
eta[2,8]                   0.19    0.01  0.50   -0.83   -0.15    0.19    0.53    1.14  3817 1.00
eta[2,9]                   0.18    0.01  0.49   -0.79   -0.15    0.18    0.53    1.15  5000 1.00
eta[2,10]                  0.05    0.01  0.50   -0.93   -0.28    0.05    0.39    1.02  3889 1.00
eta[2,11]                 -1.33    0.01  0.50   -2.30   -1.68   -1.33   -0.98   -0.36  5000 1.00
eta[2,12]                 -1.21    0.01  0.51   -2.21   -1.56   -1.21   -0.86   -0.21  5000 1.00
eta[2,13]                  0.21    0.01  0.51   -0.79   -0.13    0.21    0.54    1.19  5000 1.00
eta[2,14]                  0.21    0.01  0.49   -0.73   -0.12    0.21    0.54    1.16  5000 1.00
eta[2,15]                  0.50    0.01  0.50   -0.48    0.16    0.50    0.83    1.48  5000 1.00
eta[2,16]                  0.63    0.01  0.50   -0.36    0.29    0.63    0.96    1.59  5000 1.00
eta[2,17]                  0.10    0.01  0.51   -0.92   -0.24    0.09    0.44    1.11  5000 1.00
eta[2,18]                  0.15    0.01  0.50   -0.84   -0.19    0.15    0.48    1.14  5000 1.00
eta[2,19]                 -1.25    0.01  0.52   -2.26   -1.60   -1.24   -0.89   -0.24  5000 1.00
eta[2,20]                 -0.49    0.01  0.50   -1.47   -0.82   -0.48   -0.15    0.48  4072 1.00
eta[2,21]                 -0.54    0.01  0.50   -1.47   -0.88   -0.54   -0.20    0.44  4381 1.00
eta[2,22]                 -0.04    0.01  0.51   -1.04   -0.39   -0.04    0.31    0.96  5000 1.00
eta[2,23]                  0.00    0.01  0.50   -0.95   -0.34    0.00    0.32    0.97  5000 1.00
eta[2,24]                  1.81    0.01  0.51    0.81    1.46    1.81    2.16    2.79  3866 1.00
eta[2,25]                  0.38    0.01  0.52   -0.63    0.03    0.38    0.73    1.38  5000 1.00
eta[2,26]                  0.64    0.01  0.51   -0.36    0.29    0.63    0.99    1.63  5000 1.00
eta[2,27]                 -0.39    0.01  0.50   -1.39   -0.73   -0.40   -0.06    0.62  5000 1.00
eta[2,28]                 -0.79    0.01  0.49   -1.79   -1.12   -0.79   -0.46    0.15  5000 1.00
eta[2,29]                  0.55    0.01  0.50   -0.41    0.22    0.55    0.89    1.53  3211 1.00
eta[2,30]                 -0.51    0.01  0.51   -1.49   -0.84   -0.52   -0.17    0.49  5000 1.00
eta[2,31]                 -0.16    0.01  0.49   -1.14   -0.48   -0.15    0.18    0.77  3751 1.00
eta[2,32]                 -0.25    0.01  0.51   -1.26   -0.59   -0.24    0.10    0.75  5000 1.00
eta[3,1]                   0.48    0.01  0.69   -0.84    0.02    0.48    0.95    1.84  4437 1.00
eta[3,2]                  -0.32    0.01  0.71   -1.66   -0.82   -0.32    0.17    1.05  4165 1.00
eta[3,3]                  -0.59    0.01  0.71   -1.98   -1.06   -0.58   -0.12    0.81  3763 1.00
eta[3,4]                   0.04    0.01  0.70   -1.33   -0.43    0.05    0.52    1.41  4121 1.00
eta[3,5]                   0.20    0.01  0.72   -1.16   -0.28    0.18    0.68    1.61  3375 1.00
eta[3,6]                   0.64    0.01  0.71   -0.71    0.15    0.62    1.13    2.07  4019 1.00
eta[3,7]                  -0.82    0.01  0.70   -2.23   -1.29   -0.81   -0.34    0.51  3636 1.00
eta[3,8]                   0.09    0.01  0.71   -1.30   -0.39    0.09    0.57    1.49  3763 1.00
eta[3,9]                   0.08    0.01  0.71   -1.27   -0.41    0.07    0.56    1.47  3814 1.00
eta[3,10]                  0.41    0.01  0.71   -0.96   -0.07    0.41    0.88    1.81  3677 1.00
eta[3,11]                 -0.08    0.01  0.70   -1.46   -0.55   -0.08    0.40    1.30  3947 1.00
eta[3,12]                 -0.86    0.01  0.72   -2.25   -1.36   -0.87   -0.37    0.54  3998 1.00
eta[3,13]                  0.06    0.01  0.70   -1.35   -0.41    0.06    0.54    1.41  5000 1.00
eta[3,14]                  0.05    0.01  0.70   -1.31   -0.41    0.06    0.53    1.41  3720 1.00
eta[3,15]                 -0.11    0.01  0.72   -1.50   -0.60   -0.11    0.38    1.31  4150 1.00
eta[3,16]                 -0.12    0.01  0.70   -1.49   -0.59   -0.13    0.35    1.26  3854 1.00
eta[3,17]                  0.00    0.01  0.69   -1.34   -0.46    0.01    0.46    1.34  4035 1.00
eta[3,18]                  0.48    0.01  0.71   -0.85   -0.01    0.46    0.94    1.87  3506 1.00
eta[3,19]                 -1.80    0.01  0.70   -3.18   -2.30   -1.80   -1.32   -0.41  3979 1.00
eta[3,20]                  0.35    0.01  0.70   -1.01   -0.13    0.35    0.84    1.70  4130 1.00
eta[3,21]                  0.07    0.01  0.70   -1.30   -0.40    0.08    0.54    1.41  4236 1.00
eta[3,22]                  0.16    0.01  0.70   -1.22   -0.33    0.17    0.61    1.52  4166 1.00
eta[3,23]                 -0.83    0.01  0.72   -2.21   -1.34   -0.83   -0.35    0.56  2798 1.00
eta[3,24]                  0.08    0.01  0.72   -1.31   -0.42    0.06    0.57    1.51  3642 1.00
eta[3,25]                 -0.48    0.01  0.69   -1.85   -0.95   -0.46   -0.01    0.84  4008 1.00
eta[3,26]                  0.38    0.01  0.69   -0.99   -0.08    0.37    0.85    1.74  4548 1.00
eta[3,27]                  0.45    0.01  0.74   -0.98   -0.04    0.45    0.95    1.85  4148 1.00
eta[3,28]                 -0.39    0.01  0.71   -1.81   -0.86   -0.39    0.10    1.00  4051 1.00
eta[3,29]                  0.53    0.01  0.71   -0.86    0.05    0.51    1.00    1.94  3902 1.00
eta[3,30]                 -0.08    0.01  0.71   -1.44   -0.56   -0.07    0.40    1.30  4126 1.00
eta[3,31]                 -0.57    0.01  0.69   -1.91   -1.03   -0.57   -0.11    0.80  3831 1.00
 [ reached getOption("max.print") -- omitted 4538 rows ]

Samples were drawn using NUTS(diag_e) at Sat Dec  9 11:28:07 2017.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

$\hat{R}$ values are all 1 so MCMC simulations convergences.

Tree depth

check_treedepth(fit_gp_h)

[1] "4754 of 5000 iterations saturated the maximum tree depth of 10 (95.08%)"
[1] "  Run again with max_depth set to a larger value to avoid saturation"

Our model hit the maximum trajectory length in MCMC simulations multiple times so we need to increase our treedepth to fix this issue.

E-BFMI

check_energy(fit_gp_h)

[1] "E-BFMI indicated no pathological behavior"

Divergence

check_div(fit_gp_h)

[1] "64 of 5000 iterations ended with a divergence (1.28%)"
[1] "  Try running with larger adapt_delta to remove the divergences"

c_dark <- c("#8F272780")
green <- c("#00FF0080")
partition <- partition_div(fit_gp_h)
div_params <- partition[[1]]
nondiv_params <- partition[[2]]
par(mar = c(4, 4, 0.5, 0.5))
plot(nondiv_params$global_alpha, nondiv_params$sigma,
     col=c_dark, pch=16, cex=0.8, xlab="length_param", ylab="sigma")
points(div_params$global_alpha, div_params$sigma,
       col=green, pch=16, cex=0.8)

We tried different parameter comparisons but we were not able to find any clear problem based on the plots, why our divergences does not vanish. We will try to fix the issue by tuning the MCMC parameters adapt delta and tree depth.

Fixing the hierarchical model

fit_gp_h <- stan(file=gp_model_hier,
               data=list(y_is = y_is,
                         country_ind_is = country_ind_is,
                         year_ind_is = year_ind_is,
                         N_years_is = 15,
                         N_years_pred = N_years_pred,
                         N_is = N_is,
                         N_countries = N_countries),
               iter=1000,
               chains=10,control=list(adapt_delta=0.99,max_treedepth=15))

Convergence

print(fit_gp_h)

Inference for Stan model: debt_gp_hierarchical.
10 chains, each with iter=1000; warmup=500; thin=1; 
post-warmup draws per chain=500, total post-warmup draws=5000.

                           mean se_mean    sd    2.5%     25%     50%     75%   97.5% n_eff Rhat
alpha                      6.85    0.01  0.39    6.12    6.57    6.83    7.11    7.67  1865 1.00
length_param               2.36    0.00  0.12    2.14    2.28    2.36    2.44    2.61   669 1.01
global_alpha               3.43    0.01  0.51    2.50    3.06    3.40    3.76    4.50  5000 1.00
global_length              1.69    0.00  0.15    1.39    1.59    1.69    1.79    1.96  2133 1.00
sigma                      1.72    0.00  0.08    1.57    1.66    1.72    1.78    1.89  1425 1.00
eta[1,1]                  -0.35    0.01  0.87   -2.04   -0.92   -0.36    0.22    1.33  5000 1.00
eta[1,2]                   0.37    0.01  0.88   -1.35   -0.22    0.35    0.94    2.10  5000 1.00
eta[1,3]                   1.30    0.01  0.91   -0.49    0.69    1.30    1.90    3.13  5000 1.00
eta[1,4]                   0.96    0.01  0.88   -0.76    0.38    0.95    1.54    2.71  5000 1.00
eta[1,5]                  -0.29    0.01  0.87   -1.98   -0.87   -0.29    0.29    1.41  5000 1.00
eta[1,6]                  -0.70    0.01  0.87   -2.41   -1.28   -0.69   -0.11    1.01  5000 1.00
eta[1,7]                   0.36    0.01  0.88   -1.39   -0.23    0.36    0.95    2.11  5000 1.00
eta[1,8]                  -0.13    0.01  0.86   -1.85   -0.72   -0.13    0.46    1.54  5000 1.00
eta[1,9]                  -0.14    0.01  0.87   -1.85   -0.73   -0.15    0.45    1.57  5000 1.00
eta[1,10]                 -0.19    0.01  0.86   -1.86   -0.77   -0.19    0.39    1.48  5000 1.00
eta[1,11]                  0.29    0.01  0.87   -1.39   -0.29    0.29    0.88    2.03  5000 1.00
eta[1,12]                  0.96    0.01  0.88   -0.74    0.37    0.96    1.54    2.70  5000 1.00
eta[1,13]                 -0.17    0.01  0.88   -1.91   -0.76   -0.16    0.42    1.52  5000 1.00
eta[1,14]                 -0.17    0.01  0.86   -1.87   -0.74   -0.18    0.40    1.59  5000 1.00
eta[1,15]                  0.39    0.01  0.85   -1.27   -0.19    0.38    0.97    2.03  5000 1.00
eta[1,16]                 -0.01    0.01  0.87   -1.70   -0.59   -0.03    0.57    1.66  5000 1.00
eta[1,17]                  0.12    0.01  0.88   -1.55   -0.48    0.11    0.70    1.82  5000 1.00
eta[1,18]                 -0.72    0.01  0.86   -2.42   -1.31   -0.74   -0.12    1.00  5000 1.00
eta[1,19]                 -0.32    0.01  0.88   -1.99   -0.91   -0.33    0.26    1.46  5000 1.00
eta[1,20]                  0.03    0.01  0.89   -1.74   -0.58    0.04    0.64    1.76  5000 1.00
eta[1,21]                  0.00    0.01  0.88   -1.79   -0.58    0.00    0.58    1.69  5000 1.00
eta[1,22]                 -0.31    0.01  0.88   -2.03   -0.90   -0.30    0.26    1.40  5000 1.00
eta[1,23]                 -0.07    0.01  0.86   -1.74   -0.65   -0.08    0.53    1.62  5000 1.00
eta[1,24]                 -0.31    0.01  0.88   -2.06   -0.90   -0.32    0.28    1.39  5000 1.00
eta[1,25]                 -0.07    0.01  0.87   -1.78   -0.67   -0.08    0.52    1.60  5000 1.00
eta[1,26]                 -0.16    0.01  0.88   -1.85   -0.75   -0.16    0.44    1.59  5000 1.00
eta[1,27]                 -0.11    0.01  0.87   -1.80   -0.70   -0.11    0.47    1.58  5000 1.00
eta[1,28]                  0.33    0.01  0.86   -1.35   -0.25    0.33    0.91    2.04  5000 1.00
eta[1,29]                  0.00    0.01  0.84   -1.64   -0.57    0.00    0.57    1.64  5000 1.00
eta[1,30]                  0.49    0.01  0.87   -1.20   -0.11    0.49    1.10    2.21  5000 1.00
eta[1,31]                  0.57    0.01  0.86   -1.12   -0.02    0.56    1.14    2.26  5000 1.00
eta[1,32]                 -0.25    0.01  0.87   -1.94   -0.83   -0.26    0.33    1.46  5000 1.00
eta[2,1]                  -0.17    0.01  0.49   -1.15   -0.51   -0.18    0.16    0.78  5000 1.00
eta[2,2]                  -0.69    0.01  0.50   -1.65   -1.03   -0.69   -0.35    0.32  5000 1.00
eta[2,3]                  -1.13    0.01  0.52   -2.14   -1.47   -1.13   -0.78   -0.10  5000 1.00
eta[2,4]                   0.57    0.01  0.51   -0.42    0.22    0.57    0.91    1.59  5000 1.00
eta[2,5]                   0.36    0.01  0.50   -0.63    0.02    0.37    0.70    1.34  5000 1.00
eta[2,6]                   0.86    0.01  0.51   -0.14    0.53    0.87    1.20    1.87  5000 1.00
eta[2,7]                  -0.05    0.01  0.50   -1.04   -0.38   -0.05    0.29    0.95  5000 1.00
eta[2,8]                   0.18    0.01  0.52   -0.82   -0.17    0.18    0.54    1.20  5000 1.00
eta[2,9]                   0.19    0.01  0.51   -0.81   -0.15    0.19    0.53    1.20  5000 1.00
eta[2,10]                  0.07    0.01  0.51   -0.94   -0.26    0.07    0.40    1.06  5000 1.00
eta[2,11]                 -1.33    0.01  0.51   -2.32   -1.68   -1.33   -1.00   -0.32  5000 1.00
eta[2,12]                 -1.20    0.01  0.50   -2.19   -1.55   -1.20   -0.87   -0.23  5000 1.00
eta[2,13]                  0.21    0.01  0.50   -0.82   -0.13    0.20    0.54    1.20  5000 1.00
eta[2,14]                  0.22    0.01  0.50   -0.75   -0.12    0.22    0.56    1.22  5000 1.00
eta[2,15]                  0.48    0.01  0.51   -0.51    0.14    0.48    0.83    1.47  5000 1.00
eta[2,16]                  0.62    0.01  0.51   -0.38    0.28    0.62    0.96    1.59  5000 1.00
eta[2,17]                  0.10    0.01  0.50   -0.88   -0.25    0.09    0.43    1.11  5000 1.00
eta[2,18]                  0.13    0.01  0.51   -0.86   -0.22    0.14    0.48    1.13  5000 1.00
eta[2,19]                 -1.23    0.01  0.50   -2.21   -1.57   -1.24   -0.88   -0.25  5000 1.00
eta[2,20]                 -0.49    0.01  0.50   -1.46   -0.83   -0.49   -0.14    0.48  5000 1.00
eta[2,21]                 -0.53    0.01  0.51   -1.55   -0.86   -0.54   -0.19    0.46  5000 1.00
eta[2,22]                 -0.03    0.01  0.51   -1.01   -0.36   -0.03    0.31    0.96  5000 1.00
eta[2,23]                  0.00    0.01  0.50   -0.98   -0.35    0.00    0.33    1.00  5000 1.00
eta[2,24]                  1.79    0.01  0.50    0.81    1.46    1.79    2.13    2.77  5000 1.00
eta[2,25]                  0.37    0.01  0.51   -0.61    0.03    0.37    0.72    1.34  5000 1.00
eta[2,26]                  0.63    0.01  0.50   -0.36    0.29    0.63    0.97    1.59  5000 1.00
eta[2,27]                 -0.36    0.01  0.49   -1.31   -0.69   -0.36   -0.03    0.60  5000 1.00
eta[2,28]                 -0.81    0.01  0.49   -1.79   -1.15   -0.80   -0.47    0.13  5000 1.00
eta[2,29]                  0.54    0.01  0.50   -0.44    0.20    0.55    0.88    1.51  5000 1.00
eta[2,30]                 -0.52    0.01  0.50   -1.51   -0.85   -0.53   -0.19    0.47  5000 1.00
eta[2,31]                 -0.15    0.01  0.51   -1.16   -0.49   -0.14    0.21    0.84  5000 1.00
eta[2,32]                 -0.26    0.01  0.50   -1.24   -0.59   -0.26    0.08    0.70  5000 1.00
eta[3,1]                   0.49    0.01  0.69   -0.87    0.03    0.49    0.96    1.84  5000 1.00
eta[3,2]                  -0.32    0.01  0.71   -1.71   -0.81   -0.32    0.17    1.07  5000 1.00
eta[3,3]                  -0.62    0.01  0.71   -2.03   -1.09   -0.61   -0.16    0.75  5000 1.00
eta[3,4]                   0.05    0.01  0.71   -1.36   -0.42    0.06    0.52    1.45  5000 1.00
eta[3,5]                   0.18    0.01  0.70   -1.16   -0.31    0.17    0.65    1.53  5000 1.00
eta[3,6]                   0.64    0.01  0.68   -0.68    0.18    0.65    1.09    1.97  5000 1.00
eta[3,7]                  -0.80    0.01  0.70   -2.16   -1.28   -0.81   -0.32    0.58  5000 1.00
eta[3,8]                   0.08    0.01  0.71   -1.28   -0.42    0.08    0.57    1.44  5000 1.00
eta[3,9]                   0.06    0.01  0.70   -1.30   -0.41    0.06    0.52    1.43  5000 1.00
eta[3,10]                  0.38    0.01  0.71   -1.05   -0.08    0.39    0.84    1.80  5000 1.00
eta[3,11]                 -0.07    0.01  0.72   -1.48   -0.55   -0.06    0.42    1.28  5000 1.00
eta[3,12]                 -0.88    0.01  0.71   -2.27   -1.35   -0.88   -0.39    0.49  5000 1.00
eta[3,13]                  0.04    0.01  0.71   -1.37   -0.44    0.04    0.52    1.44  5000 1.00
eta[3,14]                  0.04    0.01  0.70   -1.32   -0.44    0.03    0.52    1.43  5000 1.00
eta[3,15]                 -0.10    0.01  0.70   -1.51   -0.58   -0.10    0.38    1.27  5000 1.00
eta[3,16]                 -0.12    0.01  0.70   -1.50   -0.57   -0.12    0.35    1.28  5000 1.00
eta[3,17]                  0.00    0.01  0.70   -1.36   -0.48   -0.01    0.45    1.35  5000 1.00
eta[3,18]                  0.48    0.01  0.73   -0.94   -0.03    0.48    0.99    1.92  5000 1.00
eta[3,19]                 -1.79    0.01  0.70   -3.15   -2.25   -1.77   -1.33   -0.39  5000 1.00
eta[3,20]                  0.37    0.01  0.68   -0.97   -0.09    0.37    0.83    1.73  5000 1.00
eta[3,21]                  0.07    0.01  0.70   -1.36   -0.41    0.08    0.54    1.42  5000 1.00
eta[3,22]                  0.15    0.01  0.71   -1.22   -0.34    0.13    0.64    1.54  5000 1.00
eta[3,23]                 -0.85    0.01  0.70   -2.22   -1.32   -0.84   -0.39    0.54  5000 1.00
eta[3,24]                  0.08    0.01  0.71   -1.32   -0.40    0.09    0.56    1.50  5000 1.00
eta[3,25]                 -0.48    0.01  0.71   -1.83   -0.98   -0.48    0.02    0.87  5000 1.00
eta[3,26]                  0.37    0.01  0.70   -1.00   -0.09    0.37    0.83    1.74  5000 1.00
eta[3,27]                  0.43    0.01  0.68   -0.91   -0.01    0.43    0.88    1.83  5000 1.00
eta[3,28]                 -0.37    0.01  0.69   -1.76   -0.83   -0.37    0.10    0.97  5000 1.00
eta[3,29]                  0.53    0.01  0.70   -0.87    0.05    0.53    1.01    1.88  5000 1.00
eta[3,30]                 -0.07    0.01  0.71   -1.46   -0.54   -0.07    0.41    1.30  5000 1.00
eta[3,31]                 -0.56    0.01  0.71   -1.95   -1.04   -0.55   -0.06    0.79  5000 1.00
 [ reached getOption("max.print") -- omitted 4538 rows ]

Samples were drawn using NUTS(diag_e) at Sat Dec  9 12:16:42 2017.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

check_div(fit_gp_h)

[1] "0 of 5000 iterations ended with a divergence (0%)"

check_treedepth(fit_gp_h)

[1] "0 of 5000 iterations saturated the maximum tree depth of 10 (0%)"

check_energy(fit_gp_h)

[1] "E-BFMI indicated no pathological behavior"

By increasing the adapt_delta and tree depth the model have no problems with divergences and does not hit the maximum trajectory threshold.

Posterior predictive checking

library(bayesplot)
sample = extract(fit_gp_h)
ppc_dens_overlay(y=data[,"FI"], yrep=sample$y_pred[,sample$country_indicators[1,] == match("FI",colnames(data))][,1:15])+ggtitle("Posterior predictive for Finland (hierarchical model)") + xlab("Debt (% of GDP)")

Hierarchical model follows the emprical data distribution quite well, better than the separate model for Finland. Variance at 40-45 (2008 financial crisis) is not as large as it is with the separate model. Moreover variance between different draws in general is lower than with separate model. As expected by increasing the complexity of our model by adding hierarchy makes posterior fit the observed phenomena better.

Comparing separate and hierarchical model

Model comparision can be achieved using psis-LOO.

library(loo)
ll_sep <- extract_log_lik(sample1$m)
l_sep <- loo(ll_sep)

Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.

ll_hier <- extract_log_lik(fit_gp_h)[,country_ind_is == match("FI",colnames(data))]
l_hier <- loo(ll_hier)

Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.

print(l_sep)

Computed from 5000 by 12 log-likelihood matrix

         Estimate  SE
elpd_loo    -30.8 4.2
p_loo         5.8 2.1
looic        61.6 8.4

Pareto k diagnostic values:
                         Count  Pct 
(-Inf, 0.5]   (good)     6     50.0%
 (0.5, 0.7]   (ok)       5     41.7%
   (0.7, 1]   (bad)      1      8.3%
   (1, Inf)   (very bad) 0      0.0%
See help('pareto-k-diagnostic') for details.

print(l_hier)

Computed from 5000 by 15 log-likelihood matrix

         Estimate  SE
elpd_loo    -30.7 0.9
p_loo         4.5 0.6
looic        61.5 1.9

Pareto k diagnostic values:
                         Count  Pct 
(-Inf, 0.5]   (good)      3    20.0%
 (0.5, 0.7]   (ok)       10    66.7%
   (0.7, 1]   (bad)       2    13.3%
   (1, Inf)   (very bad)  0     0.0%
See help('pareto-k-diagnostic') for details.

Pareto-k values are quite reasonable. Only a few values are above 0.7 for both models. Based on LOO values hierarchical model performs slightly better because elpd_loo value is higher. However we could still consider using more computationally heavy model comparison methods, because there are still a few values over 0.7, which makes use of psis-LOO for model comparison possibly questionable. Cross validation could be used instead of psis-LOO.

Problems and Future improvements

Finding a proper prior was a bit tricky, because we had only a few datapoints, at least for the separate model. If the prior is too vague, gaussian process overfits badly and if decay speed of the length scale is too slow, gaussian process smoothens the curve too much and makes the prediction too linear.

One option to increase the usefulness of the hierarchical model is to add more layers to the hierarchy. E.g we could group countries into subgroups so that nordic countries have their own “grouped” gaussian process on top of the global and country specfic gp model. There is definately more similarities between nordic countries than with e.g Finland and Italy.

Model comparison should be improved by comparing the models using cross validation. Sensitivity analysis should also be done for hierarchical model.

References

Appendix: Licenses

LS0tCnRpdGxlOiAiQmF5ZXNpYW4gZGF0YSBhbmFseXNpcyBwcm9qZWN0IgpvdXRwdXQ6CiAgaHRtbF9ub3RlYm9vazogZGVmYXVsdAogIHBkZl9kb2N1bWVudDogZGVmYXVsdAotLS0KYGBge3IsIGVjaG89RkFMU0V9CmtuaXRyOjpvcHRzX2NodW5rJHNldChlcnJvciA9IFRSVUUpCmBgYAoKCiMjIERhdGEgYW5kIHRoZSBwcm9ibGVtCgpUaGUgY2hvc2VuIGRhdGEgaXMgR292ZXJubWVudCBkZWJ0IGFzIHBlcmNlbnQgb2YgR0RQIChhdmFpbGFibGUgYXQ6IGh0dHBzOi8vd3d3Lmdvb2dsZS5jb20vcHVibGljZGF0YS9leHBsb3JlP2RzPWRzMjJhMzRrcmhxNXBfIyEgKS4gVGhlIGRhdGEgY29udGFpbnMgdGhlIHllYXJseSBkZWJ0IHByb2dyZXNzaW9uIG9mIEV1cm9wZWFuIGNvdW50cmllcyBvdXQgb2Ygd2hpY2ggYSBzdWJzZXQgb2YgY291bnRyaWVzIGlzIGFuYWx5emVkIGhlcmUuIFRoZSByYW5nZSBvZiBkYXRhIHdhcyBsaW1pdGVkIGJldHdlZW4geWVhcnMgMjAwMi0yMDE2IGZvciBjb25zaXN0ZW50IGRhdGEgYXZhaWxhYmlsaXR5LgoKVGhlIHByb2JsZW0gYWRkcmVzc2VkIGlzIHRoZSBwcmVkaWN0aW9uIG9mIHRoZSBkZWJ0IHByb2dyZXNzIGJhc2VkIG9uIGhpc3RvcnkgZGF0YSB1c2luZyBnYXVzc2lhbiBwcm9jZXNzZXMuCgoKYGBge3J9CmxpYnJhcnkoZ2dwbG90MikKbGlicmFyeShyZXNoYXBlMikKbGlicmFyeShnZ3RoZW1lcykKCmRhdGEgPC0gcmVhZC5jc3YoImRhdGEudHN2IixoZWFkZXI9VCxzZXA9Ilx0IixuYS5zdHJpbmdzID0gIjoiLHJvdy5uYW1lcyA9IDEpCgpkYXRhMiA8LSBkYXRhWywtYygxOjcpXQoKCmwgPC0gcm93bmFtZXMoZGF0YTIpCm4gPC0gZGF0YTIkbmFtZQoKZGF0YTIgPC0gYXMuZGF0YS5mcmFtZSh0KGRhdGEyWywtMV0pKQpjb2xuYW1lcyhkYXRhMikgPC0gbgpkYXRhMiRteWZhY3RvciA8LSBmYWN0b3Iocm93Lm5hbWVzKGRhdGEyKSkKY29sbmFtZXMoZGF0YTIpIDwtIGwKZGF0YTIgPC0gZGF0YTJbLC1jKDMzKV0KZGF0YTIkeWVhciA8LSBjKDIwMDM6MjAxNikKCmNuYW1lcyA8LSBjb2xuYW1lcyhkYXRhMikKZGZfbWVsdCA8LSBtZWx0KGRhdGEyWywgY25hbWVzXSwgaWQ9InllYXIiKQoKZ2cgPC0gZ2dwbG90KGRmX21lbHQsIGFlcyh4PXllYXIpKQpnZyArIGdlb21fbGluZShhZXMoeT12YWx1ZSwgY29sb3I9dmFyaWFibGUpLCBzaXplPTEpICsgc2NhbGVfY29sb3JfZGlzY3JldGUobmFtZT0iTGVnZW5kIikgKwogIHlsYWIoIiUgb2YgR0RQIikgKwogIGdndGl0bGUoIkdlbmVyYWwgZ292ZXJubWVudCBncm9zcyBkZWJ0ICUgb2YgR0RQIikgKwogIHRoZW1lX2lncmF5KCkKYGBgCgoKIyMgTW9kZWwgZGVzY3JpcHRpb24gKHNlcGFyYXRlIG1vZGVsIGZvciBlYWNoIGNvdW50cnkpCldlIGhhdmUgZGV2ZWxvcGVkIHR3byBtb2RlbHMgaGVyZS4gV2Ugc3RhcnQgd2l0aCBzZXBhcmF0ZSBtb2RlbHMgZm9yIGVhY2ggY291bnRyaWVzLiBJbiB0aGUgc2Vjb25kIHBhcnQgb2YgdGhpcyByZXBvcnQgd2UgZGV2ZWxvcCB0aGUgbW9kZWwgZnVydGhlciBhbmQgYWRkIGhpZXJhcmNoaWNhbCBjb21wb25lbnRzIGFzIHdlbGwuCgpUaGUgZGF0YSBpcyBhIHRpbWUtc2VyaWVzIGRhdGEgYnkgbmF0dXJlIHdoaWNoIG1ha2VzIEdhdXNzaWFuIFByb2Nlc3NlcyBhIG5hdHVyYWwgY2hvaWNlLiBUaGUgZ3Jvd2luZyB0cmVuZCBvZiB0aGUgZGF0YSBpcyBpbmNvcnBvcmF0ZWQgaW50byB0aGUgbW9kZWwgYXMgYSBsaW5lYXIgY29tcG9uZW50OgoKXGJlZ2lue2VxdWF0aW9ufQpcbXVfaSA9IFxtdV8wICsgXG11X2IgXGNkb3QgeF9pClxlbmR7ZXF1YXRpb259CgpUaGVuIHdlIGhhdmUgdGhlIGV4cG9uZW50aWF0ZWQgcXVhZHJhdGljIGtlcm5lbCBhcyBhIGNvdmFyaWFuY2U6CgpcYmVnaW57ZXF1YXRpb259Cktfe2ksan0gPSBcYWxwaGFeMiBcZXhwIFxCaWcoIC0gXGZyYWN7MX17MiBsXjJ9ICh4X2kgLSB4X2opXjIgXEJpZykKXGVuZHtlcXVhdGlvbn0KCkFuZCB0b2dldGhlciB0aGVzZSBjb21wb3NlIHRoZSBmdWxsIHNlcGFyYXRlIG1vZGVsIChpbiBhZGRpdGlvbiB0byB0aGUgcmFuZG9tIG5vaXNlIHNpZ21hKToKClxiZWdpbnthbGlnbn0KZiAmIFxzaW0gTXVsdGlOb3JtYWwoXG11LCBLKSBcXAp5X2kgJiBcc2ltIE5vcm1hbChmX2ksIFxzaWdtYSkKXGVuZHthbGlnbn0KCiMjIFByaW9yIGNob2ljZXMKVGhlcmUgYXJlIGluIHRvdGFsIGZvdXIgaW5pdGlhbCBwcmlvciBjaG9pY2VzIGluIHRoZSBtb2RlbC4gV2UgaGF2ZSBsaXN0ZWQgb3VyIGluaXRhbCBjaG9pY2VzIGZvciB0aGUgcHJpb3JzIGhlcmUsIGhvd2V2ZXIsIHRoZXNlIHdpbGwgYWRqdXN0ZWQgbGF0ZXIgZHVlIHRvIHNvbWUgY29udmVyZ2VuY2UgaXNzdWVzLiBUaGUgaW5pdGlhbCBwcmlvciBjaG9pY2VzIGFyZToKCiogJFxtdV8wIFxzaW0gSGFsZlxfTm9ybWFsKDAsIDEwMCkkCiogJFxtdV9iIFxzaW0gTm9ybWFsKDAsIDUpJAoqICRcYWxwaGEgXHNpbSBOb3JtYWwoMCwgMSkkCiogJGxlbmd0aFxfcGFyYW0gXHNpbSBXZWlidWxsKDEwLCAyKSQKCkluIGFkZGl0aW9uLCB0aGVyZSBhcmUgdHdvIHBhcmFtZXRlcnMgJGV0YSQgYW5kICRzaWdtYSQgZm9yIHdoaWNoIHN0YW5kYXJkIHByaW9yIGNob2ljZXMgd2VyZSB1c2VkLiAkRXRhJCBpcyBuZWVkZWQgZm9yIHRoZSBDaG9sZXNreSBkZWNvbXBvc2l0aW9uIHRvIGdlbmVyYXRlIHRoZSBjb3ZhcmlhbmNlIG1hdHJpeCwgJHNpZ21hJCBpcyB1c2VkIHRvIGluY29ycG9yYXRlIHRoZSByYW5kb20sIHVuZXhwbGFpbmFibGUgdmFyaWF0aW9uIGluIHRoZSBtb2RlbC4gQm90aCB1c2Ugc3RhbmRhcmQgcHJpb3IgY2hvaWNlcyAkTm9ybWFsKDAsIDEpJC4KCiMjIyAkXG11XzAgLCBcbXVfYiQKVGhlc2UgdHdvIHBhcmFtZXRlcnMgZm9ybSB0aGUgbGluZWFyIHRyZW5kIGZvciB0aGUgbW9kZWwuIFRoZSBpbnRlcmNlcHQgJG11XF8wJCBpcyByZXN0cmljdGVkIHRvIHBvc2l0aXZlIGFuZCBsZWZ0IHZhZ3VlLCAkbXVcXzAgXHNpbSBIYWxmXF9Ob3JtYWwoMCwgMTAwKSQsIHdoaWNoIGNvdmVycyB0aGUgcmFuZ2Ugb2YgZ292ZXJubWVudCBkZWJ0cyBpbiB5ZWFyIDIwMDIuIFRoZSBzbG9wZSBwYXJhbWV0ZXIgaXMgYWxzbyBmYWlybHkgdmFndWUgYnV0IHBvc3NpYmx5IG5vbi1wb3NpdGl2ZSBhcyB3ZWxsLCAkbXVcX2IgXHNpbSBOb3JtYWwoMCwgNSkkLgoKIyMjICRcYWxwaGEkCkFscGhhIGlzIHRoZSBzdHJlbmd0aCBwYXJhbWV0ZXIgZm9yIHRoZSBHYXVzc2lhbiBQcm9jZXNzLiBUaGUgdmFndWUgcHJpb3IgJGFscGhhIFxzaW0gTm9ybWFsKDAsMSkkIHdhcyB0aGUgaW5pdGlhbCBjaG9pY2UgYW5kIHNlZW1zIHRvIHByb3ZpZGUgc3VmZmljaWVudCBzdHJlbmd0aCBmb3IgdGhlIHZhcmlhdGlvbi4KCiMjIyBsZW5ndGhfcGFyYW0KTGVuZ3RoX3BhcmFtIGNvbnRyb2xzIHRoZSBzcGVlZCBhdCB3aGljaCB0aGUgY29ycmVsYXRpb24gZWZmZWN0IGRlY2F5cyBhcyBkaXN0YW5jZSBiZXR3ZWVuIGRhdGEgcG9pbnRzIGdyb3dzIChsb3dlciB2YWx1ZSA9IGZhc3RlciBkZWNheSkuIFNpbmNlIHRoZSBhbW91bnQgb2YgZGF0YSBwb2ludHMgaXMgbG93IGFuZCB0aGUgZGF0YSBiZWhhdmVzIGV4dHJlbWVseSBpcnJlZ3VsYXJseSwgdXNpbmcgYSB2YWd1ZSBwcmlvciBmb3IgdGhlIGxlbmd0aCBwYXJhbWV0ZXIgcmVzdWx0ZWQgaW4gbG93IHZhbHVlcyAoIDwgMS4wKSAgdmFsdWUgdGhhdCBvdmVyZml0IHRoZSBkYXRhLiBVc2luZyBwcmlvciBjaG9pY2UgJGxlbmd0aFxfcGFyYW0gXHNpbSBXZWlidWxsKDEwLCAyKSQsIHRoZSBtb2RlbCBwdXRzIGxlc3MgZW1waGFzaXMgb24gdGhlIGhpZ2ggZnJlcXVlbmN5IGZsdWN0dWF0aW9uIGFuZCBmaXQgYmVjb21lcyBzbW9vdGhlci4gV2UgdGhpbmsgdGhpcyBpcyByZWFzb25hYmxlLCBiZWNhdXNlIHRoZSBmaW5hbmNpYWwgY3Jpc2lzIGF0IDIwMDggaXMgaW1wb3NzaWJsZSB0byBleHBsYWluIGJhc2VkIG9ubHkgb24gdGhpcyBkYXRhc2V0LCBzbyB3ZSB3YW50IHRvIHByZXZlbnQgb3ZlcmZpdHRpbmcgYXMgbXVjaCBhcyBwb3NzaWJsZS4KCmBgYHtyLCBlY2hvPUZBTFNFLCBtZXNzYWdlPUZBTFNFLCB3YXJuaW5nPUZBTFNFfQp4IDwtIHNlcSgwLDUsMC4wMDEpCmcgPC0gZ2dwbG90KGRhdGEuZnJhbWUoeD14LHk9ZHdlaWJ1bGwoeCwxMCwyKSksYWVzKHg9eCx5PXkpKQpnICsgZ2VvbV9saW5lKCkrZ2d0aXRsZSgiV2VpYnVsbCgxMCwyKSIpCmBgYAoKCiMjIFN0YW4gKHNlcGFyYXRlIG1vZGVsKQpgYGB7ciwgZWNobz1GQUxTRX0Kd3JpdGVMaW5lcyhyZWFkTGluZXMoIi91LzE2L2FyYXNhbG8xL3VuaXgva291bHUvYmF5ZXMvcHJvamVjdC9kZWJ0X2dwX2xpbmVhcl93LnN0YW4iKSkKYGBgCgojIyBSIGNvZGUgKHNlcGFyYXRlIG1vZGVsKQoKVHJhaW5pbmcgb2YgdGhlIGluaXRpYWwgc2VwYXJhdGUgbW9kZWwgYW5kIHBsb3RzIGZvciBkaWZmZXJlbnQgY291bnRyaWVzLgpgYGB7ciwgbWVzc2FnZT1GQUxTRSwgd2FybmluZz1GQUxTRSwgaW5jbHVkZT1UUlVFLCByZXN1bHRzPUZBTFNFfQpsaWJyYXJ5KHJzdGFuKQpyc3Rhbl9vcHRpb25zKGF1dG9fd3JpdGUgPSBUUlVFKQpvcHRpb25zKG1jLmNvcmVzID0gcGFyYWxsZWw6OmRldGVjdENvcmVzKCkpCmRhdGEgPSByZWFkLmNzdignZGF0YS5jc3YnLCByb3cubmFtZXM9MSkKZGF0YSA9IGRhdGEuZnJhbWUodChkYXRhKSkKZGF0YVssJ3llYXInXSA9IHNlcSgyMDAyLCAyMDE2LCBieT0xKQoKZ3AgPC0gZnVuY3Rpb24oY291bnRyeV9pZCxjb3VudHJ5X25hbWUpIHsKCiAgY2MgPSBjb3VudHJ5X2lkCiAgeSA9IGRhdGFbLGNjXQogIHggPSAxOmxlbmd0aCh5KQogIAogIHNzID0gMwogIHlfaXMgPSB5WzE6KGxlbmd0aCh5KSAtIHNzKV0KICB4X2lzID0geFsxOihsZW5ndGgoeCkgLSBzcyldCiAgTl9pcyA9IGxlbmd0aCh4X2lzKQogIAogIHlfb3MgPSB0YWlsKHksIChzcyArIDEpKQogIHhfb3MgPSB0YWlsKHgsIChzcyArIDEpKQogIAogIHhfcHJlZCA9IHRhaWwoeCwgc3MpCiAgTl9wcmVkID0gbGVuZ3RoKHhfcHJlZCkKICAKICBncF9tb2RlbCA9ICdkZWJ0X2dwX2xpbmVhci5zdGFuJwogIAogIGZpdF9ncCA8LSBzdGFuKGZpbGU9Z3BfbW9kZWwsCiAgICAgICAgICAgICAgICAgZGF0YT1saXN0KHhfaXMgPSB4X2lzLAogICAgICAgICAgICAgICAgICAgICAgICAgICB5X2lzID0geV9pcywKICAgICAgICAgICAgICAgICAgICAgICAgICAgTl9pcyA9IE5faXMsCiAgICAgICAgICAgICAgICAgICAgICAgICAgIHhfcHJlZCA9IHhfcHJlZCwKICAgICAgICAgICAgICAgICAgICAgICAgICAgTl9wcmVkKSwKICAgICAgICAgICAgICAgICBpdGVyPTIwMCwKICAgICAgICAgICAgICAgICBjaGFpbnM9MykKICAKICBzYW1wbGUgPSBleHRyYWN0KGZpdF9ncCkKICBzYW1wbGVfcHJlZHMgPSBkYXRhLmZyYW1lKHNhbXBsZSR5X3ByZWQpCiAgICAKICBzYW1wbGVfcXVhbnRpbGVzIDwtIGFwcGx5KHNhbXBsZV9wcmVkcywyLGZ1bmN0aW9uKHgpe3F1YW50aWxlKHgsYygwLjA1LDAuNSwwLjk1KSl9KQogIAogIHNhbXBsZV9kYXRhIDwtIGRhdGEuZnJhbWUobG93ID0gc2FtcGxlX3F1YW50aWxlc1sxLF0sCiAgICAgICAgICAgICAgICAgICAgICAgICAgICBtZWQgPSBzYW1wbGVfcXVhbnRpbGVzWzIsXSwKICAgICAgICAgICAgICAgICAgICAgICAgICAgIGhpZ2ggPSBzYW1wbGVfcXVhbnRpbGVzWzMsXSwKICAgICAgICAgICAgICAgICAgICAgICAgICAgIHggPSB4KQogIHRyYWluaW5nX2RhdGEgPC0gZGF0YS5mcmFtZSh4PXhfaXMsIG9icz15X2lzKQogIHZhbGlkYXRpb25fZGF0YSA8LSBkYXRhLmZyYW1lKHg9eF9vcywgb2JzPXlfb3MpCiAgCiAgY29scyA8LSBjKCJHUCBtZWRpYW4iPSJibGFjayIsIkdQIDkwJSBwZXJjZW50aWxlIj0ibGlnaHRibHVlIiwiVHJ1ZSBkZWJ0Ij0icmVkIikKICBnZ3Bsb3Qoc2FtcGxlX2RhdGEsIGFlcyh4PXgsIHk9bWVkKSkgKwogICAgZ2VvbV9yaWJib24oYWVzKHg9eCx5bWluPWxvdyx5bWF4PWhpZ2gsY29sb3VyPSJHUCA5MCUgcGVyY2VudGlsZSIpLGFscGhhPTAuNixmaWxsPSJsaWdodGJsdWUiKSArCiAgICBnZW9tX2xpbmUoYWVzKHg9eCx5PW1lZCxjb2xvdXI9IkdQIG1lZGlhbiIpKSArCiAgICBnZW9tX3BvaW50KCkgKwogICAgZ2VvbV9saW5lKGRhdGE9dHJhaW5pbmdfZGF0YSwgYWVzKHg9eCwgeT1vYnMsY29sb3VyPSJUcnVlIGRlYnQiKSkgKwogICAgZ2VvbV9saW5lKGRhdGE9dmFsaWRhdGlvbl9kYXRhLCBhZXMoeD14LCB5PW9icyxjb2xvdXI9IlRydWUgZGVidCIpLCBsaW5ldHlwZT0iZGFzaGVkIikgKwogICAgeGxhYigneWVhcicpICsgCiAgICB5bGFiKCdEZWJ0ICglIG9mIEdEUCknKSArIAogICAgZ2d0aXRsZShzcHJpbnRmKCdHb3Zlcm5tZW50IGRlYnQsIGNvdW50cnkgPSBcJyVzXCcnLCBjb3VudHJ5X25hbWUpKSArCiAgICBzY2FsZV9jb2xvdXJfbWFudWFsKG5hbWU9IkNvbG91ciIsdmFsdWVzPWNvbHMpICsKICAgIHRoZW1lKHBsb3QudGl0bGUgPSBlbGVtZW50X3RleHQoaGp1c3QgPSAwLjUpKQp9CmdwKCdGSScsJ0ZpbmxhbmQnKQpncCgnU0UnLCdTd2VkZW4nKQpncCgnRUwnLCdHcmVlY2UnKQpgYGAKCiMjIENvbnZlcmdlbmNlIGRpYWdub3N0aWNzCgpBbmFseXNpcyBpcyBkb25lIHVzaW5nIGEgc2VwYXJhdGUgbW9kZWwgZ2VuZXJhdGVkIGZvciBGaW5sYW5kLgoKYGBge3IsIG1lc3NhZ2U9RkFMU0UsIHdhcm5pbmc9RkFMU0UsIGluY2x1ZGU9VFJVRSwgcmVzdWx0cz1GQUxTRX0KY2MgPSAnRkknCnkgPSBkYXRhWyxjY10KeCA9IDE6bGVuZ3RoKHkpCgpzcyA9IDMKeV9pcyA9IHlbMToobGVuZ3RoKHkpIC0gc3MpXQp4X2lzID0geFsxOihsZW5ndGgoeCkgLSBzcyldCk5faXMgPSBsZW5ndGgoeF9pcykKCnlfb3MgPSB0YWlsKHksIChzcyArIDEpKQp4X29zID0gdGFpbCh4LCAoc3MgKyAxKSkKCnhfcHJlZCA9IHRhaWwoeCwgc3MpCk5fcHJlZCA9IGxlbmd0aCh4X3ByZWQpCgpncF9tb2RlbCA9ICdkZWJ0X2dwX2xpbmVhcl93LnN0YW4nCgpmaXRfZ3AgPC0gc3RhbihmaWxlPWdwX21vZGVsLAogICAgICAgICAgICAgICBkYXRhPWxpc3QoeF9pcyA9IHhfaXMsCiAgICAgICAgICAgICAgICAgICAgICAgICB5X2lzID0geV9pcywKICAgICAgICAgICAgICAgICAgICAgICAgIE5faXMgPSBOX2lzLAogICAgICAgICAgICAgICAgICAgICAgICAgeF9wcmVkID0geF9wcmVkLAogICAgICAgICAgICAgICAgICAgICAgICAgTl9wcmVkKSwKICAgICAgICAgICAgICAgaXRlcj0xMDAwLAogICAgICAgICAgICAgICBjaGFpbnM9MTAsY29udHJvbD1saXN0KGFkYXB0X2RlbHRhPTAuOTksbWF4X3RyZWVkZXB0aD0xNSkpCmBgYAojIyMjICRcaGF0e1J9JCBkaWFnbm9zdGljcwpgYGB7cn0KZml0X2dwCmBgYAoKCkFsbCAkXGhhdHtSfSQgdmFsdWVzIGFyZSBiZWxvdyAxLjEgc28gb3VyIE1DTUMgc2ltdWxhdGlvbnMgY29udmVyZ2VuY2VzIHdlbGwuIAoKCiMjIyBUcmVlIERlcHRoCgpgYGB7cn0Kc291cmNlKCJzdGFuX3V0aWxpdHkuUiIpCmNoZWNrX3RyZWVkZXB0aChmaXRfZ3ApCmBgYApNb250ZSBDYXJsbyBzaW11bGF0aW9ucyBkaWQgbm90IGhpdCB0aGUgbWF4aW11bSB0cmFqZWN0b3J5IGxlbmd0aCwgc28gd2UgZG9uJ3QgbmVlZCB0byBpbmNyZWFzZSB0aGUgdHJlZSBkZXB0aC4KCiMjIyBFLUJGTUkKCmBgYHtyfQpjaGVja19lbmVyZ3koZml0X2dwKQpgYGAKCiMjIyBEaXZlcmdlbmNlCgpgYGB7cn0KY2hlY2tfZGl2KGZpdF9ncCkKYGBgCgpEaXZlcmdlbmNlcyBkb2VzIG5vdCB2YW5pc2ggY29tcGxldGVseSBldmVuIHRob3VnaCBhZGFwdF9kZWx0YSB2YWx1ZSBpcyAwLjk5OS4gTGV0J3MgY2hlY2sgdGhlIHBvc3RlcmlvciBvZiBsZW5ndGhfcGFyYW0gKGxlbmd0aC1zY2FsZSkgYW5kIHNpZ21hLCB3aGljaCBpcyB0aGUgdW5leHBsZWluYWJsZSB2YXJpYW5jZSBpbiB0aGUgZGF0YSwgdG8gc2VlIGlmIGRpdmVyZ2VuY2VzIGFyZSBmYWxzZSBwb3NpdGl2ZXMgb3IgY2F1c2VkIGJ5IGFuIGFjdHVhbCBwcm9ibGVtLgoKYGBge3J9CmNfZGFyayA8LSBjKCIjOEYyNzI3ODAiKQpncmVlbiA8LSBjKCIjMDBGRjAwODAiKQoKcGFydGl0aW9uIDwtIHBhcnRpdGlvbl9kaXYoZml0X2dwKQpkaXZfcGFyYW1zIDwtIHBhcnRpdGlvbltbMV1dCm5vbmRpdl9wYXJhbXMgPC0gcGFydGl0aW9uW1syXV0KCnBhcihtYXIgPSBjKDQsIDQsIDAuNSwgMC41KSkKcGxvdChub25kaXZfcGFyYW1zJGxlbmd0aF9wYXJhbSwgbm9uZGl2X3BhcmFtcyRzaWdtYSwKICAgICBjb2w9Y19kYXJrLCBwY2g9MTYsIGNleD0wLjgsIHhsYWI9Imxlbmd0aF9wYXJhbSIsIHlsYWI9InNpZ21hIiwKICAgICB4bGltPWMoMCwgNCksIHlsaW09YygwLDcpKQpwb2ludHMoZGl2X3BhcmFtcyRsZW5ndGhfcGFyYW0sIGRpdl9wYXJhbXMkc2lnbWEsCiAgICAgICBjb2w9Z3JlZW4sIHBjaD0xNiwgY2V4PTAuOCkKYGBgCgpBcyBjYW4gYmUgc2VlbiBmcm9tIHRoZSBwbG90LCBhbGwgdGhlIGRpdmVyZ2VuY2VzIGFyZSBjbHVzdGVyZWQgdG8gdGhlIGJvdHRvbSBsZWZ0IHdoZXJlIHRoZSBwb3N0ZXJpb3IgYWJydXB0cyBxdWlja2x5LiBTbyB0aGVyZSBpcyBpbmRlZWQgcHJvYmxlbSB3aXRoIG91ciBtb2RlbC4gTW9udGUgQ2FybG8gaXMgbm90IGFibGUgdG8gZXhwbG9yZSB0aGF0IGFyZWEgcHJvcGVybHkgc28gb3VyIHBvc3RlcmlvciBpcyBzbGlnaHRseSBiaWFzZWQuCgojIyMgRml4aW5nIHRoZSBtb2RlbApPbmUgaXNzdWUgd2Ugc2VlIGluIHRoZSBwb3N0ZXJpb3IgcGxvdCBhYm92ZSBpcyB0aGF0IHNvbWUgc2ltdWxhdGlvbnMgb2YgdGhlIGxlbmd0aCBwYXJhbWV0ZXIgaW4gb3VyIGV4cG9uZW50aWFsIGdhdXNzaWFuIGtlcm5lbCAgZmFsbCBiZWxvdyAxLiBUaGlzIGlzIGEgcG9zc2libGUgY2F1c2Ugb2Ygb3ZlcmZpdCB0aGF0IG1pZ2h0IHJlc3VsdCBpbiBkaXZlcmdlbmNlcywgc2luY2Ugb3VyIGRhdGEgY29udGFpbnMgb25seSB5ZWFybHkgb2JzZXJ2YXRpb25zLgoKV2UgdHJ5IHRvIGZpeCBvdXIgbW9kZWwgYnkgYWRqdXN0aW5nIG91ciBwcmlvciBmb3IgbGVuZ3RoX3BhcmFtIChsZW5ndGgtc2NhbGUpLiBBdCB0aGUgbW9tZW50IHdlIGhhdmUgJFdlaWJ1bGwoMTAsMikkLCB3aGljaCBtaWdodCBnaXZlIHVzIHZhbHVlcyBiZXR3ZWVuIDAgYW5kIDEuIFdlIGNhbiBmaXggdGhpcyBieSBtb3Zpbmcgb3VyIHByaW9yIGZ1cnRoZXIgYXdheSBmcm9tIDEgYW5kIGNsaXBwaW5nIHRoZSBwcmlvciBzbyB0aGF0IGl0IGNhbiBvbmx5IGdldCB2YWx1ZXMgaGlnaGVyIHRoYW4gMS4gT3VyIHByaW9yIGZvciBsZW5ndGgtc2NhbGUgaXMgJFdlaWJ1bGwoMTUsMykkIHdpdGggbWluaW11bSBvZiAxLiBBZnRlciBhZGp1c3RpbmcgdGhlIHByaW9yLCB0aGUgY29udmVyZ2VuY2UgYW5hbHl0aWNzIGFyZSBwZXJmb3JtZWQgYWdhaW46CgpgYGB7ciwgaW5jbHVkZT1GQUxTRSwgZWNobz1GQUxTRX0KY2MgPSAnRkknCnkgPSBkYXRhWyxjY10KeCA9IDE6bGVuZ3RoKHkpCgpzcyA9IDMKeV9pcyA9IHlbMToobGVuZ3RoKHkpIC0gc3MpXQp4X2lzID0geFsxOihsZW5ndGgoeCkgLSBzcyldCk5faXMgPSBsZW5ndGgoeF9pcykKCnlfb3MgPSB0YWlsKHksIChzcyArIDEpKQp4X29zID0gdGFpbCh4LCAoc3MgKyAxKSkKCnhfcHJlZCA9IHRhaWwoeCwgc3MpCk5fcHJlZCA9IGxlbmd0aCh4X3ByZWQpCgpncF9tb2RlbCA9ICdkZWJ0X2dwX2xpbmVhci5zdGFuJwoKZml0X2dwIDwtIHN0YW4oZmlsZT1ncF9tb2RlbCwKICAgICAgICAgICAgICAgZGF0YT1saXN0KHhfaXMgPSB4X2lzLAogICAgICAgICAgICAgICAgICAgICAgICAgeV9pcyA9IHlfaXMsCiAgICAgICAgICAgICAgICAgICAgICAgICBOX2lzID0gTl9pcywKICAgICAgICAgICAgICAgICAgICAgICAgIHhfcHJlZCA9IHhfcHJlZCwKICAgICAgICAgICAgICAgICAgICAgICAgIE5fcHJlZCksCiAgICAgICAgICAgICAgIGl0ZXI9MTAwMCwKICAgICAgICAgICAgICAgY2hhaW5zPTEwLGNvbnRyb2w9bGlzdChhZGFwdF9kZWx0YT0wLjk5LG1heF90cmVlZGVwdGg9MTUpKQpgYGAKCgpgYGB7cn0KcHJpbnQoZml0X2dwKQpjaGVja190cmVlZGVwdGgoZml0X2dwKQpjaGVja19lbmVyZ3koZml0X2dwKQpjaGVja19kaXYoZml0X2dwKQpgYGAKQ2hhbmdpbmcgdGhlIHByaW9yIG9mIHRoZSBsZW5ndGggcGFyYW1ldGVyIHNlZW1zIHRvIGZpeCBvdXIgZGl2ZXJnZW5jZSBpc3N1ZXMuCgojIyMgVmlzdWFsaXphdGlvbgoKYGBge3IsIG1lc3NhZ2U9RkFMU0UsIHdhcm5pbmc9RkFMU0UsIGluY2x1ZGU9VFJVRSwgcmVzdWx0cz1GQUxTRX0Kc2FtcGxlID0gZXh0cmFjdChmaXRfZ3ApCnNhbXBsZV9wcmVkcyA9IGRhdGEuZnJhbWUoc2FtcGxlJHlfcHJlZCkKICAKc2FtcGxlX3F1YW50aWxlcyA8LSBhcHBseShzYW1wbGVfcHJlZHMsMixmdW5jdGlvbih4KXtxdWFudGlsZSh4LGMoMC4wNSwwLjUsMC45NSkpfSkKCnNhbXBsZV9kYXRhIDwtIGRhdGEuZnJhbWUobG93ID0gc2FtcGxlX3F1YW50aWxlc1sxLF0sCiAgICAgICAgICAgICAgICAgICAgICAgICAgbWVkID0gc2FtcGxlX3F1YW50aWxlc1syLF0sCiAgICAgICAgICAgICAgICAgICAgICAgICAgaGlnaCA9IHNhbXBsZV9xdWFudGlsZXNbMyxdLAogICAgICAgICAgICAgICAgICAgICAgICAgIHggPSB4KQp0cmFpbmluZ19kYXRhIDwtIGRhdGEuZnJhbWUoeD14X2lzLCBvYnM9eV9pcykKdmFsaWRhdGlvbl9kYXRhIDwtIGRhdGEuZnJhbWUoeD14X29zLCBvYnM9eV9vcykKCmNvbHMgPC0gYygiR1AgbWVkaWFuIj0iYmxhY2siLCJHUCA5MCUgcGVyY2VudGlsZSI9ImxpZ2h0Ymx1ZSIsIlRydWUgZGVidCI9InJlZCIpCmdncGxvdChzYW1wbGVfZGF0YSwgYWVzKHg9eCwgeT1tZWQpKSArCiAgZ2VvbV9yaWJib24oYWVzKHg9eCx5bWluPWxvdyx5bWF4PWhpZ2gsY29sb3VyPSJHUCA5MCUgcGVyY2VudGlsZSIpLGFscGhhPTAuNixmaWxsPSJsaWdodGJsdWUiKSArCiAgZ2VvbV9saW5lKGFlcyh4PXgseT1tZWQsY29sb3VyPSJHUCBtZWRpYW4iKSkgKwogIGdlb21fcG9pbnQoKSArCiAgZ2VvbV9saW5lKGRhdGE9dHJhaW5pbmdfZGF0YSwgYWVzKHg9eCwgeT1vYnMsY29sb3VyPSJUcnVlIGRlYnQiKSkgKwogIGdlb21fbGluZShkYXRhPXZhbGlkYXRpb25fZGF0YSwgYWVzKHg9eCwgeT1vYnMsY29sb3VyPSJUcnVlIGRlYnQiKSwgbGluZXR5cGU9ImRhc2hlZCIpICsKICB4bGFiKCd5ZWFyJykgKyAKICB5bGFiKCdEZWJ0ICglIG9mIEdEUCknKSArIAogIGdndGl0bGUoc3ByaW50ZignR292ZXJubWVudCBkZWJ0LCBjb3VudHJ5ID0gRmlubGFuZCcpKSArCiAgc2NhbGVfY29sb3VyX21hbnVhbChuYW1lPSJDb2xvdXIiLHZhbHVlcz1jb2xzKSArCiAgdGhlbWUocGxvdC50aXRsZSA9IGVsZW1lbnRfdGV4dChoanVzdCA9IDAuNSkpCmBgYAoKVGhlIGVtcGlyaWNhbCBkYXRhIHN0YXlzIHdpdGhpbmcgdGhlIDkwJSBjb25maWRlbmNlIGludGVydmFsIGluIGFsbCBjYXNlcyBhcGFydCBmcm9tIHRoZSBZZWFyIDcgKGNvcnJlc3BvbmRpbmcgdG8geWVhciAyMDA4KS4gVGhpcyBpcyB0aGUgeWVhciBvZiBiYW5rIGNyaXNpcyBhbmQgaXMgdGhlcmVmb3JlIGFuIGFub21hbHkgdG8gdGhlIHByb2Nlc3MgYXMgZG9lcyBub3QgZm9sbG93IHRoZSBuYXR1cmFsIHByb2dyZXNzaW9uIGV4cGxhaW5lZCBieSB0aGUgbW9kZWwuCgpDaGFuZ2luZyB0aGUgcHJpb3IgZm9yIHRoZSBsZW5ndGggcGFyYW0gc21vb3RoZW5zIHRoZSBtZWRpYW4gY29tcGFyZWQgdG8gdGhlIG9yaWdpbmFsIHByaW9yIGNob2ljZS4gVGhpcyBpcyBuYXR1cmFsLCBhcyBpbmNyZWFzaW5nIHRoZSBsZW5ndGggcGFyYW1ldGVyIGluY3JlYXNlcyB0aW1lIHNjYWxlIGF0IHdoaWNoIHRoZSBjb3JyZWxhdGlvbiBkZWNheXMuIEluIG90aGVyIHdvcmRzLCB0aGUgZWZmZWN0IG9mIHRoZSB0aW1lIHN0ZXBzIGZ1cnRoZXIgYXdheSBjb250cmlidXRlIG1vcmUgdG8gdGhlIHZhbHVlLCB0aHVzIHJlZHVjaW5nIHRoZSBmcmVxdWVuY3kgb2YgdGhlIHNpZ25hbC4KClRoZSBvdXQgb2Ygc2FtcGxlIHByZWRpY3Rpb25zIChkYXNoZWQgcmVkIGxpbmUpIGZhbGwgYWxzbyB3aXRoaW4gdGhlIGNvZmlkZW5jZSBpbnRlcnZhbC4gSG93ZXZlciwgdGhlIGRhdGFzZXQgaXMgc21hbGwgYW5kIHRoZXJlZm9yZSBtYWtpbmcgYWNjdXJhdGUgcHJlZGljdGlvbnMgaXMgbm90IGEgdHJpdmlhbCB0YXNrIHdpdGhvdXQgYWRkaXRpb25hbCBrbm93bGVkZ2UuCgoKIyMgUG9zdGVyaW9yIHByZWRpY3RpdmUgY2hlY2tpbmcKCmBgYHtyfQpsaWJyYXJ5KGJheWVzcGxvdCkKc2FtcGxlID0gZXh0cmFjdChmaXRfZ3ApCnBwY19kZW5zX292ZXJsYXkoeT15X2lzLCB5cmVwPXNhbXBsZSR5X3ByZWRbLDE6KGxlbmd0aCh5KSAtIHNzKV0pK2dndGl0bGUoIlBvc3RlcmlvciBwcmVkaWN0aXZlIGZvciBGaW5sYW5kIikgKyB4bGFiKCJEZWJ0ICglIG9mIEdEUCkiKQpgYGAKClNpbXVsYXRlZCBkcmF3cyBmcm9tIHRoZSBwb3N0ZXJpb3IgZGlzdHJpYnV0aW9uIGZvbGxvdyB0aGUgZW1waXJpY2FsIGRhdGEgZGlzdHJpYnV0aW9uIHJvdWdobHkuIEh1Z2UgdmFyaWFuY2UgaW4gZHJhd3MgYXQgNDAtNDUgaXMgZHVlIHRoZSBmYWN0IHRoYXQgb3VyIGluZm9ybWF0aXZlIHByaW9yICRXZWlidWxsKDE1LDMpJCBzbW9vdGhlbnMgdGhlIGRpcCAoMjAwOCBmaW5hbmNpYWwgY3Jpc2lzKSBxdWl0ZSBoZWF2aWx5IHNvIG91ciBwb3N0ZXJpb3IgZ2l2ZXMgbW9yZSBwcm9iYWJsaXR5IG1hc3MgYXJvdW5kIDQwLTQ1IGNvbXBhcmVkIHRvIHRoZSBlbXBpcmljYWwgZGF0YS4KCgojIyBTZW5zaXRpdml0eSBhbmFseXNpcwoKV2UgYXJlIGZvY3VzaW5nIG91ciBzZW5zaXRpdml0eSBhbmFseXNpcyBvbiBkaWZmZXJlbnQgcHJpb3JzIGZvciBsZW5ndGhfcGFyYW0gd2hpY2ggaXMgdGhlIGxlbmd0aC1zY2FsZSBwYXJhbWV0ZXIgZm9yIGdhdXNzaWFuIHByb2Nlc3MuIElkZWFsbHkgdGhlIHNlbnNpdGl2aXR5IGFuYWx5c2lzIHNob3VsZCBiZSBwZXJmb3JtZWQgZm9yIHRoZSBvdGhlciBwcmlvcnMgYXMgd2VsbCwgaG93ZXZlciwgZHVlIHRvIHRpbWUgY29uc3RyYWludHMgd2UgaGF2ZSBjaG9zZSB0aGUgdGhlIGxlbmd0aCBwYXJhbWV0ZXIgYXMgY2hhbmdpbmcgaXQgaXMgZXhwZWN0ZWQgdG8gaGF2ZSB0aGUgbW9zdCBub3RhYmxlIGVmZmVjdC4KCiMjIyBMZW5ndGhfcGFyYW0KCiogVmFndWUgcHJpb3I6ICRsZW5ndGhcX3BhcmFtIFxzaW0gXEdhbW1hKDQsNCkkCiogTm9uLWluZm9ybWF0aXZlIHByaW9yOiAkbGVuZ3RoXF9wYXJhbSBcc2ltIFVuaWZvcm0oMCwxMCkkCiogSW5mb3JtYXRpdmUgcHJpb3I6ICRsZW5ndGhcX3BhcmFtIFxzaW0gV2VpYnVsbCgxNSwgMykkCmBgYHtyLCBtZXNzYWdlPUZBTFNFLCB3YXJuaW5nPUZBTFNFLCBpbmNsdWRlPVRSVUUsIHJlc3VsdHM9RkFMU0V9CgpjYyA9ICdGSScKeSA9IGRhdGFbLGNjXQp4ID0gMTpsZW5ndGgoeSkKCnNzID0gMwp5X2lzID0geVsxOihsZW5ndGgoeSkgLSBzcyldCnhfaXMgPSB4WzE6KGxlbmd0aCh4KSAtIHNzKV0KTl9pcyA9IGxlbmd0aCh4X2lzKQoKeV9vcyA9IHRhaWwoeSwgKHNzICsgMSkpCnhfb3MgPSB0YWlsKHgsIChzcyArIDEpKQoKeF9wcmVkID0gdGFpbCh4LCBzcykKTl9wcmVkID0gbGVuZ3RoKHhfcHJlZCkKCnByaW9yX2NoZWNrIDwtIGZ1bmN0aW9uKG1vZGVsX25hbWUpIHsKICBmaXRfZ3BfdG1wIDwtIHN0YW4oZmlsZT1tb2RlbF9uYW1lLAogICAgICAgICAgICAgICBkYXRhPWxpc3QoeF9pcyA9IHhfaXMsCiAgICAgICAgICAgICAgICAgICAgICAgICB5X2lzID0geV9pcywKICAgICAgICAgICAgICAgICAgICAgICAgIE5faXMgPSBOX2lzLAogICAgICAgICAgICAgICAgICAgICAgICAgeF9wcmVkID0geF9wcmVkLAogICAgICAgICAgICAgICAgICAgICAgICAgTl9wcmVkKSwKICAgICAgICAgICAgICAgaXRlcj0xMDAwLAogICAgICAgICAgICAgICBjaGFpbnM9MTAsY29udHJvbD1saXN0KGFkYXB0X2RlbHRhPTAuOTksbWF4X3RyZWVkZXB0aD0xNSkpCiAgc2FtcGxlIDwtIGV4dHJhY3QoZml0X2dwX3RtcCkKICBzYW1wbGVfcHJlZHMgPC0gZGF0YS5mcmFtZShzYW1wbGUkeV9wcmVkKQogIAogIHNhbXBsZV9xdWFudGlsZXMgPC0gYXBwbHkoc2FtcGxlX3ByZWRzLDIsZnVuY3Rpb24oeCl7cXVhbnRpbGUoeCxjKDAuMDUsMC41LDAuOTUpKX0pCiAgCiAgc2FtcGxlX2RhdGEgPC0gZGF0YS5mcmFtZShsb3cgPSBzYW1wbGVfcXVhbnRpbGVzWzEsXSwKICAgICAgICAgICAgICAgICAgICAgICAgICBtZWQgPSBzYW1wbGVfcXVhbnRpbGVzWzIsXSwKICAgICAgICAgICAgICAgICAgICAgICAgICBoaWdoID0gc2FtcGxlX3F1YW50aWxlc1szLF0sCiAgICAgICAgICAgICAgICAgICAgICAgICAgeCA9IHgpCiAgdHJhaW5pbmdfZGF0YSA8LSBkYXRhLmZyYW1lKHg9eF9pcywgb2JzPXlfaXMpCiAgdmFsaWRhdGlvbl9kYXRhIDwtIGRhdGEuZnJhbWUoeD14X29zLCBvYnM9eV9vcykKCiAgcmV0dXJuKGxpc3Qocz1zYW1wbGVfZGF0YSwgdD10cmFpbmluZ19kYXRhLCB2PXZhbGlkYXRpb25fZGF0YSxtPWZpdF9ncF90bXApKQp9CgpzYW1wbGUxIDwtIHByaW9yX2NoZWNrKCdkZWJ0X2dwX2xpbmVhci5zdGFuJykKc2FtcGxlMiA8LSBwcmlvcl9jaGVjaygnZGVidF9ncF9saW5lYXIyLnN0YW4nKSAKc2FtcGxlMyA8LSBwcmlvcl9jaGVjaygnZGVidF9ncF9saW5lYXIzLnN0YW4nKQoKI2NvbHMgPC0gYygiV2VpYnVsbCgxMCwyKSI9ImJsYWNrIiwiR2FtbWEoNCw0KSI9ImdyZWVuIiwiVW5pZm9ybSBwcmlvciI9ImJsdWUiLCJUcnVlIGRlYnQiPSJyZWQiKQpnZ3Bsb3Qoc2FtcGxlMSRzLCBhZXMoeD14LCB5PW1lZCkpICsKICBnZW9tX2xpbmUoYWVzKHg9eCx5PW1lZCxjb2xvdXI9IldlaWJ1bGwoMTUsMykiKSkgKwogIGdlb21fbGluZShkYXRhPXNhbXBsZTIkcywgYWVzKHg9eCx5PW1lZCxjb2xvdXI9IkdhbW1hKDQsNCkiKSkgKyAKICBnZW9tX2xpbmUoZGF0YT1zYW1wbGUzJHMsIGFlcyh4PXgseT1tZWQsY29sb3VyPSJVbmlmb3JtIHByaW9yIikpICsgCiAgZ2VvbV9saW5lKGRhdGE9c2FtcGxlMSR0LCBhZXMoeD14LCB5PW9icyxjb2xvdXI9IlRydWUgZGVidCIpKSArCiAgZ2VvbV9saW5lKGRhdGE9c2FtcGxlMSR2LCBhZXMoeD14LCB5PW9icyxjb2xvdXI9IlRydWUgZGVidCIpLCBsaW5ldHlwZT0iZGFzaGVkIikgKwogIHhsYWIoJ3llYXInKSArIAogIHlsYWIoJ0RlYnQgKCUgb2YgR0RQKScpICsgCiAgZ2d0aXRsZShzcHJpbnRmKCdHb3Zlcm5tZW50IGRlYnQsIGNvdW50cnkgPSBGaW5sYW5kJykpICsKICAjc2NhbGVfY29sb3VyX21hbnVhbChuYW1lPSJDb2xvdXIiLHZhbHVlcz1jb2xzKSArCiAgdGhlbWUocGxvdC50aXRsZSA9IGVsZW1lbnRfdGV4dChoanVzdCA9IDAuNSkpCgpgYGAKCkFzIGNhbiBiZSBzZWVuIGZyb20gdGhlIHBsb3QsIHByaW9yIGluZm9ybWF0aW9uIGFib3V0IGxlbmd0aC1zY2FsZSBoYXMgYW4gbm90aWNhYmxlIGVmZmVjdCBvbiB0aGUgb3V0Y29tZS4gT3VyIGNob2ljZSBmb3IgaW5mb3JtYXRpdmUgcHJpb3Igc21vb3RoZW5zIHRoZSBjdXJ2ZSB0aGUgbW9zdC4gVXNpbmcgbm9uLWluZm9ybWF0aXZlIGFuZCB2YWd1ZSBwcmlvciBvdmVyZml0cyBtb2RlbCBxdWl0ZSBtdWNoLiBCYXNlZCBvbmx5IG9uIHZpc3VhbCBhbmFseXNpcyBvdXIgaW5mb3JtYXRpdmUgcHJpb3Igc2VlbXMgdGhlIG1vc3QgcmVhc29uYWJsZSwgYmVjYXVzZSB0aGUgZGF0YSBjb250YWluIGRhdGFwb2ludHMgdGhhdCBhcmUgYmFzaWNhbGx5IGltcG9zc2libGUgdG8gZXhwbGFpbiBiYXNlZCBvbmx5IG9uIHRoaXMgZGF0YSAoZmluYW5jaWFsIGNyaXNpcykuIEhvd2V2ZXIsIG91ciBtb2RlbCBpcyBzZW5zaXRpdmUgdG8gdGhlIGNob2ljZSBvZiBwcmlvciBmb3IgbGVuZ3RoIHNjYWxlLgoKIyMgTW9kZWwgY29tcGFyaXNvbiBmb3IgZGlmZmVyZW50IHByaW9ycyAoTE9PKQoKIyMjIEluZm9ybWF0aXZlIHByaW9yCmBgYHtyfQpsaWJyYXJ5KGxvbykKbGwgPC0gZXh0cmFjdF9sb2dfbGlrKHNhbXBsZTEkbSkKbCA8LSBsb28obGwpCmwKYGBgCgojIyMgdmFndWUgcHJpb3IKYGBge3J9CmxsIDwtIGV4dHJhY3RfbG9nX2xpayhzYW1wbGUyJG0pCmwgPC0gbG9vKGxsKQpsCmBgYAoKIyMjIG5vbi1pbmZvcm1hdGl2ZSBwcmlvcgpgYGB7cn0KbGwgPC0gZXh0cmFjdF9sb2dfbGlrKHNhbXBsZTMkbSkKbCA8LSBsb28obGwpCmwKYGBgCkNvbXBhcmluZyB0aGUgbW9kZWxzIHdpdGggZGlmZmVyZW50IHByaW9yIHVzaW5nIHBzaXMtTE9PIGlzIG5vdCBwb3NzaWJsZSwgYmVjYXVzZSBtb2RlbHMgZm9yIG5vbiBpbmZvcm1hdGl2ZSBwcmlvciBhbmQgdmFndWUgcHJpb3IgZ2V0IFBhcmV0byBrIHZhbHVlcyA+MS4wLiBUaGVyZWZvcmUgd2Ugc2hvdWxkIHVzZSBtb3JlIGNvbXB1dGF0aW9uYWxseSBleHBlbmlzaXZlIG1ldGhvZHMsIHN1Y2ggYXMgc3RhbmRhcmQgY3Jvc3MtdmFsaWRhdGlvbiwgdG8gZXZhbHV0ZSB0aGUgcGVyZm9ybWFuY2Ugb2Ygb3VyIG1vZGVscy4gSG93ZXZlciBwYXJldG8tayB2YWx1ZXMgZm9yIG91ciBpbmZvcm1hdGl2ZSBwcmlvciBtb2RlbCBhcmUgcXVpdGUgcmVhc29uYmxlOiBvbmx5IDIgdmFsdWVzIGFyZSBhYm92ZSAwLjcgc28gaXQgaXMgcG9zc2libGUgdG8gY29tcGFyZSB0aGF0IG1vZGVsIHRvIG91ciBoaWVyYXJjaGljYWwgbW9kZWwuCgoKIyMgU2Vjb25kIG1vZGVsIChoaWVyYXJjaGljYWwpCgpPdXIgc2Vjb25kIG1vZGVsIHRyaWVzIHRvIGluY29ycG9yYXRlIG1vcmUgaW5mb3JtYXRpb24gdG8gZWFjaCBtb2RlbCBieSBhZGRpbmcKaGllcmFyY2hpY2FsIHN0cnVjdHVyZS4gVGhlIGlkZWEgaXMgdG8gYWRkIGdsb2JhbCB0cmVuZCBhbmQgZ2F1c3NpYW4gcHJvY2VzcyB0bwp0aGUgc2VwYXJhdGUgbW9kZWwgaW4gdGhlIHByZXZpb3VzIHNlY3Rpb24uIFRoZSBnbG9iYWwgY29tcG9uZW50cyB3aWxsIGJlCmFmZmVjdGVkIGJ5IGFsbCB0aGUgRVUgY291bnRyaWVzIGluIG91ciBkYXRhc2V0LgoKU2ltaWxhcmx5IHRvIHRoZSBzZXBhcmF0ZSBtb2RlbCwgdGhlIGNvbXBvbmVudHMgYXJlIGZvcm1lZCBhczoKClxiZWdpbntlcXVhdGlvbn0gXG11X2kgPSBcbXVfMCArIFxtdV9iIFxjZG90IHhfaSBcZW5ke2VxdWF0aW9ufQoKQW5kOgoKXGJlZ2lue2VxdWF0aW9ufSBLX3tpLGp9ID0gXGFscGhhXjIgXGV4cCBcQmlnKCAtIFxmcmFjezF9ezIgbF4yfSAoeF9pIC0geF9qKV4yClxCaWcpIFxlbmR7ZXF1YXRpb259CgpUb2dldGhlciB0aGVzZSBjb21wb3NlIHRoZSBmaW5hbApcYmVnaW57YWxpZ259IGYgJiBcc2ltIE11bHRpTm9ybWFsKFxtdV9zICsgXG11X2csIEtfYyArIEtfZykgXFwgeV9pICYgXHNpbSBOb3JtYWwoZl9pLApcc2lnbWEpIFxlbmR7YWxpZ259CgoKCiMjIFN0YW4gKGhpZXJhcmNoaWNhbCkKCmBgYHtyLCBlY2hvPUZBTFNFfQp3cml0ZUxpbmVzKHJlYWRMaW5lcygiL3UvMTYvYXJhc2FsbzEvdW5peC9rb3VsdS9iYXllcy9wcm9qZWN0L2RlYnRfZ3BfaGllcmFyY2hpY2FsLnN0YW4iKSkKYGBgCgoKIyMgTW9kZWwKCmBgYHtyLCBtZXNzYWdlPUZBTFNFLCB3YXJuaW5nPUZBTFNFLCBpbmNsdWRlPVRSVUUsIHJlc3VsdHM9RkFMU0V9CgpkYXRhID0gcmVhZC5jc3YoJ2RhdGEuY3N2Jywgcm93Lm5hbWVzPTEpCmRhdGEgPSBkYXRhLmZyYW1lKHQoZGF0YSkpCnlfaXMgPSB1bmxpc3QoZGF0YSx1c2UubmFtZXMgPSBGKQpjb3VudHJ5X2luZF9pcyA9IHJlcCgxOm5jb2woZGF0YSksZWFjaD0xNSkKeWVhcl9pbmRfaXMgPSByZXAoMToxNSxuY29sKGRhdGEpKQoKTl9pcyA9IGxlbmd0aCh5X2lzKQpOX2NvdW50cmllcyA9IG5jb2woZGF0YSkKTl95ZWFyc19wcmVkID0gMwpOX3llYXJzX2lzID0gMTUKCiMgU2ltcGxlIEdQICsgbGluZWFyIG1vZGVsCmdwX21vZGVsX2hpZXIgPSAnZGVidF9ncF9oaWVyYXJjaGljYWwuc3RhbicKCmZpdF9ncF9oIDwtIHN0YW4oZmlsZT1ncF9tb2RlbF9oaWVyLAogICAgICAgICAgICAgICBkYXRhPWxpc3QoeV9pcyA9IHlfaXMsCiAgICAgICAgICAgICAgICAgICAgICAgICBjb3VudHJ5X2luZF9pcyA9IGNvdW50cnlfaW5kX2lzLAogICAgICAgICAgICAgICAgICAgICAgICAgeWVhcl9pbmRfaXMgPSB5ZWFyX2luZF9pcywKICAgICAgICAgICAgICAgICAgICAgICAgIE5feWVhcnNfaXMgPSAxNSwKICAgICAgICAgICAgICAgICAgICAgICAgIE5feWVhcnNfcHJlZCA9IE5feWVhcnNfcHJlZCwKICAgICAgICAgICAgICAgICAgICAgICAgIE5faXMgPSBOX2lzLAogICAgICAgICAgICAgICAgICAgICAgICAgTl9jb3VudHJpZXMgPSBOX2NvdW50cmllcyksCiAgICAgICAgICAgICAgIGl0ZXI9MTAwMCwKICAgICAgICAgICAgICAgY2hhaW5zPTEwKQoKc2FtcGxlIDwtIGV4dHJhY3QoZml0X2dwX2gpCmBgYApgYGB7cn0Kc2FtcGxlX3ByZWRzIDwtIGRhdGEuZnJhbWUoc2FtcGxlJHlfcHJlZFssc2FtcGxlJGNvdW50cnlfaW5kaWNhdG9yc1sxLF0gPT0gbWF0Y2goIkZJIixjb2xuYW1lcyhkYXRhKSldKQogIApzYW1wbGVfcXVhbnRpbGVzIDwtIGFwcGx5KHNhbXBsZV9wcmVkcywyLGZ1bmN0aW9uKHgpe3F1YW50aWxlKHgsYygwLjA1LDAuNSwwLjk1KSl9KQp4ID0gMToxOApzYW1wbGVfZGF0YSA8LSBkYXRhLmZyYW1lKGxvdyA9IHNhbXBsZV9xdWFudGlsZXNbMSxdLAogICAgICAgICAgICAgICAgICAgICAgICAgIG1lZCA9IHNhbXBsZV9xdWFudGlsZXNbMixdLAogICAgICAgICAgICAgICAgICAgICAgICAgIGhpZ2ggPSBzYW1wbGVfcXVhbnRpbGVzWzMsXSwKICAgICAgICAgICAgICAgICAgICAgICAgICB4ID0geCkKCgp0cmFpbmluZ19kYXRhIDwtIGRhdGEuZnJhbWUoeD0xOjE1LCBvYnM9ZGF0YVssIkZJIl0pCiN2YWxpZGF0aW9uX2RhdGEgPC0gZGF0YS5mcmFtZSh4PTE1OjE4LCBvYnM9eV9vcykKCmNvbHMgPC0gYygiR1AgbWVkaWFuIj0iYmxhY2siLCJHUCA5MCUgcGVyY2VudGlsZSI9ImxpZ2h0Ymx1ZSIsIlRydWUgZGVidCI9InJlZCIpCmdncGxvdChzYW1wbGVfZGF0YSwgYWVzKHg9eCwgeT1tZWQpKSArCiAgZ2VvbV9yaWJib24oYWVzKHg9eCx5bWluPWxvdyx5bWF4PWhpZ2gsY29sb3VyPSJHUCA5MCUgcGVyY2VudGlsZSIpLGFscGhhPTAuNixmaWxsPSJsaWdodGJsdWUiKSArCiAgZ2VvbV9saW5lKGFlcyh4PXgseT1tZWQsY29sb3VyPSJHUCBtZWRpYW4iKSkgKwogIGdlb21fcG9pbnQoKSArCiAgZ2VvbV9saW5lKGRhdGE9dHJhaW5pbmdfZGF0YSwgYWVzKHg9eCwgeT1vYnMsY29sb3VyPSJUcnVlIGRlYnQiKSkgKwogICNnZW9tX2xpbmUoZGF0YT12YWxpZGF0aW9uX2RhdGEsIGFlcyh4PXgsIHk9b2JzLGNvbG91cj0iVHJ1ZSBkZWJ0IiksIGxpbmV0eXBlPSJkYXNoZWQiKSArCiAgeGxhYigneWVhcicpICsgCiAgeWxhYignRGVidCAoJSBvZiBHRFApJykgKyAKICBnZ3RpdGxlKHNwcmludGYoJ0dvdmVybm1lbnQgZGVidCwgY291bnRyeSA9IEZpbmxhbmQgKGhpZXJhcmNoaWNhbCBtb2RlbCknKSkgKwogIHNjYWxlX2NvbG91cl9tYW51YWwobmFtZT0iQ29sb3VyIix2YWx1ZXM9Y29scykgKwogIHRoZW1lKHBsb3QudGl0bGUgPSBlbGVtZW50X3RleHQoaGp1c3QgPSAwLjUpKQoKYGBgCgpCeSBhZGRpbmcgbW9yZSBkYXRhIHVzaW5nIGhpZXJhcmNpY2FsIGdhdXNzaWFuIHByb2Nlc3Nlcywgb3VyIG1vZGVsIGdpdmVzIHRpZ2h0ZXIgZXJyb3IgYm91bmRzIHRoYW4gdGhlIHNlcGFyYXRlIG1vZGVsIGZvciBGaW5sYW5kLiBBcyBleHBlY3RlZCB0aGUgdW5jZXJ0YWludHkgaW5jcmVhc2VzIGFmdGVyIHRoZXJlIGFyZSBubyBvYnNlcnZlZCBkYXRhIHBvaW50cy4KCiMjIENvbnZlcmdlbmNlIGRpYWdub3N0aWNzCgojIyMgJFxoYXR7Un0kIHZhbHVlcwoKYGBge3J9CmZpdF9ncF9oCmBgYAoKJFxoYXR7Un0kIHZhbHVlcyBhcmUgYWxsIDEgc28gTUNNQyBzaW11bGF0aW9ucyBjb252ZXJnZW5jZXMuCgojIyMgVHJlZSBkZXB0aAoKYGBge3J9CmNoZWNrX3RyZWVkZXB0aChmaXRfZ3BfaCkKYGBgCgpPdXIgbW9kZWwgaGl0IHRoZSBtYXhpbXVtIHRyYWplY3RvcnkgbGVuZ3RoIGluIE1DTUMgc2ltdWxhdGlvbnMgbXVsdGlwbGUgdGltZXMgc28gd2UgbmVlZCB0byBpbmNyZWFzZSBvdXIgdHJlZWRlcHRoIHRvIGZpeCB0aGlzIGlzc3VlLgoKIyMjIEUtQkZNSQoKYGBge3J9CmNoZWNrX2VuZXJneShmaXRfZ3BfaCkKYGBgCgoKIyMjIERpdmVyZ2VuY2UKCmBgYHtyfQpjaGVja19kaXYoZml0X2dwX2gpCmBgYAoKYGBge3J9CmNfZGFyayA8LSBjKCIjOEYyNzI3ODAiKQpncmVlbiA8LSBjKCIjMDBGRjAwODAiKQoKcGFydGl0aW9uIDwtIHBhcnRpdGlvbl9kaXYoZml0X2dwX2gpCmRpdl9wYXJhbXMgPC0gcGFydGl0aW9uW1sxXV0Kbm9uZGl2X3BhcmFtcyA8LSBwYXJ0aXRpb25bWzJdXQoKcGFyKG1hciA9IGMoNCwgNCwgMC41LCAwLjUpKQpwbG90KG5vbmRpdl9wYXJhbXMkZ2xvYmFsX2FscGhhLCBub25kaXZfcGFyYW1zJHNpZ21hLAogICAgIGNvbD1jX2RhcmssIHBjaD0xNiwgY2V4PTAuOCwgeGxhYj0ibGVuZ3RoX3BhcmFtIiwgeWxhYj0ic2lnbWEiKQpwb2ludHMoZGl2X3BhcmFtcyRnbG9iYWxfYWxwaGEsIGRpdl9wYXJhbXMkc2lnbWEsCiAgICAgICBjb2w9Z3JlZW4sIHBjaD0xNiwgY2V4PTAuOCkKYGBgCldlIHRyaWVkIGRpZmZlcmVudCBwYXJhbWV0ZXIgY29tcGFyaXNvbnMgYnV0IHdlIHdlcmUgbm90IGFibGUgdG8gZmluZCBhbnkgY2xlYXIgcHJvYmxlbSBiYXNlZCBvbiB0aGUgcGxvdHMsIHdoeSBvdXIgZGl2ZXJnZW5jZXMgZG9lcyBub3QgdmFuaXNoLiBXZSB3aWxsIHRyeSB0byBmaXggdGhlIGlzc3VlIGJ5IHR1bmluZyB0aGUgTUNNQyBwYXJhbWV0ZXJzIGFkYXB0IGRlbHRhIGFuZCB0cmVlIGRlcHRoLgoKCgojIyBGaXhpbmcgdGhlIGhpZXJhcmNoaWNhbCBtb2RlbAoKYGBge3IsIG1lc3NhZ2U9RkFMU0UsIHdhcm5pbmc9RkFMU0UsIGluY2x1ZGU9VFJVRSwgcmVzdWx0cz1GQUxTRX0KZml0X2dwX2ggPC0gc3RhbihmaWxlPWdwX21vZGVsX2hpZXIsCiAgICAgICAgICAgICAgIGRhdGE9bGlzdCh5X2lzID0geV9pcywKICAgICAgICAgICAgICAgICAgICAgICAgIGNvdW50cnlfaW5kX2lzID0gY291bnRyeV9pbmRfaXMsCiAgICAgICAgICAgICAgICAgICAgICAgICB5ZWFyX2luZF9pcyA9IHllYXJfaW5kX2lzLAogICAgICAgICAgICAgICAgICAgICAgICAgTl95ZWFyc19pcyA9IDE1LAogICAgICAgICAgICAgICAgICAgICAgICAgTl95ZWFyc19wcmVkID0gTl95ZWFyc19wcmVkLAogICAgICAgICAgICAgICAgICAgICAgICAgTl9pcyA9IE5faXMsCiAgICAgICAgICAgICAgICAgICAgICAgICBOX2NvdW50cmllcyA9IE5fY291bnRyaWVzKSwKICAgICAgICAgICAgICAgaXRlcj0xMDAwLAogICAgICAgICAgICAgICBjaGFpbnM9MTAsY29udHJvbD1saXN0KGFkYXB0X2RlbHRhPTAuOTksbWF4X3RyZWVkZXB0aD0xNSkpCmBgYAoKIyMjIENvbnZlcmdlbmNlCgpgYGB7cn0KcHJpbnQoZml0X2dwX2gpCmNoZWNrX2RpdihmaXRfZ3BfaCkKY2hlY2tfdHJlZWRlcHRoKGZpdF9ncF9oKQpjaGVja19lbmVyZ3koZml0X2dwX2gpCmBgYAoKQnkgaW5jcmVhc2luZyB0aGUgYWRhcHRfZGVsdGEgYW5kIHRyZWUgZGVwdGggdGhlIG1vZGVsIGhhdmUgbm8gcHJvYmxlbXMgd2l0aCBkaXZlcmdlbmNlcyBhbmQgZG9lcyBub3QgaGl0IHRoZSBtYXhpbXVtIHRyYWplY3RvcnkgdGhyZXNob2xkLgoKIyMgUG9zdGVyaW9yIHByZWRpY3RpdmUgY2hlY2tpbmcKCmBgYHtyfQpsaWJyYXJ5KGJheWVzcGxvdCkKc2FtcGxlID0gZXh0cmFjdChmaXRfZ3BfaCkKcHBjX2RlbnNfb3ZlcmxheSh5PWRhdGFbLCJGSSJdLCB5cmVwPXNhbXBsZSR5X3ByZWRbLHNhbXBsZSRjb3VudHJ5X2luZGljYXRvcnNbMSxdID09IG1hdGNoKCJGSSIsY29sbmFtZXMoZGF0YSkpXVssMToxNV0pK2dndGl0bGUoIlBvc3RlcmlvciBwcmVkaWN0aXZlIGZvciBGaW5sYW5kIChoaWVyYXJjaGljYWwgbW9kZWwpIikgKyB4bGFiKCJEZWJ0ICglIG9mIEdEUCkiKQpgYGAKCkhpZXJhcmNoaWNhbCBtb2RlbCBmb2xsb3dzIHRoZSBlbXByaWNhbCBkYXRhIGRpc3RyaWJ1dGlvbiBxdWl0ZSB3ZWxsLCBiZXR0ZXIgdGhhbiB0aGUgc2VwYXJhdGUgbW9kZWwgZm9yIEZpbmxhbmQuIFZhcmlhbmNlIGF0IDQwLTQ1ICgyMDA4IGZpbmFuY2lhbCBjcmlzaXMpIGlzIG5vdCBhcyBsYXJnZSBhcyBpdCBpcyB3aXRoIHRoZSBzZXBhcmF0ZSBtb2RlbC4gTW9yZW92ZXIgdmFyaWFuY2UgYmV0d2VlbiBkaWZmZXJlbnQgZHJhd3MgaW4gZ2VuZXJhbCBpcyBsb3dlciB0aGFuIHdpdGggc2VwYXJhdGUgbW9kZWwuIEFzIGV4cGVjdGVkIGJ5IGluY3JlYXNpbmcgdGhlIGNvbXBsZXhpdHkgb2Ygb3VyIG1vZGVsIGJ5IGFkZGluZyBoaWVyYXJjaHkgbWFrZXMgcG9zdGVyaW9yIGZpdCB0aGUgb2JzZXJ2ZWQgcGhlbm9tZW5hIGJldHRlci4gCgoKIyMgQ29tcGFyaW5nIHNlcGFyYXRlIGFuZCBoaWVyYXJjaGljYWwgbW9kZWwKCk1vZGVsIGNvbXBhcmlzaW9uIGNhbiBiZSBhY2hpZXZlZCB1c2luZyBwc2lzLUxPTy4KCmBgYHtyfQpsaWJyYXJ5KGxvbykKbGxfc2VwIDwtIGV4dHJhY3RfbG9nX2xpayhzYW1wbGUxJG0pCmxfc2VwIDwtIGxvbyhsbF9zZXApCmxsX2hpZXIgPC0gZXh0cmFjdF9sb2dfbGlrKGZpdF9ncF9oKVssY291bnRyeV9pbmRfaXMgPT0gbWF0Y2goIkZJIixjb2xuYW1lcyhkYXRhKSldCmxfaGllciA8LSBsb28obGxfaGllcikKcHJpbnQobF9zZXApCnByaW50KGxfaGllcikKYGBgCgpQYXJldG8tayB2YWx1ZXMgYXJlIHF1aXRlIHJlYXNvbmFibGUuIE9ubHkgYSBmZXcgdmFsdWVzIGFyZSBhYm92ZSAwLjcgZm9yIGJvdGggbW9kZWxzLiBCYXNlZCBvbiBMT08gdmFsdWVzIGhpZXJhcmNoaWNhbCBtb2RlbCBwZXJmb3JtcyBzbGlnaHRseSBiZXR0ZXIgYmVjYXVzZSBlbHBkX2xvbyB2YWx1ZSBpcyBoaWdoZXIuIEhvd2V2ZXIgd2UgY291bGQgc3RpbGwgY29uc2lkZXIgdXNpbmcgbW9yZSBjb21wdXRhdGlvbmFsbHkgaGVhdnkgbW9kZWwgY29tcGFyaXNvbiBtZXRob2RzLCBiZWNhdXNlIHRoZXJlIGFyZSBzdGlsbCBhIGZldyB2YWx1ZXMgb3ZlciAwLjcsIHdoaWNoIG1ha2VzIHVzZSBvZiBwc2lzLUxPTyBmb3IgbW9kZWwgY29tcGFyaXNvbiBwb3NzaWJseSBxdWVzdGlvbmFibGUuIENyb3NzIHZhbGlkYXRpb24gY291bGQgYmUgdXNlZCBpbnN0ZWFkIG9mIHBzaXMtTE9PLgoKIyMgUHJvYmxlbXMgYW5kIEZ1dHVyZSBpbXByb3ZlbWVudHMKCkZpbmRpbmcgYSBwcm9wZXIgcHJpb3Igd2FzIGEgYml0IHRyaWNreSwgYmVjYXVzZSB3ZSBoYWQgb25seSBhIGZldyBkYXRhcG9pbnRzLCBhdCBsZWFzdCBmb3IgdGhlIHNlcGFyYXRlIG1vZGVsLiBJZiB0aGUgcHJpb3IgaXMgdG9vIHZhZ3VlLCBnYXVzc2lhbiBwcm9jZXNzIG92ZXJmaXRzIGJhZGx5IGFuZCBpZiBkZWNheSBzcGVlZCBvZiB0aGUgbGVuZ3RoIHNjYWxlIGlzIHRvbyBzbG93LCBnYXVzc2lhbiBwcm9jZXNzIHNtb290aGVucyB0aGUgY3VydmUgdG9vIG11Y2ggYW5kIG1ha2VzIHRoZSBwcmVkaWN0aW9uIHRvbyBsaW5lYXIuICAKCk9uZSBvcHRpb24gdG8gaW5jcmVhc2UgdGhlIHVzZWZ1bG5lc3Mgb2YgdGhlIGhpZXJhcmNoaWNhbCBtb2RlbCBpcyB0byBhZGQgbW9yZSBsYXllcnMgdG8gdGhlIGhpZXJhcmNoeS4gRS5nIHdlIGNvdWxkIGdyb3VwIGNvdW50cmllcyBpbnRvIHN1Ymdyb3VwcyBzbyB0aGF0IG5vcmRpYyBjb3VudHJpZXMgaGF2ZSB0aGVpciBvd24gImdyb3VwZWQiIGdhdXNzaWFuIHByb2Nlc3Mgb24gdG9wIG9mIHRoZSBnbG9iYWwgYW5kIGNvdW50cnkgc3BlY2ZpYyBncCBtb2RlbC4gVGhlcmUgaXMgZGVmaW5hdGVseSBtb3JlIHNpbWlsYXJpdGllcyBiZXR3ZWVuIG5vcmRpYyBjb3VudHJpZXMgdGhhbiB3aXRoIGUuZyBGaW5sYW5kIGFuZCBJdGFseS4gIAoKTW9kZWwgY29tcGFyaXNvbiBzaG91bGQgYmUgaW1wcm92ZWQgYnkgY29tcGFyaW5nIHRoZSBtb2RlbHMgdXNpbmcgY3Jvc3MgdmFsaWRhdGlvbi4gU2Vuc2l0aXZpdHkgYW5hbHlzaXMgc2hvdWxkIGFsc28gYmUgZG9uZSBmb3IgaGllcmFyY2hpY2FsIG1vZGVsLgoKIyMgUmVmZXJlbmNlcwoKKiBodHRwOi8vbWMtc3Rhbi5vcmcvZXZlbnRzL3N0YW5jb24yMDE3LW5vdGVib29rcy9zdGFuY29uMjAxNy1zYW5kZXJzLWxlaS1tYXNzLXNob290aW5ncy5odG1sICg0LjEyLjIwMTcpCiogaHR0cHM6Ly9naXRodWIuY29tL2F2ZWh0YXJpL0JEQV9SX2RlbW9zL2Jsb2IvbWFzdGVyL2RlbW9zX3JzdGFuL3BvaXNzb25fZ3Auc3RhbiAoNC4xMi4yMDE3KQoqIGh0dHA6Ly9tYy1zdGFuLm9yZy9ldmVudHMvc3RhbmNvbjIwMTctbm90ZWJvb2tzL3N0YW5jb24yMDE3LXRyYW5ndWNjaS1oaWVyYXJjaGljYWwtZ3BzLnBkZiAoNC4xMi4yMDE3KQoqIGh0dHBzOi8vcmVwb3NpdG9yeS50dWRlbGZ0Lm5sL2lzbGFuZG9yYS9vYmplY3QvdXVpZDo3ZjM5MTZjOS03OTVjLTQwM2YtYmVjOC01ZWI2Y2FhM2UzOTgvZGF0YXN0cmVhbS9PQkovZG93bmxvYWQgKDQuMTIuMjAxNykKKiBodHRwOi8vbWMtc3Rhbi5vcmcvdXNlcnMvZG9jdW1lbnRhdGlvbi9jYXNlLXN0dWRpZXMvcnN0YW5fd29ya2Zsb3cuaHRtbCAoOC4xMi4yMDE3KQoKIyMgQXBwZW5kaXg6IExpY2Vuc2VzCiogQ29kZSAmY29weTsgMjAxNywgVHVvbWFzIFlsaXBpaGEsIE9zc2kgQXJhc2FsbywgbGljZW5zZWQgdW5kZXIgQlNELTMuCiogVGV4dCAmY29weTsgMjAxNywgVHVvbWFzIFlsaXBpaGEsIE9zc2kgQXJhc2FsbywgbGljZW5zZWQgdW5kZXIgQ0MtQlktTkMgNC4wLgoK

Bayesian data analysis project

Data and the problem

Model description (separate model for each country)

Prior choices

μ0,μbμ0,μb\mu_0 , \mu_b

αα\alpha

length_param

Stan (separate model)

R code (separate model)

Convergence diagnostics

R̂ R^\hat{R} diagnostics

Tree Depth

E-BFMI

Divergence

Fixing the model

Visualization

Posterior predictive checking

Sensitivity analysis

Length_param

Model comparison for different priors (LOO)

Informative prior

vague prior

non-informative prior

Second model (hierarchical)

Stan (hierarchical)

Model

Convergence diagnostics

R̂ R^\hat{R} values

Tree depth

E-BFMI

Divergence

Fixing the hierarchical model

Convergence

Posterior predictive checking

Comparing separate and hierarchical model

Problems and Future improvements

References

Appendix: Licenses

$\mu_0 , \mu_b$

$\alpha$

$\hat{R}$ diagnostics

$\hat{R}$ values