A linear regression is underdetermined when there are fewer observations than parameters. In this case the likelihood function does not concentrate on a compact neighborhood but rather an underdetermined hyperplane of degenerate model configurations. The corresponding posterior density has a surprisingly interesting geometry, even with weakly-informative prior densities that ensure a well-defined fit. In this short note I walk through the nature of this geometry and why underdetermined regressions are so hard to fit.

Let’s consider a basic linear regression with \(M\) slopes, \(\boldsymbol{\beta}\), a single intercept, \(\alpha\), and an unknown measurement variation, \(\sigma\). The regression will be underdetermined if there are fewer than \(M + 1\) observations to constrain all of the slopes and the intercept. Let’s see what happens when we try to fit a linear regression with many fewer observations.

First we simulate data using \(M = 200\) slopes and only \(N = 100\) observations.

writeLines(readLines("stan_programs/generate_data.stan"))
transformed data {
  int<lower=0> M = 200;
  int<lower=0> N = 100;
  real alpha = 3;
  real sigma = 1;
}

generated quantities {
  matrix[N, M] X;
  real y[N];

  vector[M] beta;
  for (m in 1:M) beta[m] = normal_rng(0, 10);

  for (n in 1:N) {
    for (m in 1:M) X[n, m] = normal_rng(0, 1);
    y[n] = normal_rng(X[n, 1:M] * beta + alpha, sigma);
  }
}
c_light <- c("#DCBCBC")
c_light_highlight <- c("#C79999")
c_mid <- c("#B97C7C")
c_mid_highlight <- c("#A25050")
c_dark <- c("#8F2727")
c_dark_highlight <- c("#7C0000")

c_light_trans <- c("#DCBCBC80")
c_dark_trans <- c("#8F272780")
c_green_trans <- c("#00FF0080")

library(rstan)
Warning: package 'rstan' was built under R version 3.5.2
Loading required package: StanHeaders
Warning: package 'StanHeaders' was built under R version 3.5.2
Loading required package: ggplot2
Warning: package 'ggplot2' was built under R version 3.5.2
rstan (Version 2.19.2, GitRev: 2e1f913d3ca3)
For execution on a local, multicore CPU with excess RAM we recommend calling
options(mc.cores = parallel::detectCores()).
To avoid recompilation of unchanged Stan programs, we recommend calling
rstan_options(auto_write = TRUE)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

util <- new.env()
source('stan_utility.R', local=util)

fit <- stan(file='stan_programs/generate_data.stan', iter=1,
            chains=1, seed=194838, algorithm="Fixed_param")

SAMPLING FOR MODEL 'generate_data' NOW (CHAIN 1).
Chain 1: Iteration: 1 / 1 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 0 seconds (Warm-up)
Chain 1:                0.002183 seconds (Sampling)
Chain 1:                0.002183 seconds (Total)
Chain 1: 
X <- extract(fit)$X[1,,]
y <- extract(fit)$y[1,]
N <- dim(X)[1]
M <- dim(X)[2]

stan_rdump(c("N", "M", "X", "y"), file="linear_regression.data.R")

We then proceed to fit the simulated data using weakly informative priors on all of the parameters to ensure a well-defined posterior distribution.

writeLines(readLines("stan_programs/linear_regression.stan"))
data {
  int<lower=1> N; // Number of observations
  int<lower=1> M; // Number of covariates
  matrix[N, M] X; // Design matrix
  vector[N] y;    // Observations
}

parameters {
  vector[M] beta;
  real alpha;
  real<lower=0> sigma;
}

model {
  // Weakly informative containment priors
  beta ~ normal(0, 10);
  alpha ~ normal(0, 2);
  sigma ~ normal(0, 2);

  // Observational model
  y ~ normal(X * beta + alpha, sigma);
}
input_data <- read_rdump("linear_regression.data.R")
fit <- stan(file='stan_programs/linear_regression.stan', data=input_data, 
            seed=4938483, refresh=0)
Warning: There were 2980 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
http://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded
Warning: There were 4 chains where the estimated Bayesian Fraction of Missing Information was low. See
http://mc-stan.org/misc/warnings.html#bfmi-low
Warning: Examine the pairs() plot to diagnose sampling problems
Warning: The largest R-hat is 1.64, indicating chains have not mixed.
Running the chains for more iterations may help. See
http://mc-stan.org/misc/warnings.html#r-hat
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
http://mc-stan.org/misc/warnings.html#bulk-ess
Warning: 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
http://mc-stan.org/misc/warnings.html#tail-ess
util$check_all_diagnostics(fit)
[1] "n_eff / iter for parameter beta[31] is 0.000894106945654838!"
[1] "n_eff / iter for parameter beta[135] is 0.000840438620335838!"
[1] "n_eff / iter for parameter beta[163] is 0.000997074034378899!"
[1] "n_eff / iter for parameter lp__ is 0.000663475092983791!"
[1] "  n_eff / iter below 0.001 indicates that the effective sample size has likely been overestimated"
[1] "Rhat for parameter beta[1] is 1.2087035666718!"
[1] "Rhat for parameter beta[3] is 1.14368775214423!"
[1] "Rhat for parameter beta[5] is 1.15030561443003!"
[1] "Rhat for parameter beta[6] is 1.14170383362239!"
[1] "Rhat for parameter beta[7] is 1.12484077022282!"
[1] "Rhat for parameter beta[8] is 1.14362022749256!"
[1] "Rhat for parameter beta[10] is 1.28566132218153!"
[1] "Rhat for parameter beta[11] is 1.1659850880562!"
[1] "Rhat for parameter beta[12] is 1.14585251957181!"
[1] "Rhat for parameter beta[15] is 1.11001386212505!"
[1] "Rhat for parameter beta[17] is 1.14341599655567!"
[1] "Rhat for parameter beta[20] is 1.17554214057238!"
[1] "Rhat for parameter beta[21] is 1.18463554992392!"
[1] "Rhat for parameter beta[23] is 1.12818017211599!"
[1] "Rhat for parameter beta[24] is 1.14777257093343!"
[1] "Rhat for parameter beta[25] is 1.13609282723697!"
[1] "Rhat for parameter beta[26] is 1.47028166339406!"
[1] "Rhat for parameter beta[30] is 1.22187805961954!"
[1] "Rhat for parameter beta[31] is 1.53764034253115!"
[1] "Rhat for parameter beta[32] is 1.20010748912073!"
[1] "Rhat for parameter beta[37] is 1.1747597536018!"
[1] "Rhat for parameter beta[39] is 1.15376333526228!"
[1] "Rhat for parameter beta[40] is 1.24704395413346!"
[1] "Rhat for parameter beta[43] is 1.20019149609663!"
[1] "Rhat for parameter beta[45] is 1.16817780559862!"
[1] "Rhat for parameter beta[49] is 1.22267090887719!"
[1] "Rhat for parameter beta[50] is 1.19580172793858!"
[1] "Rhat for parameter beta[53] is 1.14964728282493!"
[1] "Rhat for parameter beta[56] is 1.21941880631021!"
[1] "Rhat for parameter beta[61] is 1.21480649593076!"
[1] "Rhat for parameter beta[63] is 1.11408154004105!"
[1] "Rhat for parameter beta[64] is 1.36418462214846!"
[1] "Rhat for parameter beta[65] is 1.16265479740955!"
[1] "Rhat for parameter beta[66] is 1.45423594961522!"
[1] "Rhat for parameter beta[70] is 1.10518147951655!"
[1] "Rhat for parameter beta[73] is 1.12728835352287!"
[1] "Rhat for parameter beta[77] is 1.17301356009538!"
[1] "Rhat for parameter beta[79] is 1.1080801726339!"
[1] "Rhat for parameter beta[80] is 1.24725288449931!"
[1] "Rhat for parameter beta[82] is 1.12063642876362!"
[1] "Rhat for parameter beta[83] is 1.21282442807188!"
[1] "Rhat for parameter beta[84] is 1.11247567065969!"
[1] "Rhat for parameter beta[86] is 1.16639834054377!"
[1] "Rhat for parameter beta[87] is 1.11644718102863!"
[1] "Rhat for parameter beta[89] is 1.12797207212025!"
[1] "Rhat for parameter beta[90] is 1.18851418706485!"
[1] "Rhat for parameter beta[91] is 1.15849889154711!"
[1] "Rhat for parameter beta[92] is 1.41215036059405!"
[1] "Rhat for parameter beta[94] is 1.10548908065154!"
[1] "Rhat for parameter beta[96] is 1.40120055335768!"
[1] "Rhat for parameter beta[100] is 1.1050759244889!"
[1] "Rhat for parameter beta[101] is 1.1192010780607!"
[1] "Rhat for parameter beta[103] is 1.36760053617081!"
[1] "Rhat for parameter beta[106] is 1.13384709228827!"
[1] "Rhat for parameter beta[108] is 1.32223023787767!"
[1] "Rhat for parameter beta[110] is 1.28364329211358!"
[1] "Rhat for parameter beta[111] is 1.21916494782873!"
[1] "Rhat for parameter beta[113] is 1.12489631345038!"
[1] "Rhat for parameter beta[114] is 1.135617155915!"
[1] "Rhat for parameter beta[115] is 1.12350922479278!"
[1] "Rhat for parameter beta[116] is 1.15610992115841!"
[1] "Rhat for parameter beta[117] is 1.20479708534534!"
[1] "Rhat for parameter beta[123] is 1.19461961226087!"
[1] "Rhat for parameter beta[124] is 1.18432668600871!"
[1] "Rhat for parameter beta[127] is 1.21952603368409!"
[1] "Rhat for parameter beta[128] is 1.46970161291643!"
[1] "Rhat for parameter beta[130] is 1.12090069994594!"
[1] "Rhat for parameter beta[131] is 1.36023495867876!"
[1] "Rhat for parameter beta[132] is 1.32609175051137!"
[1] "Rhat for parameter beta[133] is 1.2205369811441!"
[1] "Rhat for parameter beta[134] is 1.41840062831778!"
[1] "Rhat for parameter beta[135] is 1.64979805588818!"
[1] "Rhat for parameter beta[138] is 1.18568840962408!"
[1] "Rhat for parameter beta[139] is 1.1412827168933!"
[1] "Rhat for parameter beta[140] is 1.17942892087985!"
[1] "Rhat for parameter beta[141] is 1.1185641657472!"
[1] "Rhat for parameter beta[143] is 1.1833633942917!"
[1] "Rhat for parameter beta[148] is 1.26464824429031!"
[1] "Rhat for parameter beta[152] is 1.16385605212812!"
[1] "Rhat for parameter beta[153] is 1.24079046932187!"
[1] "Rhat for parameter beta[154] is 1.18456681221102!"
[1] "Rhat for parameter beta[156] is 1.13000287868956!"
[1] "Rhat for parameter beta[157] is 1.14917799760115!"
[1] "Rhat for parameter beta[160] is 1.23496451632675!"
[1] "Rhat for parameter beta[163] is 1.44910606810777!"
[1] "Rhat for parameter beta[165] is 1.16643899896794!"
[1] "Rhat for parameter beta[166] is 1.23673899370296!"
[1] "Rhat for parameter beta[170] is 1.36163309162877!"
[1] "Rhat for parameter beta[172] is 1.10601807277243!"
[1] "Rhat for parameter beta[175] is 1.19220852600143!"
[1] "Rhat for parameter beta[177] is 1.26590570312247!"
[1] "Rhat for parameter beta[178] is 1.3463504753882!"
[1] "Rhat for parameter beta[180] is 1.16461028061245!"
[1] "Rhat for parameter beta[181] is 1.35907116238207!"
[1] "Rhat for parameter beta[182] is 1.22192280273783!"
[1] "Rhat for parameter beta[183] is 1.30914751042247!"
[1] "Rhat for parameter beta[184] is 1.18348772230506!"
[1] "Rhat for parameter beta[185] is 1.10409458801282!"
[1] "Rhat for parameter beta[186] is 1.32875951296858!"
[1] "Rhat for parameter beta[188] is 1.1853605345807!"
[1] "Rhat for parameter beta[190] is 1.26333333826128!"
[1] "Rhat for parameter beta[191] is 1.11796114375408!"
[1] "Rhat for parameter beta[192] is 1.13729647102845!"
[1] "Rhat for parameter beta[193] is 1.19000121203031!"
[1] "Rhat for parameter beta[194] is 1.12746632759695!"
[1] "Rhat for parameter beta[195] is 1.28246645252255!"
[1] "Rhat for parameter beta[198] is 1.1690307330544!"
[1] "Rhat for parameter beta[199] is 1.31418504807679!"
[1] "Rhat for parameter alpha is 1.24351520988763!"
[1] "Rhat for parameter sigma is 1.49429441152495!"
[1] "Rhat for parameter lp__ is 2.37694563782786!"
[1] "  Rhat above 1.1 indicates that the chains very likely have not mixed"
[1] "0 of 4000 iterations ended with a divergence (0%)"
[1] "2980 of 4000 iterations saturated the maximum tree depth of 10 (74.5%)"
[1] "  Run again with max_treedepth set to a larger value to avoid saturation"
[1] "Chain 1: E-FMI = 0.0281663385480791"
[1] "Chain 2: E-FMI = 0.0606133949004802"
[1] "Chain 3: E-FMI = 0.0688824339931401"
[1] "Chain 4: E-FMI = 0.150753706810018"
[1] "  E-FMI below 0.2 indicates you may need to reparameterize your model"

Whoa, that’s unpleasant. According to the split \(\hat{R}\) for each parameter the four chains we ran didn’t mix at all, and the E-FMI warnings indicate that the geometry encountered by each individual chain was not particularly pleasant. There are no divergences but the treedepth saturation warning indicates that the adapted step sizes were quite small, further suggesting an unpleasant posterior geometry. Before running with a larger maximum treedepth let’s build up some intuition about what to expect in our likelihood function, and hence our fit.

Typically as the measurement variation vanishes the likelihood function would concentrate on a single configuration of slopes and intercept, \(\tilde{\boldsymbol{\gamma}} = (\tilde{\boldsymbol{\beta}}, \tilde{\alpha})\).
Because we don’t have enough observations, however, the likelihood function in our underdetermined regression concentrates on a hyperplane of configurations. Ignoring the small contribution due to the priors, this hyperplane is defined by \[ \tilde{\boldsymbol{\gamma}} = \mathbf{R} + (\mathbf{I} - \mathbf{P}_{\mathbf{X}}) \cdot \mathbf{w} \] where \[ \mathbf{R} = \mathbf{X}^{+} \cdot \mathbf{y}, \] \[ \mathbf{P}_{\mathbf{X}} = \mathbf{X}^{+} \cdot (\mathbf{X} \cdot \mathbf{X}^{+})^{-1} \cdot \mathbf{X}, \] and \(\mathbf{w}\) is any vector. Here \(\mathbf{X}\) is the buffered design matrix that includes a column of ones to incorporate the intercept and \(^{+}\) indicates the Moore-Penrose pseudo-inverse. In practice the action of the pseudo-inverse on a vector, \(\mathbf{b}\), can be calculated by solving the linear system \[ (\mathbf{X}^{T} \cdot \mathbf{X}) \cdot (\mathbf{X}^{+} \cdot \mathbf{b}) \cdot \mathbf{X}^{T} \cdot \mathbf{b}. \]

Geometrically \(\mathbf{P}_{X}\) defines a projection operator that projects out the component of a vector orthogonal to the underdetermined hyperplane. At the same time \(\mathbf{I} - \mathbf{P}_{X}\) defines a projection operator that projects out the component of a vector along the underdetermined hyperplane. In other words, any point on the underdetermined hyperplane, \(\tilde{\boldsymbol{\gamma}}\), can be decomposed into the vector from the origin to the hyerplane, \(\mathbf{R}\), and a vector along the hyperplane,
\((\mathbf{I} - \mathbf{P}_{\mathbf{X}}) \cdot \mathbf{w}\).





Moreover, we can decompose any configuration \(\boldsymbol{\gamma}\) into the vector \(R\) and the distance from the configuration to the underdetermined hyperplane, \(\Delta \boldsymbol{\gamma}\). That distance can then be further decomposed into perpendicular and parallel components using the projection operator, \(\mathbf{P}_{\mathbf{X}}\).





As the measurement variation, \(\sigma\), increases away from zero the concentration of the likelihood function, and hence our posterior density, around the underdetermined hyperplane, weakens. Because we include \(\sigma\) as a parameter in our model, our fit will have to contend with this varying concentrations around the hyperplane.

Let’s try our fit again, only this time calculating the perpendicular and parallel components of the distance from each sampled configuration and the underdetermined hyperplane.

writeLines(readLines("stan_programs/linear_regression_w_proj.stan"))
data {
  int<lower=1> N; // Number of observations
  int<lower=1> M; // Number of covariates
  matrix[N, M] X; // Design matrix
  vector[N] y;    // Observations
}

transformed data {
  // Buffered design matrix to incorporate intercept
  matrix[N, M + 1] X_buff = append_col(X, rep_vector(1.0, N));
  matrix[M + 1, M + 1] Z = X_buff' * X_buff;
  
  // Vector from origin to underdetermined hyperplane
  vector[M + 1] R = Z \ (X_buff' * y);
  
  // Normalization of projection operator onto underdetermined hyperplane
  matrix[N, N] Pnorm = inverse(X_buff * (Z \ X_buff') );
}

parameters {
  vector[M] beta;
  real alpha;
  real<lower=0> sigma;
}

model {
  // Weakly informative containment priors
  beta ~ normal(0, 10);
  alpha ~ normal(0, 2);
  sigma ~ normal(0, 2);

  // Observational model
  // Equivalent to y ~ normal(X * beta + alpha, sigma);
  y ~ normal(X_buff * append_row(beta, alpha), sigma);
}

generated quantities {
  // Projection of parameters orthogonal to underdetermined hyperplane
  vector[M + 1] gamma_perp;

  // Projection of parameters along underdetermined hyperplane
  vector[M + 1] gamma_par;
  
  // Distance from parameters to underdetermined hyperplane
  real delta;
  
  {
    vector[N] z = Pnorm * X_buff * (append_row(beta, alpha) - R);
    gamma_perp = Z \ (X_buff' * z);
    gamma_par = (append_row(beta, alpha) - R) - gamma_perp;
    delta = sqrt(dot_self(gamma_perp));
  }
}

This time we also make sure to give our trajectories the room they need by expanding the maximum treedepth.

fit <- stan(file='stan_programs/linear_regression_w_proj.stan', data=input_data, 
            seed=4938483, control=list(max_treedepth=15), refresh=0)
Warning: There were 4 chains where the estimated Bayesian Fraction of Missing Information was low. See
http://mc-stan.org/misc/warnings.html#bfmi-low
Warning: Examine the pairs() plot to diagnose sampling problems
Warning: The largest R-hat is 1.05, indicating chains have not mixed.
Running the chains for more iterations may help. See
http://mc-stan.org/misc/warnings.html#r-hat
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
http://mc-stan.org/misc/warnings.html#bulk-ess
Warning: 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
http://mc-stan.org/misc/warnings.html#tail-ess
util$check_n_eff(fit)
[1] "n_eff / iter looks reasonable for all parameters"
util$check_rhat(fit)
[1] "Rhat looks reasonable for all parameters"
util$check_div(fit)
[1] "0 of 4000 iterations ended with a divergence (0%)"
util$check_treedepth(fit, 15)
[1] "0 of 4000 iterations saturated the maximum tree depth of 15 (0%)"
util$check_energy(fit)
[1] "Chain 1: E-FMI = 0.0757556173806369"
[1] "Chain 2: E-FMI = 0.0780362035619746"
[1] "Chain 3: E-FMI = 0.0630098579117811"
[1] "Chain 4: E-FMI = 0.0446730286729114"
[1] "  E-FMI below 0.2 indicates you may need to reparameterize your model"

There are no more effective sample size, split \(\hat{R}\), or treedepth warnings, but the energy fraction of missing information warnings persist. Let’s explore what’s happening around the underdetermined hyperplane.

As \(\sigma\) decreases our samples are drawn towards the hyperplane which manifests in a strong correlation between \(\sigma\) and the distance \(\delta\). The posterior doesn’t push all the way to \(\sigma = 0\) because the few observations that we do have are enough to contain the likelihood for the measurement variability above \(\sigma \approx 0.5\).

partition <- util$partition_div(fit)
div_params <- partition[[1]]
nondiv_params <- partition[[2]]

plot(log(nondiv_params$sigma), log(nondiv_params$delta),
     col=c_dark_trans, pch=16, cex=0.8, main="",
     xlab="log(sigma)", ylab="log(delta)", xlim=c(-3, 2))
points(log(div_params$sigma), log(div_params$delta),
       col=c_green_trans, pch=16, cex=0.8)

More troubling, the shrinkage towards the hyperplane pulls the perpendicular components of the model configuration towards zero in both directions which manifests as a funnel! The tapering of the funnel is limited by the posterior for \(\sigma\), but it is still challenging enough to require small integrator step sizes and long numerical trajectories. Had we fewer observations the funnel would collapse even further and make the fitting all the more difficult.

plot(nondiv_params$"gamma_perp[1]", log(nondiv_params$sigma),
     col=c_dark_trans, pch=16, cex=0.8, main="",
     xlab="gamma_perp[1]", ylab="log(sigma)")
points(div_params$"gamma_perp[1]", log(div_params$sigma),
       col=c_green_trans, pch=16, cex=0.8)

The parallel components also manifest a bit of a funnel, but instead of shrinking towards zero they just shrink towards the bulk of our weakly informative priors. Without these priors to contain the posterior distribution the parallel components would expand out towards infinity regardless of the value of \(\sigma\).

plot(nondiv_params$"gamma_par[1]", log(nondiv_params$sigma),
     col=c_dark_trans, pch=16, cex=0.8, main="",
     xlab="gamma_par[1]", ylab="log(sigma)")
points(div_params$"gamma_par[1]", log(div_params$sigma),
       col=c_green_trans, pch=16, cex=0.8)

Unfortunately there’s not much that we can do to ameliorate the geometry of the underdetermined regression posterior. Unlike hierarchical models, the funnel manifests directly in the observational model and hence isn’t immediately amenable to reparameterizations. A prior incorporating additional domain expertise about the minimal reasonable scale for \(\sigma\) within a specific measurement process could cut off more of the funnel. With only \(N = 100\) observations, however, the posterior already learns to suppress \(\sigma \lesssim 0.5\) which isn’t far from the true value of \(\sigma = 1\). Consequently there isn’t much room for additional information.

Ultimately the geometry of underdetermined regression posteriors is not particularly well-researched and optimal implementations are still very open areas of research!

1 Acknowledgements

I thank Frank Weber for helpful comments.

A very special thanks to everyone supporting me on Patreon: Aki Vehtari, Alan O’Donnell, Andre Zapico, Andrew Rouillard, Austin Rochford, Avraham Adler, Bo Schwartz Madsen, Bryan Yu, Cat Shark, Charles Naylor, Christopher Howlin, Colin Carroll, Daniel Simpson, David Pascall, David Roher, David Stanard, Ed Cashin, Eddie Landesberg, Elad Berkman, Eric Jonas, Ethan Goan, Finn Lindgren, Granville Matheson, Hernan Bruno, J Michael Burgess, Joel Kronander, Jonas Beltoft Gehrlein, Joshua Mayer, Justin Bois, Lars Barquist, Luiz Carvalho, Marek Kwiatkowski, Matthew Kay, Maurits van der Meer, Maxim Kesin, Michael Dillon, Michael Redman, Noah Silverman, Ole Rogeberg, Oscar Olvera, Paul Oreto, Peter Heinrich, Putra Manggala, Ravin Kumar, Riccardo Fusaroli, Richard Torkar, Robert Frost, Robin Taylor, Sam Petulla, Sam Zorowitz, Seo-young Kim, Seth Axen, Sharan Banagiri, Simon Duane, Stephen Oates, Stijn, Vladislavs Dovgalecs, and yolha.

License

A repository containing the material used in this case study is available on GitHub.

The code in this case study is copyrighted by Michael Betancourt and licensed under the new BSD (3-clause) license:

https://opensource.org/licenses/BSD-3-Clause

The text and figures in this case study are copyrighted by Michael Betancourt and licensed under the CC BY-NC 4.0 license:

https://creativecommons.org/licenses/by-nc/4.0/

Original Computing Environment

writeLines(readLines(file.path(Sys.getenv("HOME"), ".R/Makevars")))
CC=clang

CXXFLAGS=-O3 -mtune=native -march=native -Wno-unused-variable -Wno-unused-function -Wno-macro-redefined -Wno-unneeded-internal-declaration
CXX=clang++ -arch x86_64 -ftemplate-depth-256

CXX14FLAGS=-O3 -mtune=native -march=native -Wno-unused-variable -Wno-unused-function -Wno-macro-redefined -Wno-unneeded-internal-declaration -Wno-unknown-pragmas
CXX14=clang++ -arch x86_64 -ftemplate-depth-256
sessionInfo()
R version 3.5.1 (2018-07-02)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS  10.14.6

Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] rstan_2.19.2          ggplot2_3.1.1         StanHeaders_2.18.1-10

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.3         pillar_1.4.2       compiler_3.5.1     plyr_1.8.4        
 [5] prettyunits_1.0.2  base64enc_0.1-3    tools_3.5.1        digest_0.6.18     
 [9] pkgbuild_1.0.2     evaluate_0.14      tibble_2.1.3       gtable_0.2.0      
[13] pkgconfig_2.0.3    rlang_0.4.2        cli_1.0.1          parallel_3.5.1    
[17] yaml_2.2.0         xfun_0.11          loo_2.0.0          gridExtra_2.3     
[21] withr_2.1.2        dplyr_0.8.3        stringr_1.3.1      knitr_1.26        
[25] stats4_3.5.1       grid_3.5.1         tidyselect_0.2.5   glue_1.3.0        
[29] inline_0.3.15      R6_2.3.0           processx_3.2.0     rmarkdown_1.18    
[33] purrr_0.3.3        callr_3.0.0        magrittr_1.5       codetools_0.2-15  
[37] matrixStats_0.54.0 ps_1.2.0           scales_1.0.0       htmltools_0.4.0   
[41] assertthat_0.2.0   colorspace_1.3-2   stringi_1.2.4      lazyeval_0.2.1    
[45] munsell_0.5.0      crayon_1.3.4