ANLY505 - Intro to Stan

knitr::opts_chunk$set(echo = TRUE)
suppressPackageStartupMessages(library(rstan))
suppressPackageStartupMessages(library(gdata))
suppressPackageStartupMessages(library(bayesplot))
suppressPackageStartupMessages(library(parallel))
suppressPackageStartupMessages(library(modelr))
suppressPackageStartupMessages(library(dplyr))
suppressPackageStartupMessages(library(tidyverse))
suppressPackageStartupMessages(library(broom))

Intro to STAN Homework

After our Intro to Stan lecture I think it would be valuable to have you go through a similar exercise. Let’s test a second research question.

Research question: Is sea ice extent declining in the Southern Hemisphere over time? Is the same pattern happening in the Antarctic as in the Arctic? Fit a Stan model to find out!

Make sure you follow the steps we used in class.

1. Load and Inspect Data

#place the code here
SH<-read.csv("/Users/Morgane/Documents/Harrisburg University/Analytics Master Classes/ANLY505-90-O FALL Modelling Simulation and Game Theory/Intro to stan/seaice.csv")


summary(SH)

##       year       extent_north    extent_south  
##  Min.   :1979   Min.   :10.15   Min.   :10.69  
##  1st Qu.:1988   1st Qu.:10.90   1st Qu.:11.42  
##  Median :1998   Median :11.67   Median :11.67  
##  Mean   :1998   Mean   :11.46   Mean   :11.68  
##  3rd Qu.:2008   3rd Qu.:11.98   3rd Qu.:11.88  
##  Max.   :2017   Max.   :12.45   Max.   :12.78

head(SH)

##   year extent_north extent_south
## 1 1979       12.328       11.700
## 2 1980       12.337       11.230
## 3 1981       12.127       11.435
## 4 1982       12.447       11.640
## 5 1983       12.332       11.389
## 6 1984       11.910       11.454

2. Plot the data

#plot data

plot(extent_south~year,pch=20,data=SH)

3. Run a general linear model using lm()

#write the code

SouthModel1<-lm(extent_south~year,data=SH)
summary(SouthModel1)

## 
## Call:
## lm(formula = extent_south ~ year, data = SH)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.23372 -0.18142  0.01587  0.18465  0.88814 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -14.199551  10.925576  -1.300   0.2018  
## year          0.012953   0.005468   2.369   0.0232 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3843 on 37 degrees of freedom
## Multiple R-squared:  0.1317, Adjusted R-squared:  0.1082 
## F-statistic: 5.611 on 1 and 37 DF,  p-value: 0.02318

grid<-SH %>%
  data_grid(year)%>%
  add_predictions(SouthModel1)

ggplot(SH, aes(year))+
  geom_point(aes(y=extent_south))+
  geom_line(aes(y=pred),data=grid,color="red",linetype="dashed",size=1)

4. Index the data, re-run the lm(), extract summary statistics and turn the indexed data into a dataframe to pass into Stan

#write the code here

X<-I(SH$year - 1978)
Y<-SH$extent_south
N<-length(SH$year)

SouthModel2<-lm(Y~X)
summary(SouthModel2)

## 
## Call:
## lm(formula = Y ~ X)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.23372 -0.18142  0.01587  0.18465  0.88814 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 11.421555   0.125490  91.015   <2e-16 ***
## X            0.012953   0.005468   2.369   0.0232 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3843 on 37 degrees of freedom
## Multiple R-squared:  0.1317, Adjusted R-squared:  0.1082 
## F-statistic: 5.611 on 1 and 37 DF,  p-value: 0.02318

tidy(SouthModel2)

## # A tibble: 2 x 5
##   term        estimate std.error statistic  p.value
##   <chr>          <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)  11.4      0.125       91.0  4.02e-45
## 2 X             0.0130   0.00547      2.37 2.32e- 2

sigma<-summary(SouthModel2)$sigma
alpha<-summary(SouthModel2)$coeff[1]
beta<-summary(SouthModel2)$coeff[2]

sigma

## [1] 0.384331

alpha

## [1] 11.42155

beta

## [1] 0.01295304

stand_Southdata<-list(N=N,x=X,y=Y)

5. Write the Stan model

#write the code

 write("// Stan model for simple linear regression

data {
 int < lower = 1 > N; // Sample size
 vector[N] x; // Predictor
 vector[N] y; // Outcome
}

parameters {
 real alpha; // Intercept
 real beta; // Slope (regression coefficients)
 real < lower = 0 > sigma; // Error SD
}

model {
 y ~ normal(alpha + x * beta , sigma);
}

generated quantities {
} // The posterior predictive distribution",

"stan_model1.stan")


stan_model1 <- "stan_model1.stan"

6. Check to see how many cores your computer has and enable parallel computing

detectCores(all.tests = FALSE, logical = TRUE)

options(mc.cores = parallel::detectCores())

My computer has 4 cores!

7. Run the Stan model and inspect the results

#code here
fit<-stan(file=stan_model1,data=stand_Southdata,warmup=500,iter=1000,chains=4, cores=2,thin=1)
summary(fit)

## $summary
##              mean      se_mean          sd         2.5%          25%
## alpha 11.42777985 0.0046501294 0.126731755 11.175952447 11.341001412
## beta   0.01275752 0.0001883416 0.005569808  0.001608932  0.009088834
## sigma  0.39797019 0.0015103000 0.047125905  0.320277939  0.364574567
## lp__  16.30369890 0.0440279829 1.224878472 13.148042147 15.690259848
##               50%         75%       97.5%    n_eff     Rhat
## alpha 11.42783947 11.51216500 11.67229212 742.7468 1.007748
## beta   0.01282349  0.01644804  0.02412476 874.5571 1.005539
## sigma  0.39319032  0.42723240  0.49982930 973.6279 1.001542
## lp__  16.62081540 17.23246066 17.74191388 773.9777 1.002746
## 
## $c_summary
## , , chains = chain:1
## 
##          stats
## parameter        mean          sd         2.5%          25%         50%
##     alpha 11.45046942 0.121196256 11.231324111 11.366885527 11.45640197
##     beta   0.01184779 0.005390519  0.001129554  0.007947245  0.01220055
##     sigma  0.39997270 0.049026474  0.321396192  0.363640806  0.39373096
##     lp__  16.26907227 1.286965259 13.124708760 15.587607522 16.60171891
##          stats
## parameter         75%       97.5%
##     alpha 11.53344371 11.66731728
##     beta   0.01571171  0.02327274
##     sigma  0.43191140  0.50293841
##     lp__  17.24077574 17.76622946
## 
## , , chains = chain:2
## 
##          stats
## parameter       mean          sd         2.5%          25%         50%
##     alpha 11.4134484 0.129211285 11.177590175 11.325962599 11.41384719
##     beta   0.0131738 0.005545895  0.001975375  0.009406423  0.01307623
##     sigma  0.3984162 0.045423721  0.323494087  0.368543612  0.39361171
##     lp__  16.3455789 1.229415480 13.010779879 15.780386106 16.65533038
##          stats
## parameter         75%       97.5%
##     alpha 11.50042820 11.68223834
##     beta   0.01698073  0.02324555
##     sigma  0.42494548  0.49803794
##     lp__  17.28355204 17.73255504
## 
## , , chains = chain:3
## 
##          stats
## parameter        mean          sd         2.5%          25%         50%
##     alpha 11.42605838 0.121758773 11.176581425 11.349299127 11.42420125
##     beta   0.01299308 0.005548193  0.001826082  0.009271429  0.01316476
##     sigma  0.39455130 0.045460325  0.320013338  0.362519821  0.39131282
##     lp__  16.38446159 1.134581386 13.644918393 15.768048851 16.71148512
##          stats
## parameter         75%       97.5%
##     alpha 11.49784788 11.67024895
##     beta   0.01648361  0.02492298
##     sigma  0.42056674  0.49894561
##     lp__  17.22171447 17.70993558
## 
## , , chains = chain:4
## 
##          stats
## parameter       mean          sd         2.5%          25%         50%
##     alpha 11.4211432 0.131765963 11.160469610 11.329423051 11.42460322
##     beta   0.0130154 0.005705728  0.002058793  0.009385636  0.01275598
##     sigma  0.3989405 0.048438950  0.318286772  0.365853118  0.39324845
##     lp__  16.2156828 1.240183558 13.226813661 15.609484444 16.47796252
##          stats
## parameter         75%      97.5%
##     alpha 11.51247730 11.6724422
##     beta   0.01654957  0.0252751
##     sigma  0.43334723  0.4962065
##     lp__  17.12865250 17.7372495

8. Extract and inspect the posterior estimates into a list so we can plot them

#code here / I am specifying that it is the rstan package because it seems that there is a conflict with another package otherwise.
posterior_estimate <- rstan::extract(fit)
str(posterior_estimate)

## List of 4
##  $ alpha: num [1:2000(1d)] 11.5 11.6 11.5 11.5 11.5 ...
##   ..- attr(*, "dimnames")=List of 1
##   .. ..$ iterations: NULL
##  $ beta : num [1:2000(1d)] 0.0079 0.01022 0.00356 0.00316 0.00661 ...
##   ..- attr(*, "dimnames")=List of 1
##   .. ..$ iterations: NULL
##  $ sigma: num [1:2000(1d)] 0.321 0.475 0.432 0.361 0.35 ...
##   ..- attr(*, "dimnames")=List of 1
##   .. ..$ iterations: NULL
##  $ lp__ : num [1:2000(1d)] 16 15.6 15.3 15 16.8 ...
##   ..- attr(*, "dimnames")=List of 1
##   .. ..$ iterations: NULL

9. Compare your results to our results to “lm”

plot(y ~ x, pch = 20, data = stand_Southdata)
abline(SouthModel2, col = 2, lty = 2, lw = 3)
abline(mean(posterior_estimate$alpha), mean(posterior_estimate$beta), col = 4, lw = 2)

10. Plot multiple estimates from the posterior

#code here
plot(y ~ x, pch = 20, data = stand_Southdata)
for (i in 1:500) {
 abline(posterior_estimate$alpha[i], posterior_estimate$beta[i], col = "gray", lty = 1)
}
abline(mean(posterior_estimate$alpha), mean(posterior_estimate$beta), col = 6, lw = 2)

11. Change the priors and see how that affects your results

#code here
write("// Stan model for simple linear regression

data {
 int < lower = 1 > N; // Sample size
 vector[N] x; // Predictor
 vector[N] y; // Outcome
}

parameters {
 real alpha; // Intercept
 real beta; // Slope (regression coefficients)
 real < lower = 0 > sigma; // Error SD
}

model {
 alpha ~ normal(10, 0.1);
 beta ~ normal(1, 0.1);
 y ~ normal(alpha + x * beta , sigma);
}

generated quantities {}",

"stan_model2.stan")

stan_model2 <- "stan_model2.stan"

fit2 = stan(stan_model2, data = stand_Southdata, warmup = 500, iter = 1000, chains = 4, cores = 2, thin = 1)

posterior_estimate2 = rstan::extract(fit2)
str(posterior_estimate2)

## List of 4
##  $ alpha: num [1:2000(1d)] 10.2 10.1 10.2 10.2 10 ...
##   ..- attr(*, "dimnames")=List of 1
##   .. ..$ iterations: NULL
##  $ beta : num [1:2000(1d)] 0.0577 0.0603 0.0602 0.069 0.0719 ...
##   ..- attr(*, "dimnames")=List of 1
##   .. ..$ iterations: NULL
##  $ sigma: num [1:2000(1d)] 0.771 0.888 0.821 0.725 0.715 ...
##   ..- attr(*, "dimnames")=List of 1
##   .. ..$ iterations: NULL
##  $ lp__ : num [1:2000(1d)] -53 -54 -53.3 -53.1 -54.3 ...
##   ..- attr(*, "dimnames")=List of 1
##   .. ..$ iterations: NULL

plot(y ~ x, pch = 20, data = stand_Southdata)
abline(alpha, beta, col = 4, lty = 2, lw = 2)
abline(mean(posterior_estimate2$alpha), mean(posterior_estimate2$beta), col = 2, lw = 2)
abline(mean(posterior_estimate2$alpha), mean(posterior_estimate2$beta), col = 2, lw = 2)

12. What happened when you changed the priors? Does the model fit better or not?

The second model (adding more informative priors) really doesn’t fit well compared to the first one, since we restricted the possible values (we obviously did not enter valid values). Visually, we can see that by looking at the two ablines: they are going in two very different directions. The lines of the first model were one on the other.

13. Convergence diagnostics - create traceplots that show all 4 chains

We will plot using the first model, since there is no point in doing the convergence diagnostics of the second model with more priors: we alraedy know it doesn’t fit well, so it won’t converge well.

#code here

plot(posterior_estimate$alpha, type = "l")

plot(posterior_estimate$beta, type = "l")

plot(posterior_estimate$sigma, type = "l")

14. Plot parameter summaries for:

\(\alpha\), \(\beta\), \(\sigma\)

par(mfrow = c(1,3))

plot(density(posterior_estimate$alpha), main = "Alpha")
abline(v = alpha, col = 4, lty = 2)

plot(density(posterior_estimate$beta), main = "Beta")
abline(v = beta, col = 4, lty = 2)

plot(density(posterior_estimate$sigma), main = "Sigma")
abline(v=sigma, col = 4, lty = 2)

traceplot(fit)

stan_dens(fit)

stan_hist(fit)

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

plot(fit, show_density = FALSE, ci_level = 0.5, outer_level = 0.95, fill_color = "green")

## ci_level: 0.5 (50% intervals)

## outer_level: 0.95 (95% intervals)

15. What do your Stan model results indicate?

Posterior predictive checks

write("// Stan model for simple linear regression

data {
 int < lower = 1 > N; // Sample size
 vector[N] x; // Predictor
 vector[N] y; // Outcome
}

parameters {
 real alpha; // Intercept
 real beta; // Slope (regression coefficients)
 real < lower = 0 > sigma; // Error SD
}

model {
 y ~ normal(x * beta + alpha, sigma);
}

generated quantities {
 real y_rep[N];

 for (n in 1:N) {
 y_rep[n] = normal_rng(x[n] * beta + alpha, sigma);
 }

}",

"stan_model_pc.stan")

stan_model_final = "stan_model_pc.stan"

fit3 <- stan(stan_model_final, data = stand_Southdata, iter = 1000, chains = 4, cores = 2, thin = 1)
summary(fit3)

## $summary
##                  mean      se_mean          sd         2.5%          25%
## alpha     11.42380014 0.0047793839 0.133643005 11.156972507 11.332170638
## beta       0.01291032 0.0001958128 0.005814595  0.001433224  0.009287735
## sigma      0.39947399 0.0014901361 0.047677959  0.321210819  0.367468673
## y_rep[1]  11.42548127 0.0106590535 0.426010100 10.578168241 11.138296437
## y_rep[2]  11.44758895 0.0103305659 0.427127833 10.592463945 11.151115712
## y_rep[3]  11.46067187 0.0107067883 0.425168397 10.646715617 11.177226302
## y_rep[4]  11.47393002 0.0095794279 0.414310393 10.648341343 11.189953699
## y_rep[5]  11.49338718 0.0098489603 0.417134468 10.671784062 11.224319588
## y_rep[6]  11.50093099 0.0093681831 0.405509240 10.720842810 11.235724159
## y_rep[7]  11.52062259 0.0097764377 0.406538144 10.724179579 11.258699537
## y_rep[8]  11.51413912 0.0096373113 0.412004580 10.711493954 11.243306678
## y_rep[9]  11.54463307 0.0092065955 0.400199575 10.724282804 11.287776816
## y_rep[10] 11.54870447 0.0099794860 0.419169328 10.764663269 11.263693052
## y_rep[11] 11.56824017 0.0097699259 0.412253454 10.793707689 11.292078241
## y_rep[12] 11.58144788 0.0095272533 0.407369982 10.763519859 11.307980163
## y_rep[13] 11.58104063 0.0099375410 0.417827637 10.740167835 11.305579104
## y_rep[14] 11.59912632 0.0091467118 0.403829863 10.851532660 11.328367856
## y_rep[15] 11.60433987 0.0102077732 0.407340698 10.773848482 11.335057871
## y_rep[16] 11.64967134 0.0090120792 0.404383344 10.837106614 11.383912736
## y_rep[17] 11.65066262 0.0091304505 0.406981746 10.852155504 11.383794919
## y_rep[18] 11.64814981 0.0092252240 0.410035992 10.784757081 11.394536409
## y_rep[19] 11.67272273 0.0092841847 0.410533833 10.881995216 11.406458300
## y_rep[20] 11.65498129 0.0090118721 0.406164834 10.807868882 11.396986308
## y_rep[21] 11.69448318 0.0097942428 0.396008989 10.938457200 11.432751336
## y_rep[22] 11.69920356 0.0094715729 0.404539614 10.919981224 11.418989243
## y_rep[23] 11.72825085 0.0096668808 0.420000221 10.857208780 11.453109274
## y_rep[24] 11.72132950 0.0094474739 0.411985054 10.909138358 11.443805558
## y_rep[25] 11.74120128 0.0094793282 0.409753763 10.958521420 11.464699848
## y_rep[26] 11.76290511 0.0089921755 0.403202050 10.961256682 11.497904847
## y_rep[27] 11.77209189 0.0091599727 0.421809094 10.954037143 11.494947116
## y_rep[28] 11.78209764 0.0092950214 0.409192625 10.981713091 11.505276244
## y_rep[29] 11.80712305 0.0089997556 0.411715597 10.973662662 11.536328558
## y_rep[30] 11.81521950 0.0094581537 0.419242808 10.993678746 11.537613807
## y_rep[31] 11.82603717 0.0094972540 0.416432481 10.996530537 11.553456135
## y_rep[32] 11.83114762 0.0092762695 0.401777885 11.052153949 11.555993259
## y_rep[33] 11.84862333 0.0095511280 0.413624894 11.019614170 11.572642140
## y_rep[34] 11.88214673 0.0088031440 0.410106202 11.054319550 11.618711766
## y_rep[35] 11.87422527 0.0091393214 0.418217204 11.031393324 11.594673975
## y_rep[36] 11.88983964 0.0096626655 0.423403794 11.038501622 11.614605737
## y_rep[37] 11.90925868 0.0111842610 0.428473854 11.095911294 11.613042171
## y_rep[38] 11.89065074 0.0092420242 0.401077865 11.137941498 11.610552408
## y_rep[39] 11.92134212 0.0093510676 0.418128149 11.104602388 11.641514829
## lp__      16.25910289 0.0469776899 1.298328607 12.927337458 15.611341570
##                   50%        75%       97.5%     n_eff      Rhat
## alpha     11.42782839 11.5140360 11.69139987  781.8953 1.0027876
## beta       0.01284939  0.0165811  0.02486731  881.7731 1.0010376
## sigma      0.39378504  0.4281063  0.50649702 1023.7254 1.0034139
## y_rep[1]  11.42535481 11.7075934 12.28222446 1597.3589 0.9986915
## y_rep[2]  11.45607711 11.7231405 12.30400899 1709.4937 1.0007319
## y_rep[3]  11.46192271 11.7248397 12.30606658 1576.8977 1.0015730
## y_rep[4]  11.47612992 11.7500661 12.28581232 1870.5637 0.9991978
## y_rep[5]  11.49425002 11.7742864 12.33544893 1793.7891 0.9989841
## y_rep[6]  11.49551879 11.7682178 12.30610460 1873.6599 1.0001313
## y_rep[7]  11.52525964 11.7814096 12.29864317 1729.1845 0.9987473
## y_rep[8]  11.51297565 11.7814730 12.31091734 1827.6470 1.0012784
## y_rep[9]  11.55117203 11.8051343 12.29280669 1889.5358 1.0006310
## y_rep[10] 11.55174377 11.8109675 12.42584210 1764.2602 1.0004331
## y_rep[11] 11.57127121 11.8448790 12.38969593 1780.5168 0.9999775
## y_rep[12] 11.58027908 11.8425990 12.39735515 1828.2797 1.0001429
## y_rep[13] 11.58374299 11.8609033 12.41241344 1767.8135 1.0000398
## y_rep[14] 11.58863401 11.8589004 12.41999558 1949.2470 0.9987398
## y_rep[15] 11.60758986 11.8837755 12.37538142 1592.4052 0.9998865
## y_rep[16] 11.65427641 11.9149965 12.43070339 2013.4299 0.9988355
## y_rep[17] 11.64096563 11.9218122 12.43801027 1986.8517 1.0009318
## y_rep[18] 11.65542870 11.9259035 12.39952584 1975.5595 0.9982805
## y_rep[19] 11.67202504 11.9411062 12.47745441 1955.2863 1.0011827
## y_rep[20] 11.64442926 11.9387294 12.41935246 2031.3025 0.9998859
## y_rep[21] 11.68783966 11.9588145 12.46119201 1634.8140 1.0019105
## y_rep[22] 11.69513446 11.9871372 12.48795855 1824.2229 0.9988972
## y_rep[23] 11.73132396 12.0014217 12.56985983 1887.6710 0.9988175
## y_rep[24] 11.72207080 11.9900462 12.53435599 1901.6541 0.9993741
## y_rep[25] 11.73506166 12.0146117 12.55961125 1868.4900 0.9987127
## y_rep[26] 11.76831949 12.0375140 12.53176945 2010.5548 0.9988489
## y_rep[27] 11.77413404 12.0498501 12.59158604 2120.5256 0.9986560
## y_rep[28] 11.78143314 12.0506354 12.59242046 1938.0046 0.9989877
## y_rep[29] 11.80482844 12.0897174 12.62502470 2092.8264 1.0002335
## y_rep[30] 11.82994612 12.0949385 12.63432385 1964.8007 1.0004555
## y_rep[31] 11.82618034 12.0864641 12.67700559 1922.6184 1.0004277
## y_rep[32] 11.83487774 12.0942419 12.63238826 1875.9676 0.9990785
## y_rep[33] 11.84768251 12.1288545 12.66163442 1875.4436 1.0019736
## y_rep[34] 11.88083806 12.1572125 12.66205771 2170.2865 1.0025611
## y_rep[35] 11.87749166 12.1566548 12.70627136 2093.9962 1.0014053
## y_rep[36] 11.88466111 12.1718078 12.72079931 1920.0635 1.0012928
## y_rep[37] 11.91328009 12.1917054 12.77874695 1467.6891 0.9992402
## y_rep[38] 11.88761875 12.1568390 12.69000803 1883.3160 1.0005091
## y_rep[39] 11.92159081 12.1948116 12.73257556 1999.3850 0.9987932
## lp__      16.60654857 17.2322924 17.73378774  763.8111 1.0012418
## 
## $c_summary
## , , chains = chain:1
## 
##            stats
## parameter          mean          sd         2.5%          25%        50%
##   alpha     11.42365464 0.138138790 1.114921e+01 11.326691232 11.4271193
##   beta       0.01292728 0.006114776 8.022528e-04  0.009119907  0.0130393
##   sigma      0.39265482 0.045491327 3.211121e-01  0.360426126  0.3869260
##   y_rep[1]  11.43190117 0.426762067 1.055005e+01 11.138034899 11.4312635
##   y_rep[2]  11.43025045 0.411808954 1.065226e+01 11.165401044 11.4068883
##   y_rep[3]  11.45323527 0.409150606 1.066369e+01 11.200046313 11.4366230
##   y_rep[4]  11.46046639 0.407305668 1.065916e+01 11.186090588 11.4753418
##   y_rep[5]  11.47770284 0.419547823 1.062719e+01 11.197024485 11.4562681
##   y_rep[6]  11.48240854 0.385582556 1.076156e+01 11.198526055 11.4619571
##   y_rep[7]  11.50435843 0.406535712 1.072200e+01 11.239087444 11.4825560
##   y_rep[8]  11.49376795 0.382373067 1.073210e+01 11.236935348 11.5070641
##   y_rep[9]  11.52707772 0.383862481 1.074927e+01 11.287980484 11.5460239
##   y_rep[10] 11.54879804 0.399257671 1.077379e+01 11.261271694 11.5674749
##   y_rep[11] 11.57705322 0.410888355 1.080317e+01 11.310282935 11.5672374
##   y_rep[12] 11.56087502 0.407174255 1.071331e+01 11.285785896 11.5516759
##   y_rep[13] 11.55785151 0.408371574 1.076332e+01 11.275303922 11.5526006
##   y_rep[14] 11.59183128 0.387050263 1.084694e+01 11.332018525 11.5796636
##   y_rep[15] 11.61063253 0.407055842 1.079106e+01 11.345159196 11.6091542
##   y_rep[16] 11.65373644 0.386124661 1.085931e+01 11.401674308 11.6722930
##   y_rep[17] 11.64897486 0.401904309 1.085150e+01 11.370908505 11.6540532
##   y_rep[18] 11.65018996 0.382983851 1.087463e+01 11.410756231 11.6495456
##   y_rep[19] 11.64624735 0.389026910 1.090224e+01 11.367893275 11.6490733
##   y_rep[20] 11.65872747 0.397484680 1.078760e+01 11.398636202 11.6310428
##   y_rep[21] 11.71870050 0.383732321 1.098064e+01 11.475573015 11.6999555
##   y_rep[22] 11.70834565 0.389575233 1.099086e+01 11.449799853 11.7074505
##   y_rep[23] 11.72446701 0.411915993 1.084691e+01 11.443354848 11.7400603
##   y_rep[24] 11.69301285 0.401984998 1.091809e+01 11.414301970 11.6864559
##   y_rep[25] 11.73985434 0.403773288 1.098651e+01 11.459196930 11.7417760
##   y_rep[26] 11.74359909 0.415982836 1.092120e+01 11.468242766 11.7414913
##   y_rep[27] 11.77141515 0.419771286 1.097219e+01 11.490842360 11.7783671
##   y_rep[28] 11.78636177 0.404525386 1.104182e+01 11.509723793 11.7823512
##   y_rep[29] 11.79969500 0.411902488 1.098798e+01 11.529471866 11.7941989
##   y_rep[30] 11.79427813 0.401428536 1.098545e+01 11.531034360 11.8074008
##   y_rep[31] 11.81939421 0.413552177 1.102619e+01 11.553456135 11.8271837
##   y_rep[32] 11.84049521 0.404873089 1.105956e+01 11.565353343 11.8357206
##   y_rep[33] 11.85540966 0.405608333 1.101660e+01 11.599651567 11.8683319
##   y_rep[34] 11.89109581 0.396421554 1.104343e+01 11.655622137 11.8986968
##   y_rep[35] 11.89591884 0.412219363 1.107704e+01 11.630645631 11.9069712
##   y_rep[36] 11.88714984 0.419209436 1.103138e+01 11.614389457 11.8743311
##   y_rep[37] 11.91009547 0.437329022 1.108797e+01 11.624295440 11.9170918
##   y_rep[38] 11.88382215 0.368628722 1.110414e+01 11.658185818 11.8953072
##   y_rep[39] 11.94004029 0.411452923 1.116249e+01 11.668673981 11.9407035
##   lp__      16.23255404 1.309974265 1.278002e+01 15.620742508 16.5294035
##            stats
## parameter           75%       97.5%
##   alpha     11.51732016 11.69879237
##   beta       0.01690973  0.02472081
##   sigma      0.41785688  0.48877415
##   y_rep[1]  11.72597432 12.27418364
##   y_rep[2]  11.69945576 12.24883416
##   y_rep[3]  11.70351666 12.25866820
##   y_rep[4]  11.72279758 12.18596267
##   y_rep[5]  11.74541936 12.33261696
##   y_rep[6]  11.73034394 12.26467605
##   y_rep[7]  11.77462273 12.25601673
##   y_rep[8]  11.76172379 12.22489014
##   y_rep[9]  11.75284285 12.25133786
##   y_rep[10] 11.80758834 12.38604223
##   y_rep[11] 11.84025746 12.39272679
##   y_rep[12] 11.83233902 12.39171749
##   y_rep[13] 11.84083643 12.37222805
##   y_rep[14] 11.84923612 12.35727364
##   y_rep[15] 11.86727472 12.43847501
##   y_rep[16] 11.90533044 12.43048595
##   y_rep[17] 11.92025621 12.39933991
##   y_rep[18] 11.90765144 12.36574130
##   y_rep[19] 11.90676042 12.41751972
##   y_rep[20] 11.95242271 12.39242476
##   y_rep[21] 11.97309560 12.46897241
##   y_rep[22] 11.99663230 12.44186440
##   y_rep[23] 12.00173069 12.51982971
##   y_rep[24] 11.96708612 12.46187715
##   y_rep[25] 12.00781853 12.60476605
##   y_rep[26] 12.02966830 12.56756141
##   y_rep[27] 12.05347457 12.56861341
##   y_rep[28] 12.04307458 12.58088377
##   y_rep[29] 12.06411010 12.65772738
##   y_rep[30] 12.04843740 12.61605326
##   y_rep[31] 12.08105431 12.63773416
##   y_rep[32] 12.10016283 12.66979983
##   y_rep[33] 12.13149765 12.64584745
##   y_rep[34] 12.15772187 12.61093500
##   y_rep[35] 12.17949814 12.68697545
##   y_rep[36] 12.17103998 12.69308800
##   y_rep[37] 12.22670637 12.73872305
##   y_rep[38] 12.14894360 12.52136361
##   y_rep[39] 12.18718086 12.76896250
##   lp__      17.19113622 17.75481099
## 
## , , chains = chain:2
## 
##            stats
## parameter          mean          sd         2.5%          25%         50%
##   alpha     11.44260706 0.128217018 11.195727372 11.357583210 11.44461270
##   beta       0.01232166 0.005450687  0.001221125  0.009022061  0.01234454
##   sigma      0.40016689 0.046298333  0.319260064  0.369970287  0.39510927
##   y_rep[1]  11.41558432 0.415864709 10.621359093 11.142836121 11.43704118
##   y_rep[2]  11.47899226 0.413965624 10.657372122 11.196687024 11.48518130
##   y_rep[3]  11.48059794 0.415152331 10.670207795 11.213553935 11.50509665
##   y_rep[4]  11.48642421 0.426150190 10.593060062 11.224037245 11.48318533
##   y_rep[5]  11.50558542 0.421123590 10.676639621 11.240759817 11.50695015
##   y_rep[6]  11.53490214 0.404565913 10.741287642 11.255977980 11.53554344
##   y_rep[7]  11.53700778 0.402889163 10.751816056 11.286543464 11.53397936
##   y_rep[8]  11.55079016 0.408513034 10.760723932 11.298614687 11.54475420
##   y_rep[9]  11.55466579 0.404612859 10.764350567 11.301988780 11.57024715
##   y_rep[10] 11.57953270 0.433912883 10.751886348 11.282968594 11.57882827
##   y_rep[11] 11.55076708 0.389605480 10.846670332 11.279577654 11.56150891
##   y_rep[12] 11.60595730 0.410878350 10.811275093 11.337214780 11.60001046
##   y_rep[13] 11.59847406 0.408983552 10.782444337 11.332173753 11.59468918
##   y_rep[14] 11.60111243 0.416684390 10.830512966 11.330574111 11.58394532
##   y_rep[15] 11.62168642 0.404829776 10.767243879 11.371068651 11.62960001
##   y_rep[16] 11.66200159 0.398701727 10.906910700 11.405601736 11.64720973
##   y_rep[17] 11.64098278 0.415936487 10.852710325 11.372280266 11.62727831
##   y_rep[18] 11.65615003 0.425269608 10.767493643 11.398336152 11.65298517
##   y_rep[19] 11.67634520 0.430105228 10.855350643 11.416404833 11.66466621
##   y_rep[20] 11.68077171 0.387444403 10.896971541 11.436338065 11.67001333
##   y_rep[21] 11.73194781 0.388794302 10.982259332 11.481620866 11.71587415
##   y_rep[22] 11.69680923 0.421004426 10.902548406 11.413276508 11.69774413
##   y_rep[23] 11.71257958 0.416369500 10.870935744 11.442126065 11.72503759
##   y_rep[24] 11.73318334 0.427275082 10.924174819 11.460895621 11.73121593
##   y_rep[25] 11.74320838 0.408078019 11.001409759 11.455248141 11.73455570
##   y_rep[26] 11.76681604 0.413625970 10.916773429 11.490427677 11.78854206
##   y_rep[27] 11.77702295 0.423065987 10.923064082 11.492638522 11.79801338
##   y_rep[28] 11.76620725 0.398410865 10.988900649 11.495082845 11.76385750
##   y_rep[29] 11.80542213 0.407848831 10.953506335 11.534262543 11.79771454
##   y_rep[30] 11.82164946 0.423547848 11.031511070 11.535137956 11.81666557
##   y_rep[31] 11.84535167 0.406839496 11.019828199 11.567016277 11.86153102
##   y_rep[32] 11.81409188 0.391060488 11.071869476 11.550939919 11.81053912
##   y_rep[33] 11.82079956 0.402268645 10.958344060 11.547704402 11.83097723
##   y_rep[34] 11.85669323 0.426365386 11.042404129 11.564770180 11.83223107
##   y_rep[35] 11.86787201 0.420566938 11.032523487 11.591024184 11.86663877
##   y_rep[36] 11.86319532 0.444516352 10.950123272 11.589953637 11.85628337
##   y_rep[37] 11.92961155 0.420159964 11.119623007 11.631187416 11.94569486
##   y_rep[38] 11.89123563 0.410228841 11.122176386 11.617982280 11.88759328
##   y_rep[39] 11.91494117 0.416508650 11.151835167 11.611662629 11.90980596
##   lp__      16.38388788 1.256923410 13.390398368 15.706828429 16.68287789
##            stats
## parameter           75%       97.5%
##   alpha     11.51788574 11.70690923
##   beta       0.01561764  0.02329138
##   sigma      0.42509475  0.50510045
##   y_rep[1]  11.67350117 12.22813410
##   y_rep[2]  11.73945901 12.32879402
##   y_rep[3]  11.72651975 12.29574965
##   y_rep[4]  11.79331437 12.28634533
##   y_rep[5]  11.77381499 12.35944733
##   y_rep[6]  11.80079098 12.29968377
##   y_rep[7]  11.78789877 12.29541157
##   y_rep[8]  11.80635316 12.34228597
##   y_rep[9]  11.80822560 12.29818700
##   y_rep[10] 11.83951916 12.43919859
##   y_rep[11] 11.82717357 12.26702904
##   y_rep[12] 11.85830462 12.39398313
##   y_rep[13] 11.88064568 12.41281106
##   y_rep[14] 11.85442093 12.54093576
##   y_rep[15] 11.89462460 12.36627186
##   y_rep[16] 11.93939517 12.39162242
##   y_rep[17] 11.91042978 12.48484826
##   y_rep[18] 11.94177956 12.46108239
##   y_rep[19] 11.96090092 12.50629399
##   y_rep[20] 11.94520374 12.42108215
##   y_rep[21] 11.98638398 12.50638501
##   y_rep[22] 11.99367499 12.48752059
##   y_rep[23] 11.97731496 12.58926993
##   y_rep[24] 12.02829992 12.57366234
##   y_rep[25] 12.01036779 12.55752954
##   y_rep[26] 12.04085971 12.51478933
##   y_rep[27] 12.05400084 12.59561654
##   y_rep[28] 12.05349547 12.51422553
##   y_rep[29] 12.07844337 12.63288567
##   y_rep[30] 12.12206025 12.62443994
##   y_rep[31] 12.10082807 12.70463197
##   y_rep[32] 12.07487173 12.56968661
##   y_rep[33] 12.08536718 12.57695964
##   y_rep[34] 12.12892661 12.72489497
##   y_rep[35] 12.15296329 12.68634919
##   y_rep[36] 12.15060407 12.73136166
##   y_rep[37] 12.21832769 12.83852630
##   y_rep[38] 12.15528090 12.72263300
##   y_rep[39] 12.18747881 12.72494893
##   lp__      17.38019397 17.75058099
## 
## , , chains = chain:3
## 
##            stats
## parameter          mean          sd         2.5%          25%         50%
##   alpha     11.41645562 0.136304594 11.137389143 11.330590830 11.42466132
##   beta       0.01301724 0.005943822  0.001567519  0.009287735  0.01295912
##   sigma      0.40124617 0.050241702  0.325039026  0.364460739  0.39535225
##   y_rep[1]  11.44191671 0.418760381 10.619216169 11.151278786 11.40337942
##   y_rep[2]  11.45285304 0.443031438 10.529912253 11.144264926 11.47598470
##   y_rep[3]  11.48444891 0.426143329 10.689376261 11.211430093 11.48705353
##   y_rep[4]  11.47008613 0.407026168 10.697520133 11.179068228 11.46391523
##   y_rep[5]  11.48390080 0.394352795 10.696459061 11.228240694 11.49101398
##   y_rep[6]  11.49583290 0.410306667 10.678408899 11.241825039 11.49239447
##   y_rep[7]  11.52438578 0.389198614 10.776629532 11.268290720 11.53031224
##   y_rep[8]  11.49349748 0.423660001 10.684604276 11.206718803 11.48582920
##   y_rep[9]  11.56869972 0.414913008 10.807324762 11.288016224 11.56064314
##   y_rep[10] 11.54882314 0.413294482 10.776854308 11.278154574 11.55203198
##   y_rep[11] 11.56052015 0.422692772 10.730310714 11.266273438 11.56051501
##   y_rep[12] 11.56835163 0.399231090 10.795296623 11.284928161 11.56373463
##   y_rep[13] 11.59615003 0.429023481 10.711195807 11.323755418 11.61085916
##   y_rep[14] 11.59723799 0.416469680 10.864623418 11.306928178 11.59389743
##   y_rep[15] 11.60431325 0.415165062 10.742404209 11.334365296 11.60356468
##   y_rep[16] 11.64104512 0.408091492 10.809537340 11.378006564 11.63903554
##   y_rep[17] 11.66279571 0.411196766 10.784331757 11.405003939 11.64358821
##   y_rep[18] 11.64005235 0.408105267 10.702548724 11.407071348 11.65499156
##   y_rep[19] 11.69661454 0.411001973 10.907650237 11.415965336 11.69643504
##   y_rep[20] 11.62661372 0.428753884 10.742688026 11.373535243 11.63102579
##   y_rep[21] 11.66272334 0.404734609 10.904313195 11.370694660 11.66727951
##   y_rep[22] 11.68257128 0.398312649 10.903002147 11.410108538 11.66689865
##   y_rep[23] 11.74468933 0.425961257 10.902740348 11.479679919 11.74107160
##   y_rep[24] 11.73002074 0.406793253 10.911272503 11.469910246 11.75829250
##   y_rep[25] 11.74190004 0.427480949 10.914789177 11.465852559 11.74509213
##   y_rep[26] 11.77063015 0.392416921 10.954510456 11.522077317 11.76269987
##   y_rep[27] 11.75723450 0.428765535 10.817053182 11.512724475 11.76062165
##   y_rep[28] 11.77541098 0.427730551 10.937026839 11.506022434 11.78782358
##   y_rep[29] 11.82189720 0.421061566 10.962399667 11.549489612 11.82221568
##   y_rep[30] 11.81370298 0.418367670 10.993030351 11.563083933 11.82225228
##   y_rep[31] 11.82101929 0.422785588 11.008654474 11.566415467 11.81597894
##   y_rep[32] 11.82567549 0.384120408 11.040805813 11.555550117 11.83696585
##   y_rep[33] 11.86730406 0.425742406 11.020208889 11.572234049 11.85805943
##   y_rep[34] 11.86351308 0.401587540 11.087872718 11.619663514 11.86580875
##   y_rep[35] 11.84185595 0.414287071 10.969440062 11.577916178 11.84274708
##   y_rep[36] 11.90740105 0.402061347 11.196485985 11.642337948 11.91689108
##   y_rep[37] 11.87915943 0.427059429 11.085541167 11.585459460 11.86625411
##   y_rep[38] 11.89686464 0.407201657 11.187908517 11.597052595 11.87881118
##   y_rep[39] 11.91864967 0.417744797 11.126535229 11.629458409 11.89729314
##   lp__      16.16442969 1.310720306 12.915085129 15.514085842 16.54072515
##            stats
## parameter           75%       97.5%
##   alpha     11.50493416 11.68908323
##   beta       0.01665992  0.02540946
##   sigma      0.43288611  0.52135424
##   y_rep[1]  11.72706233 12.29276713
##   y_rep[2]  11.74623655 12.32970016
##   y_rep[3]  11.76153352 12.32880520
##   y_rep[4]  11.72683571 12.30637460
##   y_rep[5]  11.75050967 12.27488098
##   y_rep[6]  11.72431288 12.34604739
##   y_rep[7]  11.78219087 12.29384140
##   y_rep[8]  11.77828527 12.31756395
##   y_rep[9]  11.83833852 12.41884590
##   y_rep[10] 11.79403079 12.34483477
##   y_rep[11] 11.84340812 12.39310753
##   y_rep[12] 11.83578617 12.41524614
##   y_rep[13] 11.86718041 12.40746691
##   y_rep[14] 11.88397939 12.44592430
##   y_rep[15] 11.89360564 12.35729176
##   y_rep[16] 11.88877055 12.43331650
##   y_rep[17] 11.94263379 12.43844852
##   y_rep[18] 11.90874660 12.36265797
##   y_rep[19] 11.97422576 12.45925512
##   y_rep[20] 11.94213860 12.40496388
##   y_rep[21] 11.94378756 12.42862802
##   y_rep[22] 11.96338862 12.50219064
##   y_rep[23] 12.03319098 12.57726651
##   y_rep[24] 11.99081491 12.52705638
##   y_rep[25] 12.05159040 12.52764174
##   y_rep[26] 12.03472476 12.54212751
##   y_rep[27] 12.03380257 12.55622415
##   y_rep[28] 12.07192974 12.60334369
##   y_rep[29] 12.11839394 12.62901143
##   y_rep[30] 12.08269683 12.61923164
##   y_rep[31] 12.08050537 12.66122528
##   y_rep[32] 12.07520111 12.56144546
##   y_rep[33] 12.16146645 12.70837877
##   y_rep[34] 12.13306853 12.61337972
##   y_rep[35] 12.11118307 12.61748434
##   y_rep[36] 12.17249654 12.73694478
##   y_rep[37] 12.15611829 12.80680451
##   y_rep[38] 12.16176801 12.78659204
##   y_rep[39] 12.19274731 12.74428290
##   lp__      17.13931770 17.71049149
## 
## , , chains = chain:4
## 
##            stats
## parameter          mean          sd         2.5%          25%         50%
##   alpha     11.41248325 0.130044255 11.160783568 11.321928230 11.41052694
##   beta       0.01337509 0.005695298  0.002240218  0.009820039  0.01314768
##   sigma      0.40382807 0.047958046  0.321551857  0.371930798  0.39720481
##   y_rep[1]  11.41252286 0.442743235 10.513928174 11.110337667 11.43107221
##   y_rep[2]  11.42826006 0.438094814 10.590146760 11.138082239 11.45530115
##   y_rep[3]  11.42440536 0.447736530 10.531967393 11.114096785 11.43091908
##   y_rep[4]  11.47874334 0.417248279 10.673056750 11.181130872 11.46723263
##   y_rep[5]  11.50635966 0.433032093 10.630684897 11.241631905 11.52031411
##   y_rep[6]  11.49058038 0.420008358 10.696713649 11.218353803 11.48235980
##   y_rep[7]  11.51673836 0.427149601 10.700593725 11.247027622 11.53469727
##   y_rep[8]  11.51850089 0.430376536 10.643658101 11.225963693 11.51439783
##   y_rep[9]  11.52808907 0.396386966 10.622549886 11.268696139 11.53891441
##   y_rep[10] 11.51766398 0.428311659 10.740694882 11.231404494 11.50863763
##   y_rep[11] 11.58462023 0.425241393 10.828293338 11.299305763 11.59192760
##   y_rep[12] 11.59060758 0.411727605 10.772793298 11.310713523 11.60090667
##   y_rep[13] 11.57168693 0.424397826 10.709660011 11.280932756 11.59349331
##   y_rep[14] 11.60632358 0.395349379 10.859389261 11.338459379 11.59207264
##   y_rep[15] 11.58072726 0.402318082 10.781386542 11.298313475 11.58966850
##   y_rep[16] 11.64190220 0.424486457 10.802920967 11.372207297 11.67103107
##   y_rep[17] 11.64989714 0.399592932 10.901159659 11.396307056 11.64191683
##   y_rep[18] 11.64620689 0.423450461 10.802092996 11.351211610 11.67597026
##   y_rep[19] 11.67168385 0.410635798 10.893919055 11.407112171 11.68213366
##   y_rep[20] 11.65381224 0.409201568 10.905206273 11.365884140 11.64097489
##   y_rep[21] 11.66456105 0.402631086 10.911678186 11.369443289 11.65006788
##   y_rep[22] 11.70908808 0.409221146 10.905105497 11.417446010 11.70314845
##   y_rep[23] 11.73126750 0.426191994 10.890838865 11.463780481 11.72380865
##   y_rep[24] 11.72910107 0.411375794 10.913155699 11.434914616 11.72985712
##   y_rep[25] 11.73984235 0.400366455 10.989519402 11.473416550 11.70912079
##   y_rep[26] 11.77057517 0.390689413 11.018828942 11.519694243 11.76594413
##   y_rep[27] 11.78269497 0.416377255 11.025344218 11.494053140 11.76476276
##   y_rep[28] 11.80041055 0.405935423 11.059524325 11.521021168 11.79717273
##   y_rep[29] 11.80147788 0.406757947 11.030891169 11.520091816 11.81355937
##   y_rep[30] 11.83124744 0.433364871 10.976428449 11.532613569 11.86543398
##   y_rep[31] 11.81838351 0.422981388 10.972451118 11.544783112 11.80581132
##   y_rep[32] 11.84432788 0.426251397 11.060180618 11.553343043 11.84871464
##   y_rep[33] 11.85098006 0.420236730 11.088080935 11.569505673 11.83260858
##   y_rep[34] 11.91728480 0.413805505 11.096952198 11.639168967 11.91460871
##   y_rep[35] 11.89125429 0.424719301 11.061645037 11.601128293 11.89447395
##   y_rep[36] 11.90161237 0.426625968 11.107791740 11.620568694 11.89674697
##   y_rep[37] 11.91816826 0.428823046 11.142404743 11.603165308 11.92238181
##   y_rep[38] 11.89068054 0.417548976 11.150352506 11.586804935 11.87933831
##   y_rep[39] 11.91173735 0.427318299 10.988474321 11.667825332 11.93530347
##   lp__      16.25553995 1.309039226 12.841552593 15.624486745 16.66079005
##            stats
## parameter           75%      97.5%
##   alpha     11.51202759 11.6609327
##   beta       0.01704651  0.0241412
##   sigma      0.43200806  0.5061299
##   y_rep[1]  11.72737578 12.2356092
##   y_rep[2]  11.69226027 12.2753689
##   y_rep[3]  11.71913370 12.2825136
##   y_rep[4]  11.75074269 12.3238936
##   y_rep[5]  11.80127311 12.3446584
##   y_rep[6]  11.78791064 12.3683652
##   y_rep[7]  11.77494291 12.3839240
##   y_rep[8]  11.77687464 12.3772974
##   y_rep[9]  11.78612154 12.2601463
##   y_rep[10] 11.80090998 12.4135689
##   y_rep[11] 11.85689681 12.4305568
##   y_rep[12] 11.85587199 12.3759118
##   y_rep[13] 11.83217635 12.4194040
##   y_rep[14] 11.87224135 12.3614052
##   y_rep[15] 11.87355413 12.3166611
##   y_rep[16] 11.90513402 12.4559327
##   y_rep[17] 11.91668122 12.5057875
##   y_rep[18] 11.94210856 12.4417985
##   y_rep[19] 11.92674562 12.4791759
##   y_rep[20] 11.92104525 12.4160441
##   y_rep[21] 11.93371725 12.4387476
##   y_rep[22] 12.00054265 12.4951257
##   y_rep[23] 12.00969061 12.5685249
##   y_rep[24] 11.99484475 12.5521382
##   y_rep[25] 11.98792747 12.5565289
##   y_rep[26] 12.04112015 12.5022176
##   y_rep[27] 12.05404996 12.6233459
##   y_rep[28] 12.04543398 12.6426196
##   y_rep[29] 12.07133700 12.5880925
##   y_rep[30] 12.12504183 12.6498116
##   y_rep[31] 12.06860547 12.6654060
##   y_rep[32] 12.13199549 12.6878225
##   y_rep[33] 12.14091015 12.6962992
##   y_rep[34] 12.20034471 12.6459263
##   y_rep[35] 12.15667431 12.7392641
##   y_rep[36] 12.18039549 12.7303719
##   y_rep[37] 12.19293892 12.7898599
##   y_rep[38] 12.16833018 12.6993757
##   y_rep[39] 12.20308550 12.6697794
##   lp__      17.21072142 17.7076790

y_rep = as.matrix(fit3, pars = "y_rep")
y = stand_Southdata$y

ppc_dens_overlay(y, y_rep[1:200, ])

ppc_stat(y=y,yrep=y_rep,stat="mean")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Pretty good job! the data fits!