Assignment 3
library(lattice)
load("Feldman.rda")
Feldman <- within(Feldman, logret <- log10(retention))
xyplot(logret ~ time | group, data=Feldman, groups=subjID, type="b")
#Task 1
# Load necessary libraries
library(rstan)
# Set a seed for reproducibility
set.seed(12345)
# Load the data
load("Feldman.rda")
Feldman <- within(Feldman, logret <- log10(retention))
# Prepare data for Stan
stan_data <- list(
N = nrow(Feldman),
J = length(unique(Feldman$subjID)),
logret = Feldman$logret,
time = Feldman$time,
group = as.numeric(Feldman$group == "Lung"), # 0 for Liver, 1 for Lung
subjID = as.numeric(as.factor(Feldman$subjID))
)
# Stan model code as a string
stan_model_code <- "
data {
int<lower=0> N; // Number of observations
int<lower=0> J; // Number of subjects
vector[N] logret; // Log retention
vector[N] time; // Time
int<lower=0, upper=1> group[N]; // Treatment group (0 for Liver, 1 for Lung)
int<lower=1, upper=J> subjID[N]; // Subject ID
}
parameters {
real beta0; // Intercept for population
real beta1; // Slope for population
real beta2; // Group effect for population
real beta3; // Interaction of time and group for population
vector[J] b0; // Random intercepts for subjects
vector[J] b1; // Random slopes for subjects
real<lower=0> sigma; // Standard deviation of residuals
real<lower=0> sigma_b0; // Standard deviation of random intercepts
real<lower=0> sigma_b1; // Standard deviation of random slopes
}
model {
// Priors
beta0 ~ normal(0, 1000);
beta1 ~ normal(0, 1000);
beta2 ~ normal(0, 1000);
beta3 ~ normal(0, 1000);
b0 ~ normal(0, sigma_b0);
b1 ~ normal(0, sigma_b1);
sigma ~ cauchy(0, 2.5); // A weakly informative prior for residual SD
sigma_b0 ~ cauchy(0, 2.5); // Prior for SD of random intercepts
sigma_b1 ~ cauchy(0, 2.5); // Prior for SD of random slopes
// Likelihood
for (n in 1:N) {
logret[n] ~ normal(beta0 + b0[subjID[n]] + (beta1 + b1[subjID[n]]) * time[n] +
beta2 * group[n] + beta3 * time[n] * group[n], sigma);
}
}
"
# Compile and fit the model
stan_model <- stan(model_code = stan_model_code,
data = stan_data,
chains = 4,
iter = 3000,
warmup = 2000,
thin = 1,
control = list(adapt_delta = 0.999, max_treedepth = 30))
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 2.2e-05 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.22 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 3000 [ 0%] (Warmup)
Chain 1: Iteration: 300 / 3000 [ 10%] (Warmup)
Chain 1: Iteration: 600 / 3000 [ 20%] (Warmup)
Chain 1: Iteration: 900 / 3000 [ 30%] (Warmup)
Chain 1: Iteration: 1200 / 3000 [ 40%] (Warmup)
Chain 1: Iteration: 1500 / 3000 [ 50%] (Warmup)
Chain 1: Iteration: 1800 / 3000 [ 60%] (Warmup)
Chain 1: Iteration: 2001 / 3000 [ 66%] (Sampling)
Chain 1: Iteration: 2300 / 3000 [ 76%] (Sampling)
Chain 1: Iteration: 2600 / 3000 [ 86%] (Sampling)
Chain 1: Iteration: 2900 / 3000 [ 96%] (Sampling)
Chain 1: Iteration: 3000 / 3000 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 101.326 seconds (Warm-up)
Chain 1: 6.489 seconds (Sampling)
Chain 1: 107.815 seconds (Total)
Chain 1:
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
Chain 2:
Chain 2: Gradient evaluation took 2.2e-05 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.22 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2:
Chain 2:
Chain 2: Iteration: 1 / 3000 [ 0%] (Warmup)
Chain 2: Iteration: 300 / 3000 [ 10%] (Warmup)
Chain 2: Iteration: 600 / 3000 [ 20%] (Warmup)
Chain 2: Iteration: 900 / 3000 [ 30%] (Warmup)
Chain 2: Iteration: 1200 / 3000 [ 40%] (Warmup)
Chain 2: Iteration: 1500 / 3000 [ 50%] (Warmup)
Chain 2: Iteration: 1800 / 3000 [ 60%] (Warmup)
Chain 2: Iteration: 2001 / 3000 [ 66%] (Sampling)
Chain 2: Iteration: 2300 / 3000 [ 76%] (Sampling)
Chain 2: Iteration: 2600 / 3000 [ 86%] (Sampling)
Chain 2: Iteration: 2900 / 3000 [ 96%] (Sampling)
Chain 2: Iteration: 3000 / 3000 [100%] (Sampling)
Chain 2:
Chain 2: Elapsed Time: 17.575 seconds (Warm-up)
Chain 2: 4.132 seconds (Sampling)
Chain 2: 21.707 seconds (Total)
Chain 2:
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
Chain 3:
Chain 3: Gradient evaluation took 2.2e-05 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.22 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3:
Chain 3:
Chain 3: Iteration: 1 / 3000 [ 0%] (Warmup)
Chain 3: Iteration: 300 / 3000 [ 10%] (Warmup)
Chain 3: Iteration: 600 / 3000 [ 20%] (Warmup)
Chain 3: Iteration: 900 / 3000 [ 30%] (Warmup)
Chain 3: Iteration: 1200 / 3000 [ 40%] (Warmup)
Chain 3: Iteration: 1500 / 3000 [ 50%] (Warmup)
Chain 3: Iteration: 1800 / 3000 [ 60%] (Warmup)
Chain 3: Iteration: 2001 / 3000 [ 66%] (Sampling)
Chain 3: Iteration: 2300 / 3000 [ 76%] (Sampling)
Chain 3: Iteration: 2600 / 3000 [ 86%] (Sampling)
Chain 3: Iteration: 2900 / 3000 [ 96%] (Sampling)
Chain 3: Iteration: 3000 / 3000 [100%] (Sampling)
Chain 3:
Chain 3: Elapsed Time: 127.212 seconds (Warm-up)
Chain 3: 4.532 seconds (Sampling)
Chain 3: 131.744 seconds (Total)
Chain 3:
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
Chain 4:
Chain 4: Gradient evaluation took 2.2e-05 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.22 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4:
Chain 4:
Chain 4: Iteration: 1 / 3000 [ 0%] (Warmup)
Chain 4: Iteration: 300 / 3000 [ 10%] (Warmup)
Chain 4: Iteration: 600 / 3000 [ 20%] (Warmup)
Chain 4: Iteration: 900 / 3000 [ 30%] (Warmup)
Chain 4: Iteration: 1200 / 3000 [ 40%] (Warmup)
Chain 4: Iteration: 1500 / 3000 [ 50%] (Warmup)
Chain 4: Iteration: 1800 / 3000 [ 60%] (Warmup)
Chain 4: Iteration: 2001 / 3000 [ 66%] (Sampling)
Chain 4: Iteration: 2300 / 3000 [ 76%] (Sampling)
Chain 4: Iteration: 2600 / 3000 [ 86%] (Sampling)
Chain 4: Iteration: 2900 / 3000 [ 96%] (Sampling)
Chain 4: Iteration: 3000 / 3000 [100%] (Sampling)
Chain 4:
Chain 4: Elapsed Time: 22.397 seconds (Warm-up)
Chain 4: 3.888 seconds (Sampling)
Chain 4: 26.285 seconds (Total)
Chain 4: Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#bulk-essWarning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#tail-ess# Check model summary
print(summary(stan_model))$summary
mean se_mean sd 2.5% 25% 50%
beta0 2.024906e+00 6.337154e-04 0.022079392 1.980905e+00 2.010996e+00 2.025355e+00
beta1 -1.723271e-02 3.089716e-05 0.001398659 -1.994999e-02 -1.808812e-02 -1.724828e-02
beta2 1.251194e-03 9.712981e-04 0.031958039 -6.043907e-02 -1.816694e-02 1.425501e-03
beta3 6.532802e-03 4.832099e-05 0.002097835 2.257251e-03 5.272547e-03 6.541376e-03
b0[1] 5.583051e-03 5.316153e-04 0.020972776 -3.350320e-02 -5.197414e-03 3.304195e-03
b0[2] -1.912831e-03 6.044918e-04 0.021861535 -5.083955e-02 -1.246212e-02 -8.405477e-04
b0[3] -1.373625e-02 7.025396e-04 0.022639727 -6.838903e-02 -2.533536e-02 -9.379403e-03
b0[4] 9.677163e-03 5.812269e-04 0.022050468 -2.914484e-02 -2.496085e-03 6.059195e-03
b0[5] 7.712592e-03 5.623023e-04 0.021763489 -3.252922e-02 -3.860017e-03 4.742194e-03
b0[6] 4.308754e-03 5.028697e-04 0.021008960 -3.557722e-02 -6.432229e-03 2.565692e-03
b0[7] -1.877328e-02 7.667821e-04 0.024058029 -7.422441e-02 -3.080465e-02 -1.437509e-02
b0[8] 5.880920e-03 5.288951e-04 0.021677914 -3.440263e-02 -5.334630e-03 3.417584e-03
b1[1] 3.691544e-04 2.583407e-05 0.001312005 -2.166346e-03 -2.652701e-04 2.185243e-04
b1[2] 1.182482e-03 4.659465e-05 0.001542196 -9.862103e-04 1.008611e-04 8.431851e-04
b1[3] -5.913758e-04 2.880915e-05 0.001368114 -3.699924e-03 -1.250166e-03 -3.467007e-04
b1[4] -9.728751e-04 4.423699e-05 0.001498672 -4.698568e-03 -1.674310e-03 -6.388372e-04
b1[5] 4.808358e-04 3.130077e-05 0.001382374 -2.012959e-03 -2.046589e-04 2.824450e-04
b1[6] 2.167115e-04 2.712617e-05 0.001386970 -2.450328e-03 -4.390001e-04 9.432923e-05
b1[7] -1.004742e-04 2.946662e-05 0.001433728 -3.042147e-03 -7.587833e-04 -5.563282e-05
b1[8] -5.695279e-04 3.671897e-05 0.001533915 -4.219460e-03 -1.239789e-03 -2.790559e-04
sigma 5.938971e-02 1.490168e-04 0.006242188 4.858149e-02 5.512610e-02 5.889083e-02
sigma_b0 2.454331e-02 8.360013e-04 0.018210897 2.514949e-03 1.167690e-02 2.091970e-02
sigma_b1 1.567044e-03 5.855260e-05 0.001224813 8.085694e-05 7.046569e-04 1.330172e-03
lp__ 2.085011e+02 5.932339e-01 8.458501564 1.942033e+02 2.028592e+02 2.074780e+02
75% 97.5% n_eff Rhat
beta0 2.038634e+00 2.067328132 1213.9077 1.003245
beta1 -1.638660e-02 -0.014488189 2049.2123 1.000740
beta2 2.027020e-02 0.065844662 1082.5679 1.002633
beta3 7.804653e-03 0.010726910 1884.8245 1.000488
b0[1] 1.549822e-02 0.053355074 1556.3848 1.000261
b0[2] 9.463502e-03 0.040355086 1307.9177 1.001264
b0[3] 1.447784e-04 0.023249026 1038.4862 1.001845
b0[4] 2.021127e-02 0.063133223 1439.2767 1.002154
b0[5] 1.845263e-02 0.056496334 1498.0191 1.001445
b0[6] 1.442557e-02 0.050627699 1745.4127 1.001520
b0[7] -2.484796e-03 0.016691460 984.4104 1.001275
b0[8] 1.571402e-02 0.054261457 1679.9485 1.000631
b1[1] 1.002806e-03 0.003234237 2579.2008 1.001633
b1[2] 1.990308e-03 0.004874952 1095.4861 1.002980
b1[3] 1.373052e-04 0.001776466 2255.1921 1.000009
b1[4] -1.618865e-05 0.001256866 1147.7368 1.000921
b1[5] 1.095136e-03 0.003628077 1950.4771 1.001401
b1[6] 8.320392e-04 0.003237832 2614.3082 1.000456
b1[7] 5.168463e-04 0.002858084 2367.4054 1.000173
b1[8] 2.170931e-04 0.002034871 1745.1054 1.000129
sigma 6.311638e-02 0.073394081 1754.7021 1.001737
sigma_b0 3.241239e-02 0.069840159 474.5138 1.004469
sigma_b1 2.123070e-03 0.004685584 437.5700 1.003900
lp__ 2.134299e+02 227.437053962 203.2988 1.033377
$c_summary
, , chains = chain:1
stats
parameter mean sd 2.5% 25% 50% 75%
beta0 2.023916e+00 0.021247565 1.983844e+00 2.009602e+00 2.024043e+00 2.037828e+00
beta1 -1.718250e-02 0.001354812 -1.965397e-02 -1.804671e-02 -1.722374e-02 -1.638346e-02
beta2 2.353281e-03 0.028907082 -5.301151e-02 -1.489642e-02 2.549706e-03 2.056286e-02
beta3 6.480785e-03 0.001908195 2.653433e-03 5.368270e-03 6.483942e-03 7.695777e-03
b0[1] 5.588002e-03 0.020527559 -3.295839e-02 -5.113704e-03 3.209495e-03 1.550946e-02
b0[2] -1.727410e-03 0.021350019 -4.952946e-02 -1.136000e-02 -4.499417e-04 8.511883e-03
b0[3] -1.324625e-02 0.020984782 -6.592020e-02 -2.364539e-02 -8.938302e-03 -1.486180e-04
b0[4] 9.356560e-03 0.020899886 -2.530543e-02 -2.116736e-03 5.564948e-03 1.899126e-02
b0[5] 7.188038e-03 0.018842073 -2.710010e-02 -4.199147e-03 4.796242e-03 1.799045e-02
b0[6] 4.808964e-03 0.019527044 -3.347707e-02 -6.055155e-03 2.767500e-03 1.413783e-02
b0[7] -1.890344e-02 0.022623182 -6.725884e-02 -3.173224e-02 -1.470639e-02 -2.371555e-03
b0[8] 5.819701e-03 0.019553898 -2.932757e-02 -5.232959e-03 3.220313e-03 1.606952e-02
b1[1] 3.116791e-04 0.001233549 -2.122091e-03 -2.421363e-04 1.594822e-04 8.642716e-04
b1[2] 1.093143e-03 0.001438169 -8.620420e-04 4.133658e-05 7.429927e-04 1.972581e-03
b1[3] -5.984394e-04 0.001280431 -3.820333e-03 -1.203917e-03 -3.022774e-04 5.140340e-05
b1[4] -9.431084e-04 0.001491292 -4.778850e-03 -1.533975e-03 -5.899292e-04 -1.618865e-05
b1[5] 4.812888e-04 0.001199983 -1.631466e-03 -1.630940e-04 2.293392e-04 1.064452e-03
b1[6] 2.296150e-04 0.001193260 -2.070198e-03 -3.566442e-04 8.186107e-05 8.004733e-04
b1[7] -8.800324e-05 0.001274217 -2.794436e-03 -7.345593e-04 -4.043711e-05 4.475721e-04
b1[8] -4.918157e-04 0.001313623 -3.694407e-03 -1.084487e-03 -2.252215e-04 1.633577e-04
sigma 5.916645e-02 0.006129827 4.839651e-02 5.495943e-02 5.873390e-02 6.291879e-02
sigma_b0 2.330054e-02 0.016413703 1.809263e-03 1.176908e-02 1.996661e-02 3.083520e-02
sigma_b1 1.441528e-03 0.001156024 3.684136e-05 6.231084e-04 1.230251e-03 1.986189e-03
lp__ 2.101870e+02 9.878318739 1.955633e+02 2.038041e+02 2.083974e+02 2.146543e+02
stats
parameter 97.5%
beta0 2.065394121
beta1 -0.014578770
beta2 0.059790949
beta3 0.010073019
b0[1] 0.050630404
b0[2] 0.039392372
b0[3] 0.021446481
b0[4] 0.061248332
b0[5] 0.050113310
b0[6] 0.049327377
b0[7] 0.015736787
b0[8] 0.053039402
b1[1] 0.002991638
b1[2] 0.004250843
b1[3] 0.001595127
b1[4] 0.001139649
b1[5] 0.003321146
b1[6] 0.002947484
b1[7] 0.002680806
b1[8] 0.001711925
sigma 0.072862408
sigma_b0 0.062072301
sigma_b1 0.004436967
lp__ 236.356567610
, , chains = chain:2
stats
parameter mean sd 2.5% 25% 50% 75%
beta0 2.023642e+00 0.022035102 1.978253e+00 2.009244e+00 2.024367e+00 2.037517e+00
beta1 -1.720830e-02 0.001440684 -1.999802e-02 -1.809479e-02 -1.722909e-02 -1.635126e-02
beta2 2.738858e-03 0.032050807 -5.728715e-02 -1.787003e-02 3.475386e-03 2.277222e-02
beta3 6.451979e-03 0.002266848 1.408364e-03 5.089163e-03 6.507000e-03 7.839468e-03
b0[1] 5.696772e-03 0.020058908 -3.242438e-02 -4.665371e-03 3.276181e-03 1.532797e-02
b0[2] -1.099502e-03 0.021563893 -4.787876e-02 -1.157754e-02 -3.386867e-04 8.670033e-03
b0[3] -1.174345e-02 0.021687054 -6.374139e-02 -2.286217e-02 -7.859751e-03 9.735313e-04
b0[4] 1.002209e-02 0.022002303 -2.876378e-02 -1.806601e-03 6.464122e-03 2.049923e-02
b0[5] 7.040417e-03 0.021746983 -3.550914e-02 -4.172605e-03 4.017410e-03 1.696587e-02
b0[6] 3.136577e-03 0.021640011 -4.134884e-02 -6.823196e-03 1.454056e-03 1.381596e-02
b0[7] -1.744304e-02 0.023135424 -7.956432e-02 -2.799869e-02 -1.280421e-02 -1.622769e-03
b0[8] 6.052369e-03 0.021177628 -4.017276e-02 -4.261356e-03 3.378168e-03 1.564942e-02
b1[1] 4.273497e-04 0.001354123 -2.081547e-03 -2.765061e-04 3.212891e-04 1.164601e-03
b1[2] 1.207542e-03 0.001540121 -1.048807e-03 1.907741e-04 9.321144e-04 1.965952e-03
b1[3] -6.438184e-04 0.001448726 -3.599196e-03 -1.338500e-03 -4.857310e-04 1.766144e-04
b1[4] -1.002860e-03 0.001548088 -4.650577e-03 -1.774587e-03 -6.695506e-04 1.028872e-05
b1[5] 5.844036e-04 0.001579732 -2.123456e-03 -2.357017e-04 3.895480e-04 1.209961e-03
b1[6] 3.019068e-04 0.001505727 -2.316264e-03 -4.737739e-04 1.740262e-04 9.623464e-04
b1[7] -7.866469e-05 0.001492955 -2.853270e-03 -7.825747e-04 -7.574953e-05 5.644437e-04
b1[8] -5.430645e-04 0.001633542 -4.079783e-03 -1.310589e-03 -2.793910e-04 3.245269e-04
sigma 5.985662e-02 0.006416164 4.863158e-02 5.547812e-02 5.950033e-02 6.356547e-02
sigma_b0 2.380400e-02 0.018136586 1.882610e-03 1.097545e-02 1.995245e-02 3.199268e-02
sigma_b1 1.663089e-03 0.001264874 1.926994e-04 8.020052e-04 1.423763e-03 2.150960e-03
lp__ 2.077813e+02 7.173286807 1.947809e+02 2.029030e+02 2.072376e+02 2.123526e+02
stats
parameter 97.5%
beta0 2.064501787
beta1 -0.014356694
beta2 0.067658190
beta3 0.010789768
b0[1] 0.051671241
b0[2] 0.045471450
b0[3] 0.024763626
b0[4] 0.066189472
b0[5] 0.057111003
b0[6] 0.052535004
b0[7] 0.015817094
b0[8] 0.054263239
b1[1] 0.003210215
b1[2] 0.004863944
b1[3] 0.001829479
b1[4] 0.001231971
b1[5] 0.004155082
b1[6] 0.003515459
b1[7] 0.003072110
b1[8] 0.002317933
sigma 0.073950631
sigma_b0 0.071530716
sigma_b1 0.005042476
lp__ 223.826820642
, , chains = chain:3
stats
parameter mean sd 2.5% 25% 50% 75%
beta0 2.026352e+00 0.022406889 1.986787e+00 2.012541e+00 2.026382e+00 2.039467e+00
beta1 -1.729103e-02 0.001449595 -2.029907e-02 -1.815433e-02 -1.728846e-02 -1.634789e-02
beta2 -8.216498e-04 0.031466302 -6.006608e-02 -2.122646e-02 -7.127958e-04 1.873003e-02
beta3 6.640422e-03 0.002141272 2.732001e-03 5.278657e-03 6.533483e-03 7.878016e-03
b0[1] 5.300429e-03 0.021405536 -3.613565e-02 -6.108611e-03 3.706201e-03 1.557299e-02
b0[2] -3.396971e-03 0.022930678 -5.620785e-02 -1.384619e-02 -1.786845e-03 9.075252e-03
b0[3] -1.442444e-02 0.024087101 -7.356057e-02 -2.686358e-02 -9.447579e-03 3.589006e-04
b0[4] 8.768842e-03 0.021940383 -3.122067e-02 -2.854070e-03 6.356379e-03 2.027044e-02
b0[5] 8.588895e-03 0.020626320 -2.803840e-02 -2.931316e-03 5.777373e-03 1.808779e-02
b0[6] 5.140046e-03 0.019699470 -3.346091e-02 -5.754893e-03 3.599148e-03 1.501845e-02
b0[7] -1.874014e-02 0.023356584 -7.423815e-02 -3.062429e-02 -1.479713e-02 -3.138302e-03
b0[8] 6.368326e-03 0.021397932 -3.148328e-02 -5.530351e-03 4.303063e-03 1.571753e-02
b1[1] 4.220177e-04 0.001363886 -2.115617e-03 -3.273809e-04 2.803946e-04 1.132689e-03
b1[2] 1.323217e-03 0.001667245 -1.282056e-03 1.906619e-04 9.848508e-04 2.178291e-03
b1[3] -5.960343e-04 0.001375400 -3.582093e-03 -1.309774e-03 -3.539044e-04 1.716139e-04
b1[4] -9.652679e-04 0.001478738 -4.542183e-03 -1.746624e-03 -6.799078e-04 -1.479093e-05
b1[5] 4.282829e-04 0.001430948 -2.287915e-03 -2.170070e-04 2.966465e-04 1.083629e-03
b1[6] 1.785921e-04 0.001433667 -2.734479e-03 -4.845530e-04 8.749760e-05 8.140320e-04
b1[7] -1.693822e-04 0.001534881 -3.668340e-03 -8.253109e-04 -7.874767e-05 5.267289e-04
b1[8] -6.610783e-04 0.001568990 -4.540613e-03 -1.365382e-03 -3.588847e-04 1.854008e-04
sigma 5.914241e-02 0.006166209 4.837433e-02 5.516786e-02 5.874422e-02 6.273955e-02
sigma_b0 2.483546e-02 0.017764308 3.823109e-03 1.227866e-02 2.097083e-02 3.298380e-02
sigma_b1 1.639155e-03 0.001190491 1.672512e-04 7.945924e-04 1.414133e-03 2.197855e-03
lp__ 2.073082e+02 7.033421535 1.937804e+02 2.026406e+02 2.070232e+02 2.119727e+02
stats
parameter 97.5%
beta0 2.070421838
beta1 -0.014501114
beta2 0.063509702
beta3 0.011348731
b0[1] 0.053377272
b0[2] 0.039051211
b0[3] 0.023476760
b0[4] 0.059181066
b0[5] 0.056679166
b0[6] 0.049718199
b0[7] 0.017912351
b0[8] 0.053609960
b1[1] 0.003349958
b1[2] 0.005283464
b1[3] 0.001865914
b1[4] 0.001473093
b1[5] 0.003533512
b1[6] 0.003167517
b1[7] 0.003079449
b1[8] 0.001994734
sigma 0.072571275
sigma_b0 0.069483167
sigma_b1 0.004531129
lp__ 221.452040432
, , chains = chain:4
stats
parameter mean sd 2.5% 25% 50% 75%
beta0 2.025714e+00 0.022518067 1.975559e+00 2.012684e+00 2.026706e+00 2.039890e+00
beta1 -1.724900e-02 0.001345982 -1.983012e-02 -1.803631e-02 -1.726858e-02 -1.646632e-02
beta2 7.342865e-04 0.035033280 -6.521687e-02 -1.961574e-02 9.827236e-04 1.916991e-02
beta3 6.558019e-03 0.002056825 2.246618e-03 5.433370e-03 6.644980e-03 7.815127e-03
b0[1] 5.747001e-03 0.021879054 -3.278757e-02 -4.872654e-03 2.783084e-03 1.549895e-02
b0[2] -1.427441e-03 0.021525893 -4.997014e-02 -1.297031e-02 -1.203622e-03 1.077091e-02
b0[3] -1.553088e-02 0.023515459 -6.982805e-02 -2.796335e-02 -1.094430e-02 -3.921792e-04
b0[4] 1.056116e-02 0.023285954 -2.885982e-02 -3.181368e-03 5.977038e-03 2.125486e-02
b0[5] 8.033019e-03 0.025318392 -3.680421e-02 -4.047231e-03 4.565575e-03 2.009274e-02
b0[6] 4.149427e-03 0.022952636 -3.785304e-02 -7.420909e-03 2.884046e-03 1.455623e-02
b0[7] -2.000650e-02 0.026851352 -8.009618e-02 -3.367665e-02 -1.478011e-02 -2.614785e-03
b0[8] 5.283286e-03 0.024326780 -4.036978e-02 -5.914450e-03 2.960953e-03 1.507443e-02
b1[1] 3.155713e-04 0.001289477 -2.242831e-03 -2.318913e-04 1.431903e-04 8.872423e-04
b1[2] 1.106026e-03 0.001505426 -8.405998e-04 3.631716e-05 6.850727e-04 1.864477e-03
b1[3] -5.272109e-04 0.001362198 -3.632627e-03 -1.130751e-03 -2.677898e-04 1.689195e-04
b1[4] -9.802644e-04 0.001477067 -4.856078e-03 -1.655389e-03 -6.163007e-04 -2.555188e-05
b1[5] 4.293678e-04 0.001284636 -1.789278e-03 -2.141014e-04 2.367283e-04 9.885218e-04
b1[6] 1.567319e-04 0.001393393 -2.615906e-03 -4.070028e-04 5.286837e-05 7.551754e-04
b1[7] -6.584648e-05 0.001418973 -2.883982e-03 -6.860055e-04 -4.672827e-05 4.932632e-04
b1[8] -5.821529e-04 0.001596107 -4.443428e-03 -1.170695e-03 -2.358218e-04 1.871289e-04
sigma 5.939336e-02 0.006235641 4.910995e-02 5.508099e-02 5.859979e-02 6.309398e-02
sigma_b0 2.623323e-02 0.020213600 2.806421e-03 1.257572e-02 2.278252e-02 3.436122e-02
sigma_b1 1.524404e-03 0.001272627 7.863349e-05 5.811218e-04 1.220148e-03 2.140383e-03
lp__ 2.087279e+02 9.113391302 1.931011e+02 2.019725e+02 2.076438e+02 2.148080e+02
stats
parameter 97.5%
beta0 2.066212005
beta1 -0.014531955
beta2 0.067454840
beta3 0.010560671
b0[1] 0.057052398
b0[2] 0.039677801
b0[3] 0.023551980
b0[4] 0.067584780
b0[5] 0.060960732
b0[6] 0.055331403
b0[7] 0.017894874
b0[8] 0.056583085
b1[1] 0.003120592
b1[2] 0.004858241
b1[3] 0.001777651
b1[4] 0.001025595
b1[5] 0.003531295
b1[6] 0.003210256
b1[7] 0.002809532
b1[8] 0.002250319
sigma 0.073731301
sigma_b0 0.075002175
sigma_b1 0.004906267
lp__ 228.750721140Bayesian estimates (c)
# Extracting estimates and Monte Carlo Standard Errors (MCSE)
beta0_estimate <- fit_summary$summary["beta0", "mean"]
beta1_estimate <- fit_summary$summary["beta1", "mean"]
beta2_estimate <- fit_summary$summary["beta2", "mean"]
beta3_estimate <- fit_summary$summary["beta3", "mean"]
sigma_estimate <- fit_summary$summary["sigma", "mean"]
sigma_b0_estimate <- fit_summary$summary["sigma_b0", "mean"]
sigma_b1_estimate <- fit_summary$summary["sigma_b1", "mean"]
beta0_mcse <- fit_summary$summary["beta0", "se_mean"]
beta1_mcse <- fit_summary$summary["beta1", "se_mean"]
beta2_mcse <- fit_summary$summary["beta2", "se_mean"]
beta3_mcse <- fit_summary$summary["beta3", "se_mean"]
sigma_mcse <- fit_summary$summary["sigma", "se_mean"]
sigma_b0_mcse <- fit_summary$summary["sigma_b0", "se_mean"]
sigma_b1_mcse <- fit_summary$summary["sigma_b1", "se_mean"]
# Create a data frame to nicely display the results
estimates_df <- data.frame(
Parameter = c("beta0", "beta1", "beta2", "beta3", "sigma", "sigma_b0", "sigma_b1"),
Estimate = c(beta0_estimate, beta1_estimate, beta2_estimate, beta3_estimate, sigma_estimate, sigma_b0_estimate, sigma_b1_estimate),
MCSE = c(beta0_mcse, beta1_mcse, beta2_mcse, beta3_mcse, sigma_mcse, sigma_b0_mcse, sigma_b1_mcse)
)
# Print the data frame
print(estimates_df)NA#Task 2
# Load the required library
library(rstan)
# Set a random seed for reproducibility
set.seed(12345)
# Load the task2.rda data
load("task2.rda")
# Split the data into x (Design A) and y (Design B)
x <- z[1:30] # Measurements from Design A
y <- z[31:60] # Measurements from Design B
# Prepare the data list for STAN
data_list1 <- list(N = length(x), x = x, y = y)Model 1
# Define the STAN model for Model 1
model1_code <- "
data {
int N; // Number of observations in each group
vector[N] x; // Measurements from Design A
vector[N] y; // Measurements from Design B
}
parameters {
real mu_x; // Mean of Design A
real mu_y; // Mean of Design B
real<lower=0> sigma; // Common standard deviation
}
model {
// Priors
mu_x ~ normal(0, 10);
mu_y ~ normal(0, 10);
sigma ~ cauchy(0, 5);
// Likelihood
x ~ normal(mu_x, sigma);
y ~ normal(mu_y, sigma);
}
generated quantities {
real p_mu_y_ge_mu_x = mu_y >= mu_x ? 1 : 0; // P(mu_y >= mu_x)
real future_x = normal_rng(mu_x, sigma); // Future observation for Design A
real future_y = normal_rng(mu_y, sigma); // Future observation for Design B
real p_future_x_ge_15 = future_x >= 15 ? 1 : 0; // P(future_x >= 15)
real p_future_y_ge_15 = future_y >= 15 ? 1 : 0; // P(future_y >= 15)
real p_y_ge_x = future_y >= future_x ? 1 : 0; // P(future_y >= future_x)
}
"
# Compile and fit the model to the data
fit1 <- stan(model_code = model1_code, data = data_list1, iter = 2000, chains = 4)
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 3e-06 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.03 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 0.025 seconds (Warm-up)
Chain 1: 0.027 seconds (Sampling)
Chain 1: 0.052 seconds (Total)
Chain 1:
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
Chain 2:
Chain 2: Gradient evaluation took 6e-06 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.06 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2:
Chain 2:
Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 2:
Chain 2: Elapsed Time: 0.031 seconds (Warm-up)
Chain 2: 0.032 seconds (Sampling)
Chain 2: 0.063 seconds (Total)
Chain 2:
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
Chain 3:
Chain 3: Gradient evaluation took 3e-06 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.03 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3:
Chain 3:
Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 3:
Chain 3: Elapsed Time: 0.034 seconds (Warm-up)
Chain 3: 0.035 seconds (Sampling)
Chain 3: 0.069 seconds (Total)
Chain 3:
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
Chain 4:
Chain 4: Gradient evaluation took 3e-06 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.03 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4:
Chain 4:
Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 4:
Chain 4: Elapsed Time: 0.031 seconds (Warm-up)
Chain 4: 0.035 seconds (Sampling)
Chain 4: 0.066 seconds (Total)
Chain 4: # Print summary of the fit
print(summary(fit1))$summary
mean se_mean sd 2.5% 25% 50% 75%
mu_x 12.109596 0.007226210 0.4341484 11.275988 11.81054 12.116057 12.402912
mu_y 13.122300 0.007501895 0.4214616 12.295176 12.84016 13.121550 13.398948
sigma 2.323866 0.003935379 0.2211372 1.947507 2.16872 2.306605 2.459517
p_mu_y_ge_mu_x 0.949250 0.004404499 0.2195142 0.000000 1.00000 1.000000 1.000000
future_x 12.099714 0.038435775 2.3905116 7.441807 10.50977 12.097630 13.686327
future_y 13.120722 0.038470911 2.4100558 8.427594 11.47471 13.119818 14.779810
p_future_x_ge_15 0.107250 0.005017760 0.3094695 0.000000 0.00000 0.000000 0.000000
p_future_y_ge_15 0.220000 0.006709718 0.4142981 0.000000 0.00000 0.000000 0.000000
p_y_ge_x 0.617250 0.007833031 0.4861188 0.000000 0.00000 1.000000 1.000000
lp__ -80.825750 0.030427439 1.2615314 -84.049778 -81.41198 -80.469701 -79.907254
97.5% n_eff Rhat
mu_x 12.946857 3609.568 1.0006042
mu_y 13.949222 3156.269 0.9997819
sigma 2.799889 3157.551 1.0003712
p_mu_y_ge_mu_x 1.000000 2483.890 0.9998766
future_x 16.772769 3868.214 0.9999450
future_y 17.794491 3924.545 1.0002545
p_future_x_ge_15 1.000000 3803.785 0.9999482
p_future_y_ge_15 1.000000 3812.566 0.9995521
p_y_ge_x 1.000000 3851.458 1.0004687
lp__ -79.397422 1718.958 1.0001240
$c_summary
, , chains = chain:1
stats
parameter mean sd 2.5% 25% 50% 75% 97.5%
mu_x 12.105460 0.4284703 11.305980 11.814768 12.09956 12.396429 12.939721
mu_y 13.106572 0.4352231 12.270909 12.805558 13.09779 13.393407 13.971422
sigma 2.331369 0.2243623 1.950381 2.163347 2.31984 2.468295 2.798182
p_mu_y_ge_mu_x 0.942000 0.2338604 0.000000 1.000000 1.00000 1.000000 1.000000
future_x 12.115814 2.3767335 7.523198 10.565887 12.14594 13.677269 16.767625
future_y 13.068853 2.4319688 8.291640 11.362686 13.05781 14.739492 17.427126
p_future_x_ge_15 0.108000 0.3105357 0.000000 0.000000 0.00000 0.000000 1.000000
p_future_y_ge_15 0.216000 0.4117202 0.000000 0.000000 0.00000 0.000000 1.000000
p_y_ge_x 0.610000 0.4879940 0.000000 0.000000 1.00000 1.000000 1.000000
lp__ -80.857155 1.2597415 -84.208241 -81.450499 -80.49905 -79.972237 -79.388361
, , chains = chain:2
stats
parameter mean sd 2.5% 25% 50% 75% 97.5%
mu_x 12.09830 0.4257862 11.306493 11.793537 12.105884 12.372922 12.930948
mu_y 13.12138 0.4376483 12.268780 12.840474 13.105980 13.406914 14.020688
sigma 2.33074 0.2178657 1.946207 2.186047 2.311213 2.471041 2.820486
p_mu_y_ge_mu_x 0.95500 0.2074079 0.000000 1.000000 1.000000 1.000000 1.000000
future_x 12.02259 2.3751603 7.427263 10.449931 11.947877 13.611228 16.722901
future_y 13.20744 2.4725897 8.343695 11.591669 13.234651 14.879522 17.850643
p_future_x_ge_15 0.10300 0.3041110 0.000000 0.000000 0.000000 0.000000 1.000000
p_future_y_ge_15 0.23200 0.4223202 0.000000 0.000000 0.000000 0.000000 1.000000
p_y_ge_x 0.63100 0.4827754 0.000000 0.000000 1.000000 1.000000 1.000000
lp__ -80.83098 1.2813859 -84.023186 -81.400445 -80.466938 -79.910523 -79.382098
, , chains = chain:3
stats
parameter mean sd 2.5% 25% 50% 75% 97.5%
mu_x 12.126511 0.4403127 11.265616 11.825045 12.12593 12.448997 12.930555
mu_y 13.118803 0.4082888 12.294570 12.847254 13.12909 13.374915 13.919780
sigma 2.315601 0.2210360 1.946185 2.154761 2.30379 2.453965 2.785715
p_mu_y_ge_mu_x 0.950000 0.2180540 0.000000 1.000000 1.00000 1.000000 1.000000
future_x 12.098253 2.3275773 7.659918 10.527106 12.07948 13.565568 16.832157
future_y 13.052108 2.3613615 8.417911 11.567927 13.07064 14.690550 17.560664
p_future_x_ge_15 0.103000 0.3041110 0.000000 0.000000 0.00000 0.000000 1.000000
p_future_y_ge_15 0.208000 0.4060799 0.000000 0.000000 0.00000 0.000000 1.000000
p_y_ge_x 0.623000 0.4848774 0.000000 0.000000 1.00000 1.000000 1.000000
lp__ -80.810476 1.2359694 -83.953410 -81.326171 -80.46817 -79.911731 -79.421998
, , chains = chain:4
stats
parameter mean sd 2.5% 25% 50% 75% 97.5%
mu_x 12.108115 0.4419456 11.261580 11.80051 12.125710 12.381342 12.981745
mu_y 13.142441 0.4033995 12.353586 12.86329 13.167516 13.413659 13.905458
sigma 2.317753 0.2210954 1.945341 2.16872 2.296916 2.442714 2.812251
p_mu_y_ge_mu_x 0.950000 0.2180540 0.000000 1.00000 1.000000 1.000000 1.000000
future_x 12.162200 2.4814020 6.997038 10.53334 12.168942 13.854852 16.976909
future_y 13.154488 2.3728979 8.750903 11.42732 13.105874 14.868186 17.926128
p_future_x_ge_15 0.115000 0.3191816 0.000000 0.00000 0.000000 0.000000 1.000000
p_future_y_ge_15 0.224000 0.4171307 0.000000 0.00000 0.000000 0.000000 1.000000
p_y_ge_x 0.605000 0.4890953 0.000000 0.00000 1.000000 1.000000 1.000000
lp__ -80.804388 1.2698035 -84.048758 -81.43142 -80.421561 -79.860513 -79.400814#Monte Carlo error estimates model 1
# Assuming 'fit1' is the fitted model object from Model 1
# Extract the samples for relevant parameters
samples <- extract(fit1)
# Calculate Monte Carlo errors for key parameters
mc_error_mu_x <- sd(samples$mu_x) / sqrt(length(samples$mu_x))
mc_error_mu_y <- sd(samples$mu_y) / sqrt(length(samples$mu_y))
mc_error_sigma <- sd(samples$sigma) / sqrt(length(samples$sigma))
# Calculating Monte Carlo errors for the generated quantities
mc_error_p_mu_y_ge_mu_x <- sd(samples$p_mu_y_ge_mu_x) / sqrt(length(samples$p_mu_y_ge_mu_x))
mc_error_p_future_x_ge_15 <- sd(samples$p_future_x_ge_15) / sqrt(length(samples$p_future_x_ge_15))
mc_error_p_future_y_ge_15 <- sd(samples$p_future_y_ge_15) / sqrt(length(samples$p_future_y_ge_15))
mc_error_p_y_ge_x <- sd(samples$p_y_ge_x) / sqrt(length(samples$p_y_ge_x))
# Print the Monte Carlo errors
cat("Monte Carlo Error for mu_x: ", mc_error_mu_x, "\n")Monte Carlo Error for mu_x: 0.006864489 cat("Monte Carlo Error for mu_y: ", mc_error_mu_y, "\n")Monte Carlo Error for mu_y: 0.006663892 cat("Monte Carlo Error for sigma: ", mc_error_sigma, "\n")Monte Carlo Error for sigma: 0.003496486 cat("Monte Carlo Error for P(mu_y >= mu_x): ", mc_error_p_mu_y_ge_mu_x, "\n")Monte Carlo Error for P(mu_y >= mu_x): 0.003470824 cat("Monte Carlo Error for P(future_x >= 15): ", mc_error_p_future_x_ge_15, "\n")Monte Carlo Error for P(future_x >= 15): 0.004893143 cat("Monte Carlo Error for P(future_y >= 15): ", mc_error_p_future_y_ge_15, "\n")Monte Carlo Error for P(future_y >= 15): 0.006550628 cat("Monte Carlo Error for P(future_y >= future_x): ", mc_error_p_y_ge_x, "\n")Monte Carlo Error for P(future_y >= future_x): 0.007686214 #Model 2
# Load the required library
library(rstan)
# Set a random seed for reproducibility
set.seed(12345)
# Load the task2.rda data
load("task2.rda")
# Split the data into x (Design A) and y (Design B)
x <- z[1:30] # Measurements from Design A
y <- z[31:60] # Measurements from Design B
# Prepare the data list for STAN
data_list1 <- list(N = length(x), x = x, y = y)
# Define the STAN model for Model 2
model2_code <- "
data {
int N; // Number of observations in each group
vector[N] x; // Measurements from Design A
vector[N] y; // Measurements from Design B
}
parameters {
real mu_x; // Mean of Design A
real mu_y; // Mean of Design B
real<lower=0> sigma_x; // Standard deviation for Design A
real<lower=0> sigma_y; // Standard deviation for Design B
}
model {
// Priors
mu_x ~ normal(0, 10);
mu_y ~ normal(0, 10);
sigma_x ~ cauchy(0, 5);
sigma_y ~ cauchy(0, 5);
// Likelihood
x ~ normal(mu_x, sigma_x);
y ~ normal(mu_y, sigma_y);
}
generated quantities {
real p_mu_y_ge_mu_x = mu_y >= mu_x ? 1 : 0; // P(mu_y >= mu_x)
real future_x = normal_rng(mu_x, sigma_x); // Future observation for Design A
real future_y = normal_rng(mu_y, sigma_y); // Future observation for Design B
real p_future_x_ge_15 = future_x >= 15 ? 1 : 0; // P(future_x >= 15)
real p_future_y_ge_15 = future_y >= 15 ? 1 : 0; // P(future_y >= 15)
real p_y_ge_x = future_y >= future_x ? 1 : 0; // P(future_y >= future_x)
}
"
# Compile and fit the model to the data
fit2 <- stan(model_code = model2_code, data = data_list1, iter = 2000, chains = 4)
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 3e-06 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.03 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 0.029 seconds (Warm-up)
Chain 1: 0.031 seconds (Sampling)
Chain 1: 0.06 seconds (Total)
Chain 1:
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
Chain 2:
Chain 2: Gradient evaluation took 4e-06 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.04 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2:
Chain 2:
Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 2:
Chain 2: Elapsed Time: 0.034 seconds (Warm-up)
Chain 2: 0.035 seconds (Sampling)
Chain 2: 0.069 seconds (Total)
Chain 2:
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
Chain 3:
Chain 3: Gradient evaluation took 4e-06 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.04 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3:
Chain 3:
Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 3:
Chain 3: Elapsed Time: 0.035 seconds (Warm-up)
Chain 3: 0.035 seconds (Sampling)
Chain 3: 0.07 seconds (Total)
Chain 3:
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
Chain 4:
Chain 4: Gradient evaluation took 4e-06 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.04 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4:
Chain 4:
Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 4:
Chain 4: Elapsed Time: 0.036 seconds (Warm-up)
Chain 4: 0.037 seconds (Sampling)
Chain 4: 0.073 seconds (Total)
Chain 4: # Print summary of the fit
print(summary(fit2))$summary
mean se_mean sd 2.5% 25% 50%
mu_x 12.0873873 0.009564187 0.5742594 10.9307963 11.7283703 12.0851764
mu_y 13.1376971 0.002595867 0.1612678 12.8187453 13.0343823 13.1399305
sigma_x 3.2049201 0.007294938 0.4422150 2.4811972 2.8928799 3.1539394
sigma_y 0.8908606 0.002159571 0.1263690 0.6815854 0.8016653 0.8791285
p_mu_y_ge_mu_x 0.9635000 0.003432453 0.1875541 0.0000000 1.0000000 1.0000000
future_x 12.0344974 0.051275203 3.2479998 5.4249873 9.9468547 12.0715840
future_y 13.1290729 0.014961490 0.9260228 11.3467475 12.4810651 13.1320282
p_future_x_ge_15 0.1652500 0.005984519 0.3714525 0.0000000 0.0000000 0.0000000
p_future_y_ge_15 0.0270000 0.002595256 0.1621036 0.0000000 0.0000000 0.0000000
p_y_ge_x 0.6325000 0.007695664 0.4821845 0.0000000 0.0000000 1.0000000
lp__ -61.2986239 0.036802800 1.5088449 -65.1470333 -62.0202671 -60.9134411
75% 97.5% n_eff Rhat
mu_x 12.4582026 13.227362 3605.123 0.9996601
mu_y 13.2406317 13.454672 3859.492 1.0006405
sigma_x 3.4615044 4.238466 3674.717 0.9995315
sigma_y 0.9616113 1.176712 3424.100 0.9994299
p_mu_y_ge_mu_x 1.0000000 1.000000 2985.685 1.0003939
future_x 14.1546983 18.515708 4012.520 0.9996049
future_y 13.7310933 15.026066 3830.837 1.0015075
p_future_x_ge_15 0.0000000 1.000000 3852.547 0.9992890
p_future_y_ge_15 0.0000000 1.000000 3901.437 1.0003490
p_y_ge_x 1.0000000 1.000000 3925.855 1.0000563
lp__ -60.1927204 -59.484950 1680.844 1.0011705
$c_summary
, , chains = chain:1
stats
parameter mean sd 2.5% 25% 50% 75%
mu_x 12.0618262 0.5673623 10.8835183 11.7277750 12.0766074 12.4252098
mu_y 13.1331776 0.1646660 12.8045262 13.0207336 13.1402309 13.2425538
sigma_x 3.2028210 0.4601299 2.4411314 2.8752060 3.1464202 3.4535795
sigma_y 0.8913321 0.1247024 0.6815789 0.8059835 0.8829016 0.9584666
p_mu_y_ge_mu_x 0.9670000 0.1787259 0.0000000 1.0000000 1.0000000 1.0000000
future_x 12.0812455 3.2437850 5.8750285 10.0191604 12.1004265 14.2445906
future_y 13.0644958 0.9546033 11.2411941 12.3944501 13.0434197 13.6839271
p_future_x_ge_15 0.1720000 0.3775693 0.0000000 0.0000000 0.0000000 0.0000000
p_future_y_ge_15 0.0250000 0.1562031 0.0000000 0.0000000 0.0000000 0.0000000
p_y_ge_x 0.6140000 0.4870742 0.0000000 0.0000000 1.0000000 1.0000000
lp__ -61.3445182 1.5556847 -65.2476260 -62.0123723 -60.9933998 -60.2398509
stats
parameter 97.5%
mu_x 13.165243
mu_y 13.442718
sigma_x 4.287033
sigma_y 1.194365
p_mu_y_ge_mu_x 1.000000
future_x 18.313629
future_y 14.974145
p_future_x_ge_15 1.000000
p_future_y_ge_15 0.025000
p_y_ge_x 1.000000
lp__ -59.476897
, , chains = chain:2
stats
parameter mean sd 2.5% 25% 50% 75%
mu_x 12.0905916 0.5587560 10.9297835 11.7500248 12.0935158 12.4536961
mu_y 13.1381992 0.1577947 12.8361780 13.0405861 13.1383621 13.2338021
sigma_x 3.2034535 0.4251043 2.4999273 2.9042186 3.1701513 3.4464788
sigma_y 0.8879969 0.1223379 0.6751883 0.8025271 0.8777372 0.9603933
p_mu_y_ge_mu_x 0.9610000 0.1936918 0.0000000 1.0000000 1.0000000 1.0000000
future_x 12.0168065 3.2470207 5.2666500 9.8934287 12.1404908 14.0826510
future_y 13.1762614 0.9368065 11.5041366 12.5267401 13.1629944 13.7682123
p_future_x_ge_15 0.1570000 0.3639828 0.0000000 0.0000000 0.0000000 0.0000000
p_future_y_ge_15 0.0380000 0.1912919 0.0000000 0.0000000 0.0000000 0.0000000
p_y_ge_x 0.6440000 0.4790548 0.0000000 0.0000000 1.0000000 1.0000000
lp__ -61.2184265 1.4286025 -64.8630636 -61.9529615 -60.8542615 -60.1340541
stats
parameter 97.5%
mu_x 13.255993
mu_y 13.462679
sigma_x 4.132462
sigma_y 1.144946
p_mu_y_ge_mu_x 1.000000
future_x 18.732160
future_y 15.165517
p_future_x_ge_15 1.000000
p_future_y_ge_15 1.000000
p_y_ge_x 1.000000
lp__ -59.443004
, , chains = chain:3
stats
parameter mean sd 2.5% 25% 50% 75%
mu_x 12.094046 0.6254649 10.9075314 11.6987743 12.0901153 12.5013090
mu_y 13.146378 0.1529068 12.8216157 13.0456315 13.1427020 13.2464259
sigma_x 3.211512 0.4339086 2.4849010 2.8954571 3.1581959 3.4792768
sigma_y 0.892668 0.1245718 0.6950184 0.8026061 0.8775821 0.9653594
p_mu_y_ge_mu_x 0.960000 0.1960572 0.0000000 1.0000000 1.0000000 1.0000000
future_x 12.067505 3.2239922 5.2410948 9.9929280 12.1220685 14.2399180
future_y 13.090796 0.9068534 11.2582005 12.4722174 13.1249024 13.6688216
p_future_x_ge_15 0.164000 0.3704608 0.0000000 0.0000000 0.0000000 0.0000000
p_future_y_ge_15 0.022000 0.1467567 0.0000000 0.0000000 0.0000000 0.0000000
p_y_ge_x 0.629000 0.4833142 0.0000000 0.0000000 1.0000000 1.0000000
lp__ -61.296828 1.5067779 -65.1119570 -62.0645516 -60.9326630 -60.1974732
stats
parameter 97.5%
mu_x 13.306862
mu_y 13.445955
sigma_x 4.186975
sigma_y 1.163835
p_mu_y_ge_mu_x 1.000000
future_x 18.551050
future_y 14.971524
p_future_x_ge_15 1.000000
p_future_y_ge_15 0.000000
p_y_ge_x 1.000000
lp__ -59.528866
, , chains = chain:4
stats
parameter mean sd 2.5% 25% 50% 75%
mu_x 12.1030855 0.5420477 11.0289156 11.7366790 12.0812883 12.4631382
mu_y 13.1330331 0.1691000 12.8063531 13.0166566 13.1371581 13.2409395
sigma_x 3.2018935 0.4494826 2.4965580 2.8894807 3.1405916 3.4726175
sigma_y 0.8914456 0.1337050 0.6758908 0.7973931 0.8805194 0.9646941
p_mu_y_ge_mu_x 0.9660000 0.1813198 0.0000000 1.0000000 1.0000000 1.0000000
future_x 11.9724323 3.2806757 5.6578624 9.8574356 11.9377166 14.1467849
future_y 13.1847382 0.9002330 11.4711785 12.5839162 13.1866199 13.7825411
p_future_x_ge_15 0.1680000 0.3740534 0.0000000 0.0000000 0.0000000 0.0000000
p_future_y_ge_15 0.0230000 0.1499783 0.0000000 0.0000000 0.0000000 0.0000000
p_y_ge_x 0.6430000 0.4793545 0.0000000 0.0000000 1.0000000 1.0000000
lp__ -61.3347230 1.5401280 -65.2594251 -62.0452706 -60.8909779 -60.2034907
stats
parameter 97.5%
mu_x 13.180050
mu_y 13.479074
sigma_x 4.263246
sigma_y 1.218050
p_mu_y_ge_mu_x 1.000000
future_x 18.278019
future_y 14.956147
p_future_x_ge_15 1.000000
p_future_y_ge_15 0.000000
p_y_ge_x 1.000000
lp__ -59.502256#Monte Carlo Error Model 2
# Assuming 'fit2' is the fitted model object from Model 2
# Extract the samples for relevant parameters
samples2 <- extract(fit2)
# Calculate Monte Carlo errors for key parameters
mc_error_mu_x2 <- sd(samples2$mu_x) / sqrt(length(samples2$mu_x))
mc_error_mu_y2 <- sd(samples2$mu_y) / sqrt(length(samples2$mu_y))
mc_error_sigma_x2 <- sd(samples2$sigma_x) / sqrt(length(samples2$sigma_x))
mc_error_sigma_y2 <- sd(samples2$sigma_y) / sqrt(length(samples2$sigma_y))
# Calculating Monte Carlo errors for the generated quantities
mc_error_p_mu_y_ge_mu_x2 <- sd(samples2$p_mu_y_ge_mu_x) / sqrt(length(samples2$p_mu_y_ge_mu_x))
mc_error_p_future_x_ge_152 <- sd(samples2$p_future_x_ge_15) / sqrt(length(samples2$p_future_x_ge_15))
mc_error_p_future_y_ge_152 <- sd(samples2$p_future_y_ge_15) / sqrt(length(samples2$p_future_y_ge_15))
mc_error_p_y_ge_x2 <- sd(samples2$p_y_ge_x) / sqrt(length(samples2$p_y_ge_x))
# Print the Monte Carlo errors
cat("Monte Carlo Error for mu_x (Model 2): ", mc_error_mu_x2, "\n")Monte Carlo Error for mu_x (Model 2): 0.009079838 cat("Monte Carlo Error for mu_y (Model 2): ", mc_error_mu_y2, "\n")Monte Carlo Error for mu_y (Model 2): 0.002549868 cat("Monte Carlo Error for sigma_x (Model 2): ", mc_error_sigma_x2, "\n")Monte Carlo Error for sigma_x (Model 2): 0.006992034 cat("Monte Carlo Error for sigma_y (Model 2): ", mc_error_sigma_y2, "\n")Monte Carlo Error for sigma_y (Model 2): 0.00199807 cat("Monte Carlo Error for P(mu_y >= mu_x) (Model 2): ", mc_error_p_mu_y_ge_mu_x2, "\n")Monte Carlo Error for P(mu_y >= mu_x) (Model 2): 0.002965491 cat("Monte Carlo Error for P(future_x >= 15) (Model 2): ", mc_error_p_future_x_ge_152, "\n")Monte Carlo Error for P(future_x >= 15) (Model 2): 0.005873179 cat("Monte Carlo Error for P(future_y >= 15) (Model 2): ", mc_error_p_future_y_ge_152, "\n")Monte Carlo Error for P(future_y >= 15) (Model 2): 0.002563083 cat("Monte Carlo Error for P(future_y >= future_x) (Model 2): ", mc_error_p_y_ge_x2, "\n")Monte Carlo Error for P(future_y >= future_x) (Model 2): 0.007624006 Task 3)
3a)
set.seed(123)
# Required Libraries
library(rstan)
# Load Data
bb <- read.csv("meta.csv", stringsAsFactors = TRUE)
# Create Survivals Variable
bb <- within(bb, Survivals <- Total - Deaths)
# Summarized Table
summary_table <- xtabs(cbind(Deaths, Survivals) ~ Treatment, data=bb)
print(summary_table)
Treatment Deaths Survivals
control 985 8864
treated 826 9615# Extracting Data for the Model
m1 <- sum(bb$Total[bb$Treatment == "control"]) # Total in control
y1 <- sum(bb$Deaths[bb$Treatment == "control"]) # Deaths in control
m2 <- sum(bb$Total[bb$Treatment == "treated"]) # Total in treated
y2 <- sum(bb$Deaths[bb$Treatment == "treated"]) # Deaths in treated
# Define the Stan Model
stan_model_code <- "
data {
int<lower=0> m1; // Total in control
int<lower=0> y1; // Deaths in control
int<lower=0> m2; // Total in treated
int<lower=0> y2; // Deaths in treated
}
parameters {
real<lower=0, upper=1> p1; // Probability of death in control
real<lower=0, upper=1> p2; // Probability of death in treated
}
model {
// Uniform priors
p1 ~ beta(1, 1);
p2 ~ beta(1, 1);
// Likelihood
y1 ~ binomial(m1, p1);
y2 ~ binomial(m2, p2);
}
generated quantities {
real psi = (p2 / (1 - p2)) / (p1 / (1 - p1)); // Odds ratio
}
"
# Data for Stan Model
data_list <- list(
m1 = m1,
y1 = y1,
m2 = m2,
y2 = y2
)
# Model Fitting
fit <- stan(model_code = stan_model_code, data = data_list, iter = 2000, chains = 4)
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 4e-06 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.04 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 0.049 seconds (Warm-up)
Chain 1: 0.035 seconds (Sampling)
Chain 1: 0.084 seconds (Total)
Chain 1:
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
Chain 2:
Chain 2: Gradient evaluation took 5e-06 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.05 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2:
Chain 2:
Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 2:
Chain 2: Elapsed Time: 0.067 seconds (Warm-up)
Chain 2: 0.034 seconds (Sampling)
Chain 2: 0.101 seconds (Total)
Chain 2:
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
Chain 3:
Chain 3: Gradient evaluation took 5e-06 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.05 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3:
Chain 3:
Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 3:
Chain 3: Elapsed Time: 0.035 seconds (Warm-up)
Chain 3: 0.034 seconds (Sampling)
Chain 3: 0.069 seconds (Total)
Chain 3:
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
Chain 4:
Chain 4: Gradient evaluation took 2e-05 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.2 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4:
Chain 4:
Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 4:
Chain 4: Elapsed Time: 0.034 seconds (Warm-up)
Chain 4: 0.035 seconds (Sampling)
Chain 4: 0.069 seconds (Total)
Chain 4: print(fit, pars = "psi")Inference for Stan model: anon_model.
4 chains, each with iter=2000; warmup=1000; thin=1;
post-warmup draws per chain=1000, total post-warmup draws=4000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
psi 0.77 0 0.04 0.7 0.75 0.77 0.8 0.85 3389 1
Samples were drawn using NUTS(diag_e) at Sun Nov 12 20:57:04 2023.
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).Monte Carlo Error
# Required Library
library(rstan)
# Assuming 'fit' is the fitted model object from the previous Stan model run
# Use the monitor function to get diagnostics, including MCSE
diagnostics <- monitor(fit, warmup = 0)Inference for the input samples (4 chains: each with iter = 2000; warmup = 0):
For each parameter, Bulk_ESS and Tail_ESS are crude measures of
effective sample size for bulk and tail quantities respectively (an ESS > 100
per chain is considered good), and Rhat is the potential scale reduction
factor on rank normalized split chains (at convergence, Rhat <= 1.05).# Print the Monte Carlo Standard Errors
print(diagnostics$mcse)NULLlibrary(rstan)
# Stan model code
stan_model_code <- "
data {
int<lower=0> n; // number of observations
int<lower=0> deaths[n]; // number of deaths for each observation
int<lower=0> total[n]; // total number of participants for each observation
int<lower=0, upper=1> treatment[n]; // treatment indicator
}
parameters {
real beta0; // intercept
real beta1; // coefficient for treatment effect
}
model {
for (i in 1:n) {
deaths[i] ~ binomial_logit(total[i], beta0 + beta1 * treatment[i]);
}
// Priors
beta0 ~ normal(0, 10);
beta1 ~ normal(0, 10);
}
generated quantities {
real odds_ratio = exp(beta1);
}
"
# Data preparation (assuming bb is already loaded)
bb <- within(bb, x <- ifelse(Treatment=="treated", 1, 0))
# Data for Stan model
data_list <- list(
n = nrow(bb),
deaths = bb$Deaths,
total = bb$Total,
treatment = bb$x
)
# Set random seed for reproducibility
set.seed(123)
# Fit the model
fit <- stan(model_code = stan_model_code, data = data_list, iter = 2000, chains = 4)
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 7e-06 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.07 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 0.056 seconds (Warm-up)
Chain 1: 0.054 seconds (Sampling)
Chain 1: 0.11 seconds (Total)
Chain 1:
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
Chain 2:
Chain 2: Gradient evaluation took 6e-06 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.06 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2:
Chain 2:
Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 2:
Chain 2: Elapsed Time: 0.056 seconds (Warm-up)
Chain 2: 0.052 seconds (Sampling)
Chain 2: 0.108 seconds (Total)
Chain 2:
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
Chain 3:
Chain 3: Gradient evaluation took 6e-06 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.06 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3:
Chain 3:
Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 3:
Chain 3: Elapsed Time: 0.054 seconds (Warm-up)
Chain 3: 0.057 seconds (Sampling)
Chain 3: 0.111 seconds (Total)
Chain 3:
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
Chain 4:
Chain 4: Gradient evaluation took 6e-06 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.06 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4:
Chain 4:
Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 4:
Chain 4: Elapsed Time: 0.056 seconds (Warm-up)
Chain 4: 0.053 seconds (Sampling)
Chain 4: 0.109 seconds (Total)
Chain 4: # Print the summary for beta1 and odds_ratio
print(fit, pars = c("beta1", "odds_ratio"))Inference for Stan model: anon_model.
4 chains, each with iter=2000; warmup=1000; thin=1;
post-warmup draws per chain=1000, total post-warmup draws=4000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
beta1 -0.26 0 0.05 -0.36 -0.29 -0.26 -0.23 -0.16 1247 1
odds_ratio 0.77 0 0.04 0.70 0.75 0.77 0.80 0.85 1247 1
Samples were drawn using NUTS(diag_e) at Sun Nov 12 20:58:01 2023.
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).Monte Carlo error
# Required Library
library(rstan)
# Assuming 'fit' is the fitted model object from the previous Stan model run
# Use the monitor function to get diagnostics, including MCSE
diagnostics <- monitor(fit, warmup = 0)Inference for the input samples (4 chains: each with iter = 2000; warmup = 0):
For each parameter, Bulk_ESS and Tail_ESS are crude measures of
effective sample size for bulk and tail quantities respectively (an ESS > 100
per chain is considered good), and Rhat is the potential scale reduction
factor on rank normalized split chains (at convergence, Rhat <= 1.05).# Print the Monte Carlo Standard Errors
print(diagnostics$mcse)NULL# Load RStan library
library(rstan)
# Set random seed for reproducibility
set.seed(12345)
# Load your data
bb <- read.csv("meta.csv", stringsAsFactors = TRUE)
# Define the indicator variable for Treatment
bb <- within(bb, x <- ifelse(Treatment=="treated", 1, 0))
# Prepare data for Stan
data_list <- list(
N = nrow(bb),
J = length(unique(bb$Study)), # Total number of studies
group = as.numeric(factor(bb$Study)), # Convert study identifiers to numeric group index
y = bb$Deaths,
m = bb$Total,
x = bb$x
)
# Stan model
stan_model <- "
data {
int<lower=0> N; // Total number of observations
int<lower=0> J; // Total number of groups
int<lower=1, upper=J> group[N]; // Group index for each observation
int<lower=0> y[N]; // Response variable (number of deaths)
int<lower=0> m[N]; // Total participants for each observation
int<lower=0, upper=1> x[N]; // Treatment indicator
}
parameters {
real beta1; // Global treatment effect
real mu_beta0; // Global intercept
real<lower=0> sigma_beta0; // Standard deviation of group-specific intercepts
vector[J] beta0_raw; // Standardized group-specific intercepts
}
transformed parameters {
vector[J] beta0; // Group-specific intercepts
beta0 = mu_beta0 + sigma_beta0 * beta0_raw; // Non-centered parameterization
}
model {
// Priors
beta1 ~ normal(0, 10);
mu_beta0 ~ normal(0, 10);
sigma_beta0 ~ uniform(0, 10);
beta0_raw ~ normal(0, 1);
// Likelihood
for (i in 1:N) {
y[i] ~ binomial_logit(m[i], beta0[group[i]] + beta1 * x[i]);
}
}
generated quantities {
real odds_ratio;
odds_ratio = exp(beta1); // Odds ratio for treatment effect
}
"
# Fit the model
fit <- stan(model_code = stan_model, data = data_list, iter = 2000, chains = 4)
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 9e-06 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.09 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 0.203 seconds (Warm-up)
Chain 1: 0.147 seconds (Sampling)
Chain 1: 0.35 seconds (Total)
Chain 1:
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
Chain 2:
Chain 2: Gradient evaluation took 9e-06 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.09 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2:
Chain 2:
Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 2:
Chain 2: Elapsed Time: 0.221 seconds (Warm-up)
Chain 2: 0.191 seconds (Sampling)
Chain 2: 0.412 seconds (Total)
Chain 2:
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
Chain 3:
Chain 3: Gradient evaluation took 8e-06 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.08 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3:
Chain 3:
Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 3:
Chain 3: Elapsed Time: 0.246 seconds (Warm-up)
Chain 3: 0.145 seconds (Sampling)
Chain 3: 0.391 seconds (Total)
Chain 3:
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
Chain 4:
Chain 4: Gradient evaluation took 9e-06 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.09 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4:
Chain 4:
Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 4:
Chain 4: Elapsed Time: 0.238 seconds (Warm-up)
Chain 4: 0.154 seconds (Sampling)
Chain 4: 0.392 seconds (Total)
Chain 4: # Extract and compute the Bayesian estimate for the odds ratio
odds_ratio <- extract(fit)$odds_ratio
bayesian_estimate <- mean(odds_ratio)
# Print the Bayesian estimate for the odds ratio
cat("Bayesian estimate for the odds ratio (e^beta1):", bayesian_estimate, "\n")Bayesian estimate for the odds ratio (e^beta1): 0.772011 # Extract the samples for beta1
beta1_samples <- extract(fit)$beta1
# Calculate the odds ratio for each sample
odds_ratio_samples <- exp(beta1_samples)
# Calculate the Monte Carlo error for the odds ratio
mc_error <- sd(odds_ratio_samples) / sqrt(length(odds_ratio_samples))
# Print the Monte Carlo error
cat("Monte Carlo error for the odds ratio estimate:", mc_error, "\n")Monte Carlo error for the odds ratio estimate: 0.0006043269 # Load RStan library
library(rstan)
# Set random seed for reproducibility
set.seed(12345)
# Load your data
bb <- read.csv("meta.csv", stringsAsFactors = TRUE)
# Define the indicator variable for Treatment
bb <- within(bb, x <- ifelse(Treatment == "treated", 1, 0))
# Prepare data for Stan
data_list <- list(
N = nrow(bb),
J = length(unique(bb$Study)),
group = as.numeric(factor(bb$Study)),
y = bb$Deaths,
m = bb$Total,
x = bb$x
)
# Stan model string
stan_model <- "
data {
// Data block defines the data types and their constraints
int<lower=0> N; // Total number of observations
int<lower=0> J; // Total number of groups/studies
int<lower=1, upper=J> group[N]; // Group index for each observation, constrained between 1 and J
int<lower=0> y[N]; // Number of deaths (response variable) for each observation
int<lower=0> m[N]; // Total number of participants for each observation
int<lower=0, upper=1> x[N]; // Treatment indicator (0 or 1) for each observation
}
parameters {
// Parameters block defines the model parameters
real mu_beta0; // Global mean for the intercepts
real<lower=0> sigma_beta0; // Standard deviation for the group-specific intercepts (must be positive)
vector[J] beta0_raw; // Raw, non-centered group-specific intercepts
real mu_beta1; // Global mean for the treatment effects
real<lower=0> sigma_beta1; // Standard deviation for the group-specific treatment effects (must be positive)
vector[J] beta1_raw; // Raw, non-centered group-specific treatment effects
}
transformed parameters {
// Transformed parameters block for intermediate calculations
vector[J] beta0; // Actual group-specific intercepts
vector[J] beta1; // Actual group-specific treatment effects
beta0 = mu_beta0 + sigma_beta0 * beta0_raw; // Centering the group-specific intercepts
beta1 = mu_beta1 + sigma_beta1 * beta1_raw; // Centering the group-specific treatment effects
}
model {
// Model block defines the likelihood and priors
// Priors for global means and standard deviations
mu_beta0 ~ normal(0, 10); // Weakly informative normal prior for global mean intercept
sigma_beta0 ~ uniform(0, 10); // Uniform prior for standard deviation of intercepts
beta0_raw ~ normal(0, 1); // Standard normal prior for raw intercepts
mu_beta1 ~ normal(0, 10); // Weakly informative normal prior for global mean treatment effect
sigma_beta1 ~ uniform(0, 10); // Uniform prior for standard deviation of treatment effects
beta1_raw ~ normal(0, 1); // Standard normal prior for raw treatment effects
// Likelihood of the data given the model
for (i in 1:N) {
y[i] ~ binomial_logit(m[i], beta0[group[i]] + beta1[group[i]] * x[i]); // Binomial logit model for each observation
}
}
generated quantities {
// Generated quantities block for post-processing
real exp_mu_beta1 = exp(mu_beta1); // Exponential of mu_beta1 to get the odds ratio
}
"
# Fit the model using Stan
fit <- stan(model_code = stan_model, data = data_list, iter = 2000, chains = 4)
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 1.5e-05 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.15 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 0.403 seconds (Warm-up)
Chain 1: 0.306 seconds (Sampling)
Chain 1: 0.709 seconds (Total)
Chain 1:
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
Chain 2:
Chain 2: Gradient evaluation took 1e-05 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.1 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2:
Chain 2:
Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 2:
Chain 2: Elapsed Time: 0.415 seconds (Warm-up)
Chain 2: 0.309 seconds (Sampling)
Chain 2: 0.724 seconds (Total)
Chain 2:
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
Chain 3:
Chain 3: Gradient evaluation took 1.2e-05 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.12 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3:
Chain 3:
Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 3:
Chain 3: Elapsed Time: 0.397 seconds (Warm-up)
Chain 3: 0.295 seconds (Sampling)
Chain 3: 0.692 seconds (Total)
Chain 3:
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
Chain 4:
Chain 4: Gradient evaluation took 1.2e-05 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.12 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4:
Chain 4:
Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 4:
Chain 4: Elapsed Time: 0.373 seconds (Warm-up)
Chain 4: 0.304 seconds (Sampling)
Chain 4: 0.677 seconds (Total)
Chain 4: # Extract the posterior samples for exp(mu_beta1)
exp_mu_beta1_samples <- extract(fit)$exp_mu_beta1
# Calculate the Bayesian estimate for exp(mu_beta1)
bayesian_estimate_exp_mu_beta1 <- mean(exp_mu_beta1_samples)
# Print the Bayesian estimate for exp(mu_beta1)
cat("Bayesian estimate for exp(mu_beta1):", bayesian_estimate_exp_mu_beta1, "\n")Bayesian estimate for exp(mu_beta1): 0.7778829 # Calculate the standard deviation of the sampled posterior values
std_dev_exp_mu_beta1 <- sd(exp_mu_beta1_samples)
# Calculate the number of samples
num_samples <- length(exp_mu_beta1_samples)
# Calculate the Monte Carlo error
monte_carlo_error_exp_mu_beta1 <- std_dev_exp_mu_beta1 / sqrt(num_samples)
# Print the Monte Carlo error for exp(mu_beta1)
cat("Monte Carlo error for exp(mu_beta1):", monte_carlo_error_exp_mu_beta1, "\n")Monte Carlo error for exp(mu_beta1): 0.0007773792
library(ggplot2)
# Extract model parameters
posterior_samples <- extract(fit)
mu_beta1_samples <- posterior_samples$mu_beta1
mu_beta0_samples <- posterior_samples$mu_beta0
sigma_beta1_samples <- posterior_samples$sigma_beta1
# Simulate expected number of deaths for a new study with 100 participants in each group
n_participants <- 100
n_simulations <- length(mu_beta1_samples)
# Initialize vectors to store simulated deaths
simulated_deaths_control <- numeric(n_simulations)
simulated_deaths_treated <- numeric(n_simulations)
for (i in 1:n_simulations) {
# Compute probabilities
prob_control <- plogis(mu_beta0_samples[i])
prob_treated <- plogis(mu_beta0_samples[i] + mu_beta1_samples[i])
# Simulate deaths
simulated_deaths_control[i] <- rbinom(1, n_participants, prob_control)
simulated_deaths_treated[i] <- rbinom(1, n_participants, prob_treated)
}
# Compute average simulated deaths
avg_deaths_control <- mean(simulated_deaths_control)
avg_deaths_treated <- mean(simulated_deaths_treated)
cat("Average simulated deaths in control group:", avg_deaths_control, "\n")Average simulated deaths in control group: 9.888 cat("Average simulated deaths in treated group:", avg_deaths_treated, "\n")Average simulated deaths in treated group: 7.9225 # Generate plot of probability mass function for the number of deaths in the treated group
df_plot <- data.frame(deaths = simulated_deaths_treated)
ggplot(df_plot, aes(x = deaths)) +
geom_histogram(binwidth = 1, boundary = 0, fill = "blue", color = "black") +
labs(title = "Probability Mass Function for Number of Deaths in Treated Group",
x = "Number of Deaths",
y = "Frequency") +
theme_minimal()
LS0tCnRpdGxlOiAiQXNzaWdubWVudCAzIgpTdHVkZW50SUQ6ICIyMTU5NTM1OSIKb3V0cHV0OiBodG1sX25vdGVib29rCmVkaXRvcl9vcHRpb25zOiAKICBjaHVua19vdXRwdXRfdHlwZTogaW5saW5lCi0tLQoKCgpgYGB7cn0KbGlicmFyeShsYXR0aWNlKQoKbG9hZCgiRmVsZG1hbi5yZGEiKQpGZWxkbWFuIDwtIHdpdGhpbihGZWxkbWFuLCBsb2dyZXQgPC0gbG9nMTAocmV0ZW50aW9uKSkKeHlwbG90KGxvZ3JldCB+IHRpbWUgfCBncm91cCwgZGF0YT1GZWxkbWFuLCBncm91cHM9c3ViaklELCB0eXBlPSJiIikKYGBgCgoKCgojVGFzayAxCgpgYGB7cn0KIyBMb2FkIG5lY2Vzc2FyeSBsaWJyYXJpZXMKbGlicmFyeShyc3RhbikKCiMgU2V0IGEgc2VlZCBmb3IgcmVwcm9kdWNpYmlsaXR5CnNldC5zZWVkKDEyMzQ1KQoKIyBMb2FkIHRoZSBkYXRhCmxvYWQoIkZlbGRtYW4ucmRhIikKRmVsZG1hbiA8LSB3aXRoaW4oRmVsZG1hbiwgbG9ncmV0IDwtIGxvZzEwKHJldGVudGlvbikpCgojIFByZXBhcmUgZGF0YSBmb3IgU3RhbgpzdGFuX2RhdGEgPC0gbGlzdCgKICBOID0gbnJvdyhGZWxkbWFuKSwKICBKID0gbGVuZ3RoKHVuaXF1ZShGZWxkbWFuJHN1YmpJRCkpLAogIGxvZ3JldCA9IEZlbGRtYW4kbG9ncmV0LAogIHRpbWUgPSBGZWxkbWFuJHRpbWUsCiAgZ3JvdXAgPSBhcy5udW1lcmljKEZlbGRtYW4kZ3JvdXAgPT0gIkx1bmciKSwgIyAwIGZvciBMaXZlciwgMSBmb3IgTHVuZwogIHN1YmpJRCA9IGFzLm51bWVyaWMoYXMuZmFjdG9yKEZlbGRtYW4kc3ViaklEKSkKKQoKIyBTdGFuIG1vZGVsIGNvZGUgYXMgYSBzdHJpbmcKc3Rhbl9tb2RlbF9jb2RlIDwtICIKZGF0YSB7CiAgaW50PGxvd2VyPTA+IE47IC8vIE51bWJlciBvZiBvYnNlcnZhdGlvbnMKICBpbnQ8bG93ZXI9MD4gSjsgLy8gTnVtYmVyIG9mIHN1YmplY3RzCiAgdmVjdG9yW05dIGxvZ3JldDsgLy8gTG9nIHJldGVudGlvbgogIHZlY3RvcltOXSB0aW1lOyAvLyBUaW1lCiAgaW50PGxvd2VyPTAsIHVwcGVyPTE+IGdyb3VwW05dOyAvLyBUcmVhdG1lbnQgZ3JvdXAgKDAgZm9yIExpdmVyLCAxIGZvciBMdW5nKQogIGludDxsb3dlcj0xLCB1cHBlcj1KPiBzdWJqSURbTl07IC8vIFN1YmplY3QgSUQKfQpwYXJhbWV0ZXJzIHsKICByZWFsIGJldGEwOyAvLyBJbnRlcmNlcHQgZm9yIHBvcHVsYXRpb24KICByZWFsIGJldGExOyAvLyBTbG9wZSBmb3IgcG9wdWxhdGlvbgogIHJlYWwgYmV0YTI7IC8vIEdyb3VwIGVmZmVjdCBmb3IgcG9wdWxhdGlvbgogIHJlYWwgYmV0YTM7IC8vIEludGVyYWN0aW9uIG9mIHRpbWUgYW5kIGdyb3VwIGZvciBwb3B1bGF0aW9uCiAgdmVjdG9yW0pdIGIwOyAvLyBSYW5kb20gaW50ZXJjZXB0cyBmb3Igc3ViamVjdHMKICB2ZWN0b3JbSl0gYjE7IC8vIFJhbmRvbSBzbG9wZXMgZm9yIHN1YmplY3RzCiAgcmVhbDxsb3dlcj0wPiBzaWdtYTsgLy8gU3RhbmRhcmQgZGV2aWF0aW9uIG9mIHJlc2lkdWFscwogIHJlYWw8bG93ZXI9MD4gc2lnbWFfYjA7IC8vIFN0YW5kYXJkIGRldmlhdGlvbiBvZiByYW5kb20gaW50ZXJjZXB0cwogIHJlYWw8bG93ZXI9MD4gc2lnbWFfYjE7IC8vIFN0YW5kYXJkIGRldmlhdGlvbiBvZiByYW5kb20gc2xvcGVzCn0KbW9kZWwgewogIC8vIFByaW9ycwogIGJldGEwIH4gbm9ybWFsKDAsIDEwMDApOwogIGJldGExIH4gbm9ybWFsKDAsIDEwMDApOwogIGJldGEyIH4gbm9ybWFsKDAsIDEwMDApOwogIGJldGEzIH4gbm9ybWFsKDAsIDEwMDApOwogIGIwIH4gbm9ybWFsKDAsIHNpZ21hX2IwKTsKICBiMSB+IG5vcm1hbCgwLCBzaWdtYV9iMSk7CiAgc2lnbWEgfiBjYXVjaHkoMCwgMi41KTsgLy8gQSB3ZWFrbHkgaW5mb3JtYXRpdmUgcHJpb3IgZm9yIHJlc2lkdWFsIFNECiAgc2lnbWFfYjAgfiBjYXVjaHkoMCwgMi41KTsgLy8gUHJpb3IgZm9yIFNEIG9mIHJhbmRvbSBpbnRlcmNlcHRzCiAgc2lnbWFfYjEgfiBjYXVjaHkoMCwgMi41KTsgLy8gUHJpb3IgZm9yIFNEIG9mIHJhbmRvbSBzbG9wZXMKCiAgLy8gTGlrZWxpaG9vZAogIGZvciAobiBpbiAxOk4pIHsKICAgIGxvZ3JldFtuXSB+IG5vcm1hbChiZXRhMCArIGIwW3N1YmpJRFtuXV0gKyAoYmV0YTEgKyBiMVtzdWJqSURbbl1dKSAqIHRpbWVbbl0gKyAKICAgICAgICAgICAgICAgICAgICAgICBiZXRhMiAqIGdyb3VwW25dICsgYmV0YTMgKiB0aW1lW25dICogZ3JvdXBbbl0sIHNpZ21hKTsKICB9CgoKfQoiCgoKIyBDb21waWxlIGFuZCBmaXQgdGhlIG1vZGVsCnN0YW5fbW9kZWwgPC0gc3Rhbihtb2RlbF9jb2RlID0gc3Rhbl9tb2RlbF9jb2RlLCAKICAgICAgICAgICAgICAgICAgIGRhdGEgPSBzdGFuX2RhdGEsIAogICAgICAgICAgICAgICAgICAgY2hhaW5zID0gNCwgCiAgICAgICAgICAgICAgICAgICBpdGVyID0gMzAwMCwgCiAgICAgICAgICAgICAgICAgICB3YXJtdXAgPSAyMDAwLCAKICAgICAgICAgICAgICAgICAgIHRoaW4gPSAxLCAKICAgICAgICAgICAgICAgICAgIGNvbnRyb2wgPSBsaXN0KGFkYXB0X2RlbHRhID0gMC45OTksIG1heF90cmVlZGVwdGggPSAzMCkpCgojIENoZWNrIG1vZGVsIHN1bW1hcnkKcHJpbnQoc3VtbWFyeShzdGFuX21vZGVsKSkKYGBgCgpCYXllc2lhbiBlc3RpbWF0ZXMgKGMpCgpgYGB7cn0KIyBFeHRyYWN0aW5nIGVzdGltYXRlcyBhbmQgTW9udGUgQ2FybG8gU3RhbmRhcmQgRXJyb3JzIChNQ1NFKQpiZXRhMF9lc3RpbWF0ZSA8LSBmaXRfc3VtbWFyeSRzdW1tYXJ5WyJiZXRhMCIsICJtZWFuIl0KYmV0YTFfZXN0aW1hdGUgPC0gZml0X3N1bW1hcnkkc3VtbWFyeVsiYmV0YTEiLCAibWVhbiJdCmJldGEyX2VzdGltYXRlIDwtIGZpdF9zdW1tYXJ5JHN1bW1hcnlbImJldGEyIiwgIm1lYW4iXQpiZXRhM19lc3RpbWF0ZSA8LSBmaXRfc3VtbWFyeSRzdW1tYXJ5WyJiZXRhMyIsICJtZWFuIl0Kc2lnbWFfZXN0aW1hdGUgPC0gZml0X3N1bW1hcnkkc3VtbWFyeVsic2lnbWEiLCAibWVhbiJdCnNpZ21hX2IwX2VzdGltYXRlIDwtIGZpdF9zdW1tYXJ5JHN1bW1hcnlbInNpZ21hX2IwIiwgIm1lYW4iXQpzaWdtYV9iMV9lc3RpbWF0ZSA8LSBmaXRfc3VtbWFyeSRzdW1tYXJ5WyJzaWdtYV9iMSIsICJtZWFuIl0KCmJldGEwX21jc2UgPC0gZml0X3N1bW1hcnkkc3VtbWFyeVsiYmV0YTAiLCAic2VfbWVhbiJdCmJldGExX21jc2UgPC0gZml0X3N1bW1hcnkkc3VtbWFyeVsiYmV0YTEiLCAic2VfbWVhbiJdCmJldGEyX21jc2UgPC0gZml0X3N1bW1hcnkkc3VtbWFyeVsiYmV0YTIiLCAic2VfbWVhbiJdCmJldGEzX21jc2UgPC0gZml0X3N1bW1hcnkkc3VtbWFyeVsiYmV0YTMiLCAic2VfbWVhbiJdCnNpZ21hX21jc2UgPC0gZml0X3N1bW1hcnkkc3VtbWFyeVsic2lnbWEiLCAic2VfbWVhbiJdCnNpZ21hX2IwX21jc2UgPC0gZml0X3N1bW1hcnkkc3VtbWFyeVsic2lnbWFfYjAiLCAic2VfbWVhbiJdCnNpZ21hX2IxX21jc2UgPC0gZml0X3N1bW1hcnkkc3VtbWFyeVsic2lnbWFfYjEiLCAic2VfbWVhbiJdCgojIENyZWF0ZSBhIGRhdGEgZnJhbWUgdG8gbmljZWx5IGRpc3BsYXkgdGhlIHJlc3VsdHMKZXN0aW1hdGVzX2RmIDwtIGRhdGEuZnJhbWUoCiAgUGFyYW1ldGVyID0gYygiYmV0YTAiLCAiYmV0YTEiLCAiYmV0YTIiLCAiYmV0YTMiLCAic2lnbWEiLCAic2lnbWFfYjAiLCAic2lnbWFfYjEiKSwKICBFc3RpbWF0ZSA9IGMoYmV0YTBfZXN0aW1hdGUsIGJldGExX2VzdGltYXRlLCBiZXRhMl9lc3RpbWF0ZSwgYmV0YTNfZXN0aW1hdGUsIHNpZ21hX2VzdGltYXRlLCBzaWdtYV9iMF9lc3RpbWF0ZSwgc2lnbWFfYjFfZXN0aW1hdGUpLAogIE1DU0UgPSBjKGJldGEwX21jc2UsIGJldGExX21jc2UsIGJldGEyX21jc2UsIGJldGEzX21jc2UsIHNpZ21hX21jc2UsIHNpZ21hX2IwX21jc2UsIHNpZ21hX2IxX21jc2UpCikKCiMgUHJpbnQgdGhlIGRhdGEgZnJhbWUKcHJpbnQoZXN0aW1hdGVzX2RmKQoKYGBgCgoKCgojVGFzayAyCgpgYGB7cn0KIyBMb2FkIHRoZSByZXF1aXJlZCBsaWJyYXJ5CmxpYnJhcnkocnN0YW4pCgojIFNldCBhIHJhbmRvbSBzZWVkIGZvciByZXByb2R1Y2liaWxpdHkKc2V0LnNlZWQoMTIzNDUpCgojIExvYWQgdGhlIHRhc2syLnJkYSBkYXRhCmxvYWQoInRhc2syLnJkYSIpCgojIFNwbGl0IHRoZSBkYXRhIGludG8geCAoRGVzaWduIEEpIGFuZCB5IChEZXNpZ24gQikKeCA8LSB6WzE6MzBdICAjIE1lYXN1cmVtZW50cyBmcm9tIERlc2lnbiBBCnkgPC0gelszMTo2MF0gIyBNZWFzdXJlbWVudHMgZnJvbSBEZXNpZ24gQgoKIyBQcmVwYXJlIHRoZSBkYXRhIGxpc3QgZm9yIFNUQU4KZGF0YV9saXN0MSA8LSBsaXN0KE4gPSBsZW5ndGgoeCksIHggPSB4LCB5ID0geSkKYGBgCgpNb2RlbCAxCgpgYGB7cn0KIyBEZWZpbmUgdGhlIFNUQU4gbW9kZWwgZm9yIE1vZGVsIDEKbW9kZWwxX2NvZGUgPC0gIgpkYXRhIHsKICBpbnQgTjsgICAgICAgICAvLyBOdW1iZXIgb2Ygb2JzZXJ2YXRpb25zIGluIGVhY2ggZ3JvdXAKICB2ZWN0b3JbTl0geDsgICAvLyBNZWFzdXJlbWVudHMgZnJvbSBEZXNpZ24gQQogIHZlY3RvcltOXSB5OyAgIC8vIE1lYXN1cmVtZW50cyBmcm9tIERlc2lnbiBCCn0KcGFyYW1ldGVycyB7CiAgcmVhbCBtdV94OyAgICAgICAgICAgICAvLyBNZWFuIG9mIERlc2lnbiBBCiAgcmVhbCBtdV95OyAgICAgICAgICAgICAvLyBNZWFuIG9mIERlc2lnbiBCCiAgcmVhbDxsb3dlcj0wPiBzaWdtYTsgICAvLyBDb21tb24gc3RhbmRhcmQgZGV2aWF0aW9uCn0KbW9kZWwgewogIC8vIFByaW9ycwogIG11X3ggfiBub3JtYWwoMCwgMTApOwogIG11X3kgfiBub3JtYWwoMCwgMTApOwogIHNpZ21hIH4gY2F1Y2h5KDAsIDUpOwoKICAvLyBMaWtlbGlob29kCiAgeCB+IG5vcm1hbChtdV94LCBzaWdtYSk7CiAgeSB+IG5vcm1hbChtdV95LCBzaWdtYSk7Cn0KZ2VuZXJhdGVkIHF1YW50aXRpZXMgewogIHJlYWwgcF9tdV95X2dlX211X3ggPSBtdV95ID49IG11X3ggPyAxIDogMDsgIC8vIFAobXVfeSA+PSBtdV94KQogIHJlYWwgZnV0dXJlX3ggPSBub3JtYWxfcm5nKG11X3gsIHNpZ21hKTsgICAgIC8vIEZ1dHVyZSBvYnNlcnZhdGlvbiBmb3IgRGVzaWduIEEKICByZWFsIGZ1dHVyZV95ID0gbm9ybWFsX3JuZyhtdV95LCBzaWdtYSk7ICAgICAvLyBGdXR1cmUgb2JzZXJ2YXRpb24gZm9yIERlc2lnbiBCCiAgcmVhbCBwX2Z1dHVyZV94X2dlXzE1ID0gZnV0dXJlX3ggPj0gMTUgPyAxIDogMDsgLy8gUChmdXR1cmVfeCA+PSAxNSkKICByZWFsIHBfZnV0dXJlX3lfZ2VfMTUgPSBmdXR1cmVfeSA+PSAxNSA/IDEgOiAwOyAvLyBQKGZ1dHVyZV95ID49IDE1KQogIHJlYWwgcF95X2dlX3ggPSBmdXR1cmVfeSA+PSBmdXR1cmVfeCA/IDEgOiAwOyAgIC8vIFAoZnV0dXJlX3kgPj0gZnV0dXJlX3gpCn0KIgoKIyBDb21waWxlIGFuZCBmaXQgdGhlIG1vZGVsIHRvIHRoZSBkYXRhCmZpdDEgPC0gc3Rhbihtb2RlbF9jb2RlID0gbW9kZWwxX2NvZGUsIGRhdGEgPSBkYXRhX2xpc3QxLCBpdGVyID0gMjAwMCwgY2hhaW5zID0gNCkKCiMgUHJpbnQgc3VtbWFyeSBvZiB0aGUgZml0CnByaW50KHN1bW1hcnkoZml0MSkpCgpgYGAKCiNNb250ZSBDYXJsbyBlcnJvciBlc3RpbWF0ZXMgbW9kZWwgMQoKYGBge3J9CiMgQXNzdW1pbmcgJ2ZpdDEnIGlzIHRoZSBmaXR0ZWQgbW9kZWwgb2JqZWN0IGZyb20gTW9kZWwgMQoKIyBFeHRyYWN0IHRoZSBzYW1wbGVzIGZvciByZWxldmFudCBwYXJhbWV0ZXJzCnNhbXBsZXMgPC0gZXh0cmFjdChmaXQxKQoKIyBDYWxjdWxhdGUgTW9udGUgQ2FybG8gZXJyb3JzIGZvciBrZXkgcGFyYW1ldGVycwptY19lcnJvcl9tdV94IDwtIHNkKHNhbXBsZXMkbXVfeCkgLyBzcXJ0KGxlbmd0aChzYW1wbGVzJG11X3gpKQptY19lcnJvcl9tdV95IDwtIHNkKHNhbXBsZXMkbXVfeSkgLyBzcXJ0KGxlbmd0aChzYW1wbGVzJG11X3kpKQptY19lcnJvcl9zaWdtYSA8LSBzZChzYW1wbGVzJHNpZ21hKSAvIHNxcnQobGVuZ3RoKHNhbXBsZXMkc2lnbWEpKQoKIyBDYWxjdWxhdGluZyBNb250ZSBDYXJsbyBlcnJvcnMgZm9yIHRoZSBnZW5lcmF0ZWQgcXVhbnRpdGllcwptY19lcnJvcl9wX211X3lfZ2VfbXVfeCA8LSBzZChzYW1wbGVzJHBfbXVfeV9nZV9tdV94KSAvIHNxcnQobGVuZ3RoKHNhbXBsZXMkcF9tdV95X2dlX211X3gpKQptY19lcnJvcl9wX2Z1dHVyZV94X2dlXzE1IDwtIHNkKHNhbXBsZXMkcF9mdXR1cmVfeF9nZV8xNSkgLyBzcXJ0KGxlbmd0aChzYW1wbGVzJHBfZnV0dXJlX3hfZ2VfMTUpKQptY19lcnJvcl9wX2Z1dHVyZV95X2dlXzE1IDwtIHNkKHNhbXBsZXMkcF9mdXR1cmVfeV9nZV8xNSkgLyBzcXJ0KGxlbmd0aChzYW1wbGVzJHBfZnV0dXJlX3lfZ2VfMTUpKQptY19lcnJvcl9wX3lfZ2VfeCA8LSBzZChzYW1wbGVzJHBfeV9nZV94KSAvIHNxcnQobGVuZ3RoKHNhbXBsZXMkcF95X2dlX3gpKQoKIyBQcmludCB0aGUgTW9udGUgQ2FybG8gZXJyb3JzCmNhdCgiTW9udGUgQ2FybG8gRXJyb3IgZm9yIG11X3g6ICIsIG1jX2Vycm9yX211X3gsICJcbiIpCmNhdCgiTW9udGUgQ2FybG8gRXJyb3IgZm9yIG11X3k6ICIsIG1jX2Vycm9yX211X3ksICJcbiIpCmNhdCgiTW9udGUgQ2FybG8gRXJyb3IgZm9yIHNpZ21hOiAiLCBtY19lcnJvcl9zaWdtYSwgIlxuIikKY2F0KCJNb250ZSBDYXJsbyBFcnJvciBmb3IgUChtdV95ID49IG11X3gpOiAiLCBtY19lcnJvcl9wX211X3lfZ2VfbXVfeCwgIlxuIikKY2F0KCJNb250ZSBDYXJsbyBFcnJvciBmb3IgUChmdXR1cmVfeCA+PSAxNSk6ICIsIG1jX2Vycm9yX3BfZnV0dXJlX3hfZ2VfMTUsICJcbiIpCmNhdCgiTW9udGUgQ2FybG8gRXJyb3IgZm9yIFAoZnV0dXJlX3kgPj0gMTUpOiAiLCBtY19lcnJvcl9wX2Z1dHVyZV95X2dlXzE1LCAiXG4iKQpjYXQoIk1vbnRlIENhcmxvIEVycm9yIGZvciBQKGZ1dHVyZV95ID49IGZ1dHVyZV94KTogIiwgbWNfZXJyb3JfcF95X2dlX3gsICJcbiIpCmBgYAoKCiNNb2RlbCAyCgpgYGB7cn0KIyBMb2FkIHRoZSByZXF1aXJlZCBsaWJyYXJ5CmxpYnJhcnkocnN0YW4pCgojIFNldCBhIHJhbmRvbSBzZWVkIGZvciByZXByb2R1Y2liaWxpdHkKc2V0LnNlZWQoMTIzNDUpCgojIExvYWQgdGhlIHRhc2syLnJkYSBkYXRhCmxvYWQoInRhc2syLnJkYSIpCgojIFNwbGl0IHRoZSBkYXRhIGludG8geCAoRGVzaWduIEEpIGFuZCB5IChEZXNpZ24gQikKeCA8LSB6WzE6MzBdICAjIE1lYXN1cmVtZW50cyBmcm9tIERlc2lnbiBBCnkgPC0gelszMTo2MF0gIyBNZWFzdXJlbWVudHMgZnJvbSBEZXNpZ24gQgoKIyBQcmVwYXJlIHRoZSBkYXRhIGxpc3QgZm9yIFNUQU4KZGF0YV9saXN0MSA8LSBsaXN0KE4gPSBsZW5ndGgoeCksIHggPSB4LCB5ID0geSkKCiMgRGVmaW5lIHRoZSBTVEFOIG1vZGVsIGZvciBNb2RlbCAyCm1vZGVsMl9jb2RlIDwtICIKZGF0YSB7CiAgaW50IE47ICAgICAgICAgLy8gTnVtYmVyIG9mIG9ic2VydmF0aW9ucyBpbiBlYWNoIGdyb3VwCiAgdmVjdG9yW05dIHg7ICAgLy8gTWVhc3VyZW1lbnRzIGZyb20gRGVzaWduIEEKICB2ZWN0b3JbTl0geTsgICAvLyBNZWFzdXJlbWVudHMgZnJvbSBEZXNpZ24gQgp9CnBhcmFtZXRlcnMgewogIHJlYWwgbXVfeDsgICAgICAgICAgICAgLy8gTWVhbiBvZiBEZXNpZ24gQQogIHJlYWwgbXVfeTsgICAgICAgICAgICAgLy8gTWVhbiBvZiBEZXNpZ24gQgogIHJlYWw8bG93ZXI9MD4gc2lnbWFfeDsgLy8gU3RhbmRhcmQgZGV2aWF0aW9uIGZvciBEZXNpZ24gQQogIHJlYWw8bG93ZXI9MD4gc2lnbWFfeTsgLy8gU3RhbmRhcmQgZGV2aWF0aW9uIGZvciBEZXNpZ24gQgp9Cm1vZGVsIHsKICAvLyBQcmlvcnMKICBtdV94IH4gbm9ybWFsKDAsIDEwKTsKICBtdV95IH4gbm9ybWFsKDAsIDEwKTsKICBzaWdtYV94IH4gY2F1Y2h5KDAsIDUpOwogIHNpZ21hX3kgfiBjYXVjaHkoMCwgNSk7CgogIC8vIExpa2VsaWhvb2QKICB4IH4gbm9ybWFsKG11X3gsIHNpZ21hX3gpOwogIHkgfiBub3JtYWwobXVfeSwgc2lnbWFfeSk7Cn0KZ2VuZXJhdGVkIHF1YW50aXRpZXMgewogIHJlYWwgcF9tdV95X2dlX211X3ggPSBtdV95ID49IG11X3ggPyAxIDogMDsgIC8vIFAobXVfeSA+PSBtdV94KQogIHJlYWwgZnV0dXJlX3ggPSBub3JtYWxfcm5nKG11X3gsIHNpZ21hX3gpOyAgIC8vIEZ1dHVyZSBvYnNlcnZhdGlvbiBmb3IgRGVzaWduIEEKICByZWFsIGZ1dHVyZV95ID0gbm9ybWFsX3JuZyhtdV95LCBzaWdtYV95KTsgICAvLyBGdXR1cmUgb2JzZXJ2YXRpb24gZm9yIERlc2lnbiBCCiAgcmVhbCBwX2Z1dHVyZV94X2dlXzE1ID0gZnV0dXJlX3ggPj0gMTUgPyAxIDogMDsgLy8gUChmdXR1cmVfeCA+PSAxNSkKICByZWFsIHBfZnV0dXJlX3lfZ2VfMTUgPSBmdXR1cmVfeSA+PSAxNSA/IDEgOiAwOyAvLyBQKGZ1dHVyZV95ID49IDE1KQogIHJlYWwgcF95X2dlX3ggPSBmdXR1cmVfeSA+PSBmdXR1cmVfeCA/IDEgOiAwOyAgIC8vIFAoZnV0dXJlX3kgPj0gZnV0dXJlX3gpCn0KIgoKIyBDb21waWxlIGFuZCBmaXQgdGhlIG1vZGVsIHRvIHRoZSBkYXRhCmZpdDIgPC0gc3Rhbihtb2RlbF9jb2RlID0gbW9kZWwyX2NvZGUsIGRhdGEgPSBkYXRhX2xpc3QxLCBpdGVyID0gMjAwMCwgY2hhaW5zID0gNCkKCiMgUHJpbnQgc3VtbWFyeSBvZiB0aGUgZml0CnByaW50KHN1bW1hcnkoZml0MikpCmBgYAoKI01vbnRlIENhcmxvIEVycm9yIE1vZGVsIDIKCmBgYHtyfQojIEFzc3VtaW5nICdmaXQyJyBpcyB0aGUgZml0dGVkIG1vZGVsIG9iamVjdCBmcm9tIE1vZGVsIDIKCiMgRXh0cmFjdCB0aGUgc2FtcGxlcyBmb3IgcmVsZXZhbnQgcGFyYW1ldGVycwpzYW1wbGVzMiA8LSBleHRyYWN0KGZpdDIpCgojIENhbGN1bGF0ZSBNb250ZSBDYXJsbyBlcnJvcnMgZm9yIGtleSBwYXJhbWV0ZXJzCm1jX2Vycm9yX211X3gyIDwtIHNkKHNhbXBsZXMyJG11X3gpIC8gc3FydChsZW5ndGgoc2FtcGxlczIkbXVfeCkpCm1jX2Vycm9yX211X3kyIDwtIHNkKHNhbXBsZXMyJG11X3kpIC8gc3FydChsZW5ndGgoc2FtcGxlczIkbXVfeSkpCm1jX2Vycm9yX3NpZ21hX3gyIDwtIHNkKHNhbXBsZXMyJHNpZ21hX3gpIC8gc3FydChsZW5ndGgoc2FtcGxlczIkc2lnbWFfeCkpCm1jX2Vycm9yX3NpZ21hX3kyIDwtIHNkKHNhbXBsZXMyJHNpZ21hX3kpIC8gc3FydChsZW5ndGgoc2FtcGxlczIkc2lnbWFfeSkpCgojIENhbGN1bGF0aW5nIE1vbnRlIENhcmxvIGVycm9ycyBmb3IgdGhlIGdlbmVyYXRlZCBxdWFudGl0aWVzCm1jX2Vycm9yX3BfbXVfeV9nZV9tdV94MiA8LSBzZChzYW1wbGVzMiRwX211X3lfZ2VfbXVfeCkgLyBzcXJ0KGxlbmd0aChzYW1wbGVzMiRwX211X3lfZ2VfbXVfeCkpCm1jX2Vycm9yX3BfZnV0dXJlX3hfZ2VfMTUyIDwtIHNkKHNhbXBsZXMyJHBfZnV0dXJlX3hfZ2VfMTUpIC8gc3FydChsZW5ndGgoc2FtcGxlczIkcF9mdXR1cmVfeF9nZV8xNSkpCm1jX2Vycm9yX3BfZnV0dXJlX3lfZ2VfMTUyIDwtIHNkKHNhbXBsZXMyJHBfZnV0dXJlX3lfZ2VfMTUpIC8gc3FydChsZW5ndGgoc2FtcGxlczIkcF9mdXR1cmVfeV9nZV8xNSkpCm1jX2Vycm9yX3BfeV9nZV94MiA8LSBzZChzYW1wbGVzMiRwX3lfZ2VfeCkgLyBzcXJ0KGxlbmd0aChzYW1wbGVzMiRwX3lfZ2VfeCkpCgojIFByaW50IHRoZSBNb250ZSBDYXJsbyBlcnJvcnMKY2F0KCJNb250ZSBDYXJsbyBFcnJvciBmb3IgbXVfeCAoTW9kZWwgMik6ICIsIG1jX2Vycm9yX211X3gyLCAiXG4iKQpjYXQoIk1vbnRlIENhcmxvIEVycm9yIGZvciBtdV95IChNb2RlbCAyKTogIiwgbWNfZXJyb3JfbXVfeTIsICJcbiIpCmNhdCgiTW9udGUgQ2FybG8gRXJyb3IgZm9yIHNpZ21hX3ggKE1vZGVsIDIpOiAiLCBtY19lcnJvcl9zaWdtYV94MiwgIlxuIikKY2F0KCJNb250ZSBDYXJsbyBFcnJvciBmb3Igc2lnbWFfeSAoTW9kZWwgMik6ICIsIG1jX2Vycm9yX3NpZ21hX3kyLCAiXG4iKQpjYXQoIk1vbnRlIENhcmxvIEVycm9yIGZvciBQKG11X3kgPj0gbXVfeCkgKE1vZGVsIDIpOiAiLCBtY19lcnJvcl9wX211X3lfZ2VfbXVfeDIsICJcbiIpCmNhdCgiTW9udGUgQ2FybG8gRXJyb3IgZm9yIFAoZnV0dXJlX3ggPj0gMTUpIChNb2RlbCAyKTogIiwgbWNfZXJyb3JfcF9mdXR1cmVfeF9nZV8xNTIsICJcbiIpCmNhdCgiTW9udGUgQ2FybG8gRXJyb3IgZm9yIFAoZnV0dXJlX3kgPj0gMTUpIChNb2RlbCAyKTogIiwgbWNfZXJyb3JfcF9mdXR1cmVfeV9nZV8xNTIsICJcbiIpCmNhdCgiTW9udGUgQ2FybG8gRXJyb3IgZm9yIFAoZnV0dXJlX3kgPj0gZnV0dXJlX3gpIChNb2RlbCAyKTogIiwgbWNfZXJyb3JfcF95X2dlX3gyLCAiXG4iKQpgYGAKCgpUYXNrIDMpCgozYSkKCmBgYHtyfQpzZXQuc2VlZCgxMjMpCiMgUmVxdWlyZWQgTGlicmFyaWVzCmxpYnJhcnkocnN0YW4pCgojIExvYWQgRGF0YQpiYiA8LSByZWFkLmNzdigibWV0YS5jc3YiLCBzdHJpbmdzQXNGYWN0b3JzID0gVFJVRSkKCiMgQ3JlYXRlIFN1cnZpdmFscyBWYXJpYWJsZQpiYiA8LSB3aXRoaW4oYmIsIFN1cnZpdmFscyA8LSBUb3RhbCAtIERlYXRocykKCiMgU3VtbWFyaXplZCBUYWJsZQpzdW1tYXJ5X3RhYmxlIDwtIHh0YWJzKGNiaW5kKERlYXRocywgU3Vydml2YWxzKSB+IFRyZWF0bWVudCwgZGF0YT1iYikKcHJpbnQoc3VtbWFyeV90YWJsZSkKCiMgRXh0cmFjdGluZyBEYXRhIGZvciB0aGUgTW9kZWwKbTEgPC0gc3VtKGJiJFRvdGFsW2JiJFRyZWF0bWVudCA9PSAiY29udHJvbCJdKSAgICMgVG90YWwgaW4gY29udHJvbAp5MSA8LSBzdW0oYmIkRGVhdGhzW2JiJFRyZWF0bWVudCA9PSAiY29udHJvbCJdKSAgIyBEZWF0aHMgaW4gY29udHJvbAptMiA8LSBzdW0oYmIkVG90YWxbYmIkVHJlYXRtZW50ID09ICJ0cmVhdGVkIl0pICAgIyBUb3RhbCBpbiB0cmVhdGVkCnkyIDwtIHN1bShiYiREZWF0aHNbYmIkVHJlYXRtZW50ID09ICJ0cmVhdGVkIl0pICAjIERlYXRocyBpbiB0cmVhdGVkCgojIERlZmluZSB0aGUgU3RhbiBNb2RlbApzdGFuX21vZGVsX2NvZGUgPC0gIgpkYXRhIHsKICAgIGludDxsb3dlcj0wPiBtMTsgLy8gVG90YWwgaW4gY29udHJvbAogICAgaW50PGxvd2VyPTA+IHkxOyAvLyBEZWF0aHMgaW4gY29udHJvbAogICAgaW50PGxvd2VyPTA+IG0yOyAvLyBUb3RhbCBpbiB0cmVhdGVkCiAgICBpbnQ8bG93ZXI9MD4geTI7IC8vIERlYXRocyBpbiB0cmVhdGVkCn0KCnBhcmFtZXRlcnMgewogICAgcmVhbDxsb3dlcj0wLCB1cHBlcj0xPiBwMTsgLy8gUHJvYmFiaWxpdHkgb2YgZGVhdGggaW4gY29udHJvbAogICAgcmVhbDxsb3dlcj0wLCB1cHBlcj0xPiBwMjsgLy8gUHJvYmFiaWxpdHkgb2YgZGVhdGggaW4gdHJlYXRlZAp9Cgptb2RlbCB7CiAgICAvLyBVbmlmb3JtIHByaW9ycwogICAgcDEgfiBiZXRhKDEsIDEpOwogICAgcDIgfiBiZXRhKDEsIDEpOwoKICAgIC8vIExpa2VsaWhvb2QKICAgIHkxIH4gYmlub21pYWwobTEsIHAxKTsKICAgIHkyIH4gYmlub21pYWwobTIsIHAyKTsKfQoKZ2VuZXJhdGVkIHF1YW50aXRpZXMgewogICAgcmVhbCBwc2kgPSAocDIgLyAoMSAtIHAyKSkgLyAocDEgLyAoMSAtIHAxKSk7IC8vIE9kZHMgcmF0aW8KfQoiCgojIERhdGEgZm9yIFN0YW4gTW9kZWwKZGF0YV9saXN0IDwtIGxpc3QoCiAgICBtMSA9IG0xLAogICAgeTEgPSB5MSwKICAgIG0yID0gbTIsCiAgICB5MiA9IHkyCikKCiMgTW9kZWwgRml0dGluZwpmaXQgPC0gc3Rhbihtb2RlbF9jb2RlID0gc3Rhbl9tb2RlbF9jb2RlLCBkYXRhID0gZGF0YV9saXN0LCBpdGVyID0gMjAwMCwgY2hhaW5zID0gNCkKcHJpbnQoZml0LCBwYXJzID0gInBzaSIpCgpgYGAKTW9udGUgQ2FybG8gRXJyb3IKCmBgYHtyfQojIFJlcXVpcmVkIExpYnJhcnkKbGlicmFyeShyc3RhbikKCiMgQXNzdW1pbmcgJ2ZpdCcgaXMgdGhlIGZpdHRlZCBtb2RlbCBvYmplY3QgZnJvbSB0aGUgcHJldmlvdXMgU3RhbiBtb2RlbCBydW4KIyBVc2UgdGhlIG1vbml0b3IgZnVuY3Rpb24gdG8gZ2V0IGRpYWdub3N0aWNzLCBpbmNsdWRpbmcgTUNTRQpkaWFnbm9zdGljcyA8LSBtb25pdG9yKGZpdCwgd2FybXVwID0gMCkKCiMgUHJpbnQgdGhlIE1vbnRlIENhcmxvIFN0YW5kYXJkIEVycm9ycwpwcmludChkaWFnbm9zdGljcyRtY3NlKQoKYGBgCgoKCmIpCgooMSkKCmBgYHtyfQpsaWJyYXJ5KHJzdGFuKQoKIyBTdGFuIG1vZGVsIGNvZGUKc3Rhbl9tb2RlbF9jb2RlIDwtICIKZGF0YSB7CiAgICBpbnQ8bG93ZXI9MD4gbjsgICAgICAgIC8vIG51bWJlciBvZiBvYnNlcnZhdGlvbnMKICAgIGludDxsb3dlcj0wPiBkZWF0aHNbbl07IC8vIG51bWJlciBvZiBkZWF0aHMgZm9yIGVhY2ggb2JzZXJ2YXRpb24KICAgIGludDxsb3dlcj0wPiB0b3RhbFtuXTsgIC8vIHRvdGFsIG51bWJlciBvZiBwYXJ0aWNpcGFudHMgZm9yIGVhY2ggb2JzZXJ2YXRpb24KICAgIGludDxsb3dlcj0wLCB1cHBlcj0xPiB0cmVhdG1lbnRbbl07IC8vIHRyZWF0bWVudCBpbmRpY2F0b3IKfQoKcGFyYW1ldGVycyB7CiAgICByZWFsIGJldGEwOyAgICAgICAgLy8gaW50ZXJjZXB0CiAgICByZWFsIGJldGExOyAgICAgICAgLy8gY29lZmZpY2llbnQgZm9yIHRyZWF0bWVudCBlZmZlY3QKfQoKbW9kZWwgewogICAgZm9yIChpIGluIDE6bikgewogICAgICAgIGRlYXRoc1tpXSB+IGJpbm9taWFsX2xvZ2l0KHRvdGFsW2ldLCBiZXRhMCArIGJldGExICogdHJlYXRtZW50W2ldKTsKICAgIH0KICAgIAogICAgLy8gUHJpb3JzCiAgICBiZXRhMCB+IG5vcm1hbCgwLCAxMCk7CiAgICBiZXRhMSB+IG5vcm1hbCgwLCAxMCk7Cn0KCmdlbmVyYXRlZCBxdWFudGl0aWVzIHsKICAgIHJlYWwgb2Rkc19yYXRpbyA9IGV4cChiZXRhMSk7Cn0KIgoKIyBEYXRhIHByZXBhcmF0aW9uIChhc3N1bWluZyBiYiBpcyBhbHJlYWR5IGxvYWRlZCkKYmIgPC0gd2l0aGluKGJiLCB4IDwtIGlmZWxzZShUcmVhdG1lbnQ9PSJ0cmVhdGVkIiwgMSwgMCkpCgojIERhdGEgZm9yIFN0YW4gbW9kZWwKZGF0YV9saXN0IDwtIGxpc3QoCiAgICBuID0gbnJvdyhiYiksCiAgICBkZWF0aHMgPSBiYiREZWF0aHMsCiAgICB0b3RhbCA9IGJiJFRvdGFsLAogICAgdHJlYXRtZW50ID0gYmIkeAopCgojIFNldCByYW5kb20gc2VlZCBmb3IgcmVwcm9kdWNpYmlsaXR5CnNldC5zZWVkKDEyMykKCiMgRml0IHRoZSBtb2RlbApmaXQgPC0gc3Rhbihtb2RlbF9jb2RlID0gc3Rhbl9tb2RlbF9jb2RlLCBkYXRhID0gZGF0YV9saXN0LCBpdGVyID0gMjAwMCwgY2hhaW5zID0gNCkKCiMgUHJpbnQgdGhlIHN1bW1hcnkgZm9yIGJldGExIGFuZCBvZGRzX3JhdGlvCnByaW50KGZpdCwgcGFycyA9IGMoImJldGExIiwgIm9kZHNfcmF0aW8iKSkKYGBgCgpNb250ZSBDYXJsbyBlcnJvcgoKYGBge3J9CiMgUmVxdWlyZWQgTGlicmFyeQpsaWJyYXJ5KHJzdGFuKQoKIyBBc3N1bWluZyAnZml0JyBpcyB0aGUgZml0dGVkIG1vZGVsIG9iamVjdCBmcm9tIHRoZSBwcmV2aW91cyBTdGFuIG1vZGVsIHJ1bgojIFVzZSB0aGUgbW9uaXRvciBmdW5jdGlvbiB0byBnZXQgZGlhZ25vc3RpY3MsIGluY2x1ZGluZyBNQ1NFCmRpYWdub3N0aWNzIDwtIG1vbml0b3IoZml0LCB3YXJtdXAgPSAwKQoKIyBQcmludCB0aGUgTW9udGUgQ2FybG8gU3RhbmRhcmQgRXJyb3JzCnByaW50KGRpYWdub3N0aWNzJG1jc2UpCmBgYAoKCgooMikKCmBgYHtyfQojIExvYWQgUlN0YW4gbGlicmFyeQpsaWJyYXJ5KHJzdGFuKQoKIyBTZXQgcmFuZG9tIHNlZWQgZm9yIHJlcHJvZHVjaWJpbGl0eQpzZXQuc2VlZCgxMjM0NSkKCiMgTG9hZCBkYXRhIApiYiA8LSByZWFkLmNzdigibWV0YS5jc3YiLCBzdHJpbmdzQXNGYWN0b3JzID0gVFJVRSkKCiMgRGVmaW5lIHRoZSBpbmRpY2F0b3IgdmFyaWFibGUgZm9yIFRyZWF0bWVudApiYiA8LSB3aXRoaW4oYmIsIHggPC0gaWZlbHNlKFRyZWF0bWVudD09InRyZWF0ZWQiLCAxLCAwKSkKCiMgUHJlcGFyZSBkYXRhIGZvciBTdGFuCmRhdGFfbGlzdCA8LSBsaXN0KAogIE4gPSBucm93KGJiKSwKICBKID0gbGVuZ3RoKHVuaXF1ZShiYiRTdHVkeSkpLCAgIyBUb3RhbCBudW1iZXIgb2Ygc3R1ZGllcwogIGdyb3VwID0gYXMubnVtZXJpYyhmYWN0b3IoYmIkU3R1ZHkpKSwgICMgQ29udmVydCBzdHVkeSBpZGVudGlmaWVycyB0byBudW1lcmljIGdyb3VwIGluZGV4CiAgeSA9IGJiJERlYXRocywKICBtID0gYmIkVG90YWwsCiAgeCA9IGJiJHgKKQoKIyBTdGFuIG1vZGVsCnN0YW5fbW9kZWwgPC0gIgpkYXRhIHsKICBpbnQ8bG93ZXI9MD4gTjsgIC8vIFRvdGFsIG51bWJlciBvZiBvYnNlcnZhdGlvbnMKICBpbnQ8bG93ZXI9MD4gSjsgIC8vIFRvdGFsIG51bWJlciBvZiBncm91cHMKICBpbnQ8bG93ZXI9MSwgdXBwZXI9Sj4gZ3JvdXBbTl07ICAvLyBHcm91cCBpbmRleCBmb3IgZWFjaCBvYnNlcnZhdGlvbgogIGludDxsb3dlcj0wPiB5W05dOyAgLy8gUmVzcG9uc2UgdmFyaWFibGUgKG51bWJlciBvZiBkZWF0aHMpCiAgaW50PGxvd2VyPTA+IG1bTl07ICAvLyBUb3RhbCBwYXJ0aWNpcGFudHMgZm9yIGVhY2ggb2JzZXJ2YXRpb24KICBpbnQ8bG93ZXI9MCwgdXBwZXI9MT4geFtOXTsgIC8vIFRyZWF0bWVudCBpbmRpY2F0b3IKfQoKcGFyYW1ldGVycyB7CiAgcmVhbCBiZXRhMTsgIC8vIEdsb2JhbCB0cmVhdG1lbnQgZWZmZWN0CiAgcmVhbCBtdV9iZXRhMDsgIC8vIEdsb2JhbCBpbnRlcmNlcHQKICByZWFsPGxvd2VyPTA+IHNpZ21hX2JldGEwOyAgLy8gU3RhbmRhcmQgZGV2aWF0aW9uIG9mIGdyb3VwLXNwZWNpZmljIGludGVyY2VwdHMKICB2ZWN0b3JbSl0gYmV0YTBfcmF3OyAgLy8gU3RhbmRhcmRpemVkIGdyb3VwLXNwZWNpZmljIGludGVyY2VwdHMKfQoKdHJhbnNmb3JtZWQgcGFyYW1ldGVycyB7CiAgdmVjdG9yW0pdIGJldGEwOyAgLy8gR3JvdXAtc3BlY2lmaWMgaW50ZXJjZXB0cwogIGJldGEwID0gbXVfYmV0YTAgKyBzaWdtYV9iZXRhMCAqIGJldGEwX3JhdzsgIC8vIE5vbi1jZW50ZXJlZCBwYXJhbWV0ZXJpemF0aW9uCn0KCm1vZGVsIHsKICAvLyBQcmlvcnMKICBiZXRhMSB+IG5vcm1hbCgwLCAxMCk7CiAgbXVfYmV0YTAgfiBub3JtYWwoMCwgMTApOwogIHNpZ21hX2JldGEwIH4gdW5pZm9ybSgwLCAxMCk7CiAgYmV0YTBfcmF3IH4gbm9ybWFsKDAsIDEpOwogIAogIC8vIExpa2VsaWhvb2QKICBmb3IgKGkgaW4gMTpOKSB7CiAgICB5W2ldIH4gYmlub21pYWxfbG9naXQobVtpXSwgYmV0YTBbZ3JvdXBbaV1dICsgYmV0YTEgKiB4W2ldKTsKICB9Cn0KCmdlbmVyYXRlZCBxdWFudGl0aWVzIHsKICByZWFsIG9kZHNfcmF0aW87CiAgb2Rkc19yYXRpbyA9IGV4cChiZXRhMSk7ICAvLyBPZGRzIHJhdGlvIGZvciB0cmVhdG1lbnQgZWZmZWN0Cn0KIgoKIyBGaXQgdGhlIG1vZGVsCmZpdCA8LSBzdGFuKG1vZGVsX2NvZGUgPSBzdGFuX21vZGVsLCBkYXRhID0gZGF0YV9saXN0LCBpdGVyID0gMjAwMCwgY2hhaW5zID0gNCkKCiMgRXh0cmFjdCBhbmQgY29tcHV0ZSB0aGUgQmF5ZXNpYW4gZXN0aW1hdGUgZm9yIHRoZSBvZGRzIHJhdGlvCm9kZHNfcmF0aW8gPC0gZXh0cmFjdChmaXQpJG9kZHNfcmF0aW8KYmF5ZXNpYW5fZXN0aW1hdGUgPC0gbWVhbihvZGRzX3JhdGlvKQoKIyBQcmludCB0aGUgQmF5ZXNpYW4gZXN0aW1hdGUgZm9yIHRoZSBvZGRzIHJhdGlvCmNhdCgiQmF5ZXNpYW4gZXN0aW1hdGUgZm9yIHRoZSBvZGRzIHJhdGlvIChlXmJldGExKToiLCBiYXllc2lhbl9lc3RpbWF0ZSwgIlxuIikKCiMgRXh0cmFjdCB0aGUgc2FtcGxlcyBmb3IgYmV0YTEKYmV0YTFfc2FtcGxlcyA8LSBleHRyYWN0KGZpdCkkYmV0YTEKCiMgQ2FsY3VsYXRlIHRoZSBvZGRzIHJhdGlvIGZvciBlYWNoIHNhbXBsZQpvZGRzX3JhdGlvX3NhbXBsZXMgPC0gZXhwKGJldGExX3NhbXBsZXMpCgojIENhbGN1bGF0ZSB0aGUgTW9udGUgQ2FybG8gZXJyb3IgZm9yIHRoZSBvZGRzIHJhdGlvCm1jX2Vycm9yIDwtIHNkKG9kZHNfcmF0aW9fc2FtcGxlcykgLyBzcXJ0KGxlbmd0aChvZGRzX3JhdGlvX3NhbXBsZXMpKQoKIyBQcmludCB0aGUgTW9udGUgQ2FybG8gZXJyb3IKY2F0KCJNb250ZSBDYXJsbyBlcnJvciBmb3IgdGhlIG9kZHMgcmF0aW8gZXN0aW1hdGU6IiwgbWNfZXJyb3IsICJcbiIpCgpgYGAKCgoKCigzKQoKaS4KCmBgYHtyfQojIExvYWQgUlN0YW4gbGlicmFyeQpsaWJyYXJ5KHJzdGFuKQoKIyBTZXQgcmFuZG9tIHNlZWQgZm9yIHJlcHJvZHVjaWJpbGl0eQpzZXQuc2VlZCgxMjM0NSkKCiMgTG9hZCB5b3VyIGRhdGEKYmIgPC0gcmVhZC5jc3YoIm1ldGEuY3N2Iiwgc3RyaW5nc0FzRmFjdG9ycyA9IFRSVUUpCiMgRGVmaW5lIHRoZSBpbmRpY2F0b3IgdmFyaWFibGUgZm9yIFRyZWF0bWVudApiYiA8LSB3aXRoaW4oYmIsIHggPC0gaWZlbHNlKFRyZWF0bWVudCA9PSAidHJlYXRlZCIsIDEsIDApKQoKIyBQcmVwYXJlIGRhdGEgZm9yIFN0YW4KZGF0YV9saXN0IDwtIGxpc3QoCiAgTiA9IG5yb3coYmIpLAogIEogPSBsZW5ndGgodW5pcXVlKGJiJFN0dWR5KSksCiAgZ3JvdXAgPSBhcy5udW1lcmljKGZhY3RvcihiYiRTdHVkeSkpLAogIHkgPSBiYiREZWF0aHMsCiAgbSA9IGJiJFRvdGFsLAogIHggPSBiYiR4CikKCiMgU3RhbiBtb2RlbCBzdHJpbmcKc3Rhbl9tb2RlbCA8LSAiCmRhdGEgewogIC8vIERhdGEgYmxvY2sgZGVmaW5lcyB0aGUgZGF0YSB0eXBlcyBhbmQgdGhlaXIgY29uc3RyYWludHMKICBpbnQ8bG93ZXI9MD4gTjsgIC8vIFRvdGFsIG51bWJlciBvZiBvYnNlcnZhdGlvbnMKICBpbnQ8bG93ZXI9MD4gSjsgIC8vIFRvdGFsIG51bWJlciBvZiBncm91cHMvc3R1ZGllcwogIGludDxsb3dlcj0xLCB1cHBlcj1KPiBncm91cFtOXTsgIC8vIEdyb3VwIGluZGV4IGZvciBlYWNoIG9ic2VydmF0aW9uLCBjb25zdHJhaW5lZCBiZXR3ZWVuIDEgYW5kIEoKICBpbnQ8bG93ZXI9MD4geVtOXTsgIC8vIE51bWJlciBvZiBkZWF0aHMgKHJlc3BvbnNlIHZhcmlhYmxlKSBmb3IgZWFjaCBvYnNlcnZhdGlvbgogIGludDxsb3dlcj0wPiBtW05dOyAgLy8gVG90YWwgbnVtYmVyIG9mIHBhcnRpY2lwYW50cyBmb3IgZWFjaCBvYnNlcnZhdGlvbgogIGludDxsb3dlcj0wLCB1cHBlcj0xPiB4W05dOyAgLy8gVHJlYXRtZW50IGluZGljYXRvciAoMCBvciAxKSBmb3IgZWFjaCBvYnNlcnZhdGlvbgp9CgpwYXJhbWV0ZXJzIHsKICAvLyBQYXJhbWV0ZXJzIGJsb2NrIGRlZmluZXMgdGhlIG1vZGVsIHBhcmFtZXRlcnMKICByZWFsIG11X2JldGEwOyAgLy8gR2xvYmFsIG1lYW4gZm9yIHRoZSBpbnRlcmNlcHRzCiAgcmVhbDxsb3dlcj0wPiBzaWdtYV9iZXRhMDsgIC8vIFN0YW5kYXJkIGRldmlhdGlvbiBmb3IgdGhlIGdyb3VwLXNwZWNpZmljIGludGVyY2VwdHMgKG11c3QgYmUgcG9zaXRpdmUpCiAgdmVjdG9yW0pdIGJldGEwX3JhdzsgIC8vIFJhdywgbm9uLWNlbnRlcmVkIGdyb3VwLXNwZWNpZmljIGludGVyY2VwdHMKCiAgcmVhbCBtdV9iZXRhMTsgIC8vIEdsb2JhbCBtZWFuIGZvciB0aGUgdHJlYXRtZW50IGVmZmVjdHMKICByZWFsPGxvd2VyPTA+IHNpZ21hX2JldGExOyAgLy8gU3RhbmRhcmQgZGV2aWF0aW9uIGZvciB0aGUgZ3JvdXAtc3BlY2lmaWMgdHJlYXRtZW50IGVmZmVjdHMgKG11c3QgYmUgcG9zaXRpdmUpCiAgdmVjdG9yW0pdIGJldGExX3JhdzsgIC8vIFJhdywgbm9uLWNlbnRlcmVkIGdyb3VwLXNwZWNpZmljIHRyZWF0bWVudCBlZmZlY3RzCn0KCnRyYW5zZm9ybWVkIHBhcmFtZXRlcnMgewogIC8vIFRyYW5zZm9ybWVkIHBhcmFtZXRlcnMgYmxvY2sgZm9yIGludGVybWVkaWF0ZSBjYWxjdWxhdGlvbnMKICB2ZWN0b3JbSl0gYmV0YTA7ICAvLyBBY3R1YWwgZ3JvdXAtc3BlY2lmaWMgaW50ZXJjZXB0cwogIHZlY3RvcltKXSBiZXRhMTsgIC8vIEFjdHVhbCBncm91cC1zcGVjaWZpYyB0cmVhdG1lbnQgZWZmZWN0cwogIGJldGEwID0gbXVfYmV0YTAgKyBzaWdtYV9iZXRhMCAqIGJldGEwX3JhdzsgIC8vIENlbnRlcmluZyB0aGUgZ3JvdXAtc3BlY2lmaWMgaW50ZXJjZXB0cwogIGJldGExID0gbXVfYmV0YTEgKyBzaWdtYV9iZXRhMSAqIGJldGExX3JhdzsgIC8vIENlbnRlcmluZyB0aGUgZ3JvdXAtc3BlY2lmaWMgdHJlYXRtZW50IGVmZmVjdHMKfQoKbW9kZWwgewogIC8vIE1vZGVsIGJsb2NrIGRlZmluZXMgdGhlIGxpa2VsaWhvb2QgYW5kIHByaW9ycwogIC8vIFByaW9ycyBmb3IgZ2xvYmFsIG1lYW5zIGFuZCBzdGFuZGFyZCBkZXZpYXRpb25zCiAgbXVfYmV0YTAgfiBub3JtYWwoMCwgMTApOyAgLy8gV2Vha2x5IGluZm9ybWF0aXZlIG5vcm1hbCBwcmlvciBmb3IgZ2xvYmFsIG1lYW4gaW50ZXJjZXB0CiAgc2lnbWFfYmV0YTAgfiB1bmlmb3JtKDAsIDEwKTsgIC8vIFVuaWZvcm0gcHJpb3IgZm9yIHN0YW5kYXJkIGRldmlhdGlvbiBvZiBpbnRlcmNlcHRzCiAgYmV0YTBfcmF3IH4gbm9ybWFsKDAsIDEpOyAgLy8gU3RhbmRhcmQgbm9ybWFsIHByaW9yIGZvciByYXcgaW50ZXJjZXB0cwoKICBtdV9iZXRhMSB+IG5vcm1hbCgwLCAxMCk7ICAvLyBXZWFrbHkgaW5mb3JtYXRpdmUgbm9ybWFsIHByaW9yIGZvciBnbG9iYWwgbWVhbiB0cmVhdG1lbnQgZWZmZWN0CiAgc2lnbWFfYmV0YTEgfiB1bmlmb3JtKDAsIDEwKTsgIC8vIFVuaWZvcm0gcHJpb3IgZm9yIHN0YW5kYXJkIGRldmlhdGlvbiBvZiB0cmVhdG1lbnQgZWZmZWN0cwogIGJldGExX3JhdyB+IG5vcm1hbCgwLCAxKTsgIC8vIFN0YW5kYXJkIG5vcm1hbCBwcmlvciBmb3IgcmF3IHRyZWF0bWVudCBlZmZlY3RzCgogIC8vIExpa2VsaWhvb2Qgb2YgdGhlIGRhdGEgZ2l2ZW4gdGhlIG1vZGVsCiAgZm9yIChpIGluIDE6TikgewogICAgeVtpXSB+IGJpbm9taWFsX2xvZ2l0KG1baV0sIGJldGEwW2dyb3VwW2ldXSArIGJldGExW2dyb3VwW2ldXSAqIHhbaV0pOyAgLy8gQmlub21pYWwgbG9naXQgbW9kZWwgZm9yIGVhY2ggb2JzZXJ2YXRpb24KICB9Cn0KCmdlbmVyYXRlZCBxdWFudGl0aWVzIHsKICAvLyBHZW5lcmF0ZWQgcXVhbnRpdGllcyBibG9jayBmb3IgcG9zdC1wcm9jZXNzaW5nCiAgcmVhbCBleHBfbXVfYmV0YTEgPSBleHAobXVfYmV0YTEpOyAgLy8gRXhwb25lbnRpYWwgb2YgbXVfYmV0YTEgdG8gZ2V0IHRoZSBvZGRzIHJhdGlvCn0KIgoKIyBGaXQgdGhlIG1vZGVsIHVzaW5nIFN0YW4KZml0IDwtIHN0YW4obW9kZWxfY29kZSA9IHN0YW5fbW9kZWwsIGRhdGEgPSBkYXRhX2xpc3QsIGl0ZXIgPSAyMDAwLCBjaGFpbnMgPSA0KQoKIyBFeHRyYWN0IHRoZSBwb3N0ZXJpb3Igc2FtcGxlcyBmb3IgZXhwKG11X2JldGExKQpleHBfbXVfYmV0YTFfc2FtcGxlcyA8LSBleHRyYWN0KGZpdCkkZXhwX211X2JldGExCgojIENhbGN1bGF0ZSB0aGUgQmF5ZXNpYW4gZXN0aW1hdGUgZm9yIGV4cChtdV9iZXRhMSkKYmF5ZXNpYW5fZXN0aW1hdGVfZXhwX211X2JldGExIDwtIG1lYW4oZXhwX211X2JldGExX3NhbXBsZXMpCgojIFByaW50IHRoZSBCYXllc2lhbiBlc3RpbWF0ZSBmb3IgZXhwKG11X2JldGExKQpjYXQoIkJheWVzaWFuIGVzdGltYXRlIGZvciBleHAobXVfYmV0YTEpOiIsIGJheWVzaWFuX2VzdGltYXRlX2V4cF9tdV9iZXRhMSwgIlxuIikKCiMgQ2FsY3VsYXRlIHRoZSBzdGFuZGFyZCBkZXZpYXRpb24gb2YgdGhlIHNhbXBsZWQgcG9zdGVyaW9yIHZhbHVlcwpzdGRfZGV2X2V4cF9tdV9iZXRhMSA8LSBzZChleHBfbXVfYmV0YTFfc2FtcGxlcykKCiMgQ2FsY3VsYXRlIHRoZSBudW1iZXIgb2Ygc2FtcGxlcwpudW1fc2FtcGxlcyA8LSBsZW5ndGgoZXhwX211X2JldGExX3NhbXBsZXMpCgojIENhbGN1bGF0ZSB0aGUgTW9udGUgQ2FybG8gZXJyb3IKbW9udGVfY2FybG9fZXJyb3JfZXhwX211X2JldGExIDwtIHN0ZF9kZXZfZXhwX211X2JldGExIC8gc3FydChudW1fc2FtcGxlcykKCiMgUHJpbnQgdGhlIE1vbnRlIENhcmxvIGVycm9yIGZvciBleHAobXVfYmV0YTEpCmNhdCgiTW9udGUgQ2FybG8gZXJyb3IgZm9yIGV4cChtdV9iZXRhMSk6IiwgbW9udGVfY2FybG9fZXJyb3JfZXhwX211X2JldGExLCAiXG4iKQoKCmBgYAoKYGBge3J9CgpsaWJyYXJ5KGdncGxvdDIpIAoKIyBFeHRyYWN0IG1vZGVsIHBhcmFtZXRlcnMKcG9zdGVyaW9yX3NhbXBsZXMgPC0gZXh0cmFjdChmaXQpCm11X2JldGExX3NhbXBsZXMgPC0gcG9zdGVyaW9yX3NhbXBsZXMkbXVfYmV0YTEKbXVfYmV0YTBfc2FtcGxlcyA8LSBwb3N0ZXJpb3Jfc2FtcGxlcyRtdV9iZXRhMApzaWdtYV9iZXRhMV9zYW1wbGVzIDwtIHBvc3Rlcmlvcl9zYW1wbGVzJHNpZ21hX2JldGExCgojIFNpbXVsYXRlIGV4cGVjdGVkIG51bWJlciBvZiBkZWF0aHMgZm9yIGEgbmV3IHN0dWR5IHdpdGggMTAwIHBhcnRpY2lwYW50cyBpbiBlYWNoIGdyb3VwCm5fcGFydGljaXBhbnRzIDwtIDEwMApuX3NpbXVsYXRpb25zIDwtIGxlbmd0aChtdV9iZXRhMV9zYW1wbGVzKQoKIyBJbml0aWFsaXplIHZlY3RvcnMgdG8gc3RvcmUgc2ltdWxhdGVkIGRlYXRocwpzaW11bGF0ZWRfZGVhdGhzX2NvbnRyb2wgPC0gbnVtZXJpYyhuX3NpbXVsYXRpb25zKQpzaW11bGF0ZWRfZGVhdGhzX3RyZWF0ZWQgPC0gbnVtZXJpYyhuX3NpbXVsYXRpb25zKQoKZm9yIChpIGluIDE6bl9zaW11bGF0aW9ucykgewogICMgQ29tcHV0ZSBwcm9iYWJpbGl0aWVzCiAgcHJvYl9jb250cm9sIDwtIHBsb2dpcyhtdV9iZXRhMF9zYW1wbGVzW2ldKQogIHByb2JfdHJlYXRlZCA8LSBwbG9naXMobXVfYmV0YTBfc2FtcGxlc1tpXSArIG11X2JldGExX3NhbXBsZXNbaV0pCgogICMgU2ltdWxhdGUgZGVhdGhzCiAgc2ltdWxhdGVkX2RlYXRoc19jb250cm9sW2ldIDwtIHJiaW5vbSgxLCBuX3BhcnRpY2lwYW50cywgcHJvYl9jb250cm9sKQogIHNpbXVsYXRlZF9kZWF0aHNfdHJlYXRlZFtpXSA8LSByYmlub20oMSwgbl9wYXJ0aWNpcGFudHMsIHByb2JfdHJlYXRlZCkKfQoKIyBDb21wdXRlIGF2ZXJhZ2Ugc2ltdWxhdGVkIGRlYXRocwphdmdfZGVhdGhzX2NvbnRyb2wgPC0gbWVhbihzaW11bGF0ZWRfZGVhdGhzX2NvbnRyb2wpCmF2Z19kZWF0aHNfdHJlYXRlZCA8LSBtZWFuKHNpbXVsYXRlZF9kZWF0aHNfdHJlYXRlZCkKCmNhdCgiQXZlcmFnZSBzaW11bGF0ZWQgZGVhdGhzIGluIGNvbnRyb2wgZ3JvdXA6IiwgYXZnX2RlYXRoc19jb250cm9sLCAiXG4iKQpjYXQoIkF2ZXJhZ2Ugc2ltdWxhdGVkIGRlYXRocyBpbiB0cmVhdGVkIGdyb3VwOiIsIGF2Z19kZWF0aHNfdHJlYXRlZCwgIlxuIikKCiMgR2VuZXJhdGUgcGxvdCBvZiBwcm9iYWJpbGl0eSBtYXNzIGZ1bmN0aW9uIGZvciB0aGUgbnVtYmVyIG9mIGRlYXRocyBpbiB0aGUgdHJlYXRlZCBncm91cApkZl9wbG90IDwtIGRhdGEuZnJhbWUoZGVhdGhzID0gc2ltdWxhdGVkX2RlYXRoc190cmVhdGVkKQpnZ3Bsb3QoZGZfcGxvdCwgYWVzKHggPSBkZWF0aHMpKSArCiAgZ2VvbV9oaXN0b2dyYW0oYmlud2lkdGggPSAxLCBib3VuZGFyeSA9IDAsIGZpbGwgPSAiYmx1ZSIsIGNvbG9yID0gImJsYWNrIikgKwogIGxhYnModGl0bGUgPSAiUHJvYmFiaWxpdHkgTWFzcyBGdW5jdGlvbiBmb3IgTnVtYmVyIG9mIERlYXRocyBpbiBUcmVhdGVkIEdyb3VwIiwKICAgICAgIHggPSAiTnVtYmVyIG9mIERlYXRocyIsCiAgICAgICB5ID0gIkZyZXF1ZW5jeSIpICsKICB0aGVtZV9taW5pbWFsKCkKCmBgYAoKCg==