manual longer term

[bob-carpenter]

This is the issue for things that would be nice to get into the manual in the future. Rather than closing this issue, we'll just tick things off as part of the regular round of "next manual" issues.

[bob-carpenter]

Originally from Jeffrey Arnold here: stan-dev/stan#1081 (comment)

Regarding the Kalman filter examples, I have examples and fully implemented Kalman filtering (of several flavors: with / without missing values, as batch / sequentially) with backward sampling in the generated quantities block here: https://github.com/jrnold/ssmodels-in-stan. I could probably write up that section.

[bob-carpenter]

• reparameterize Gaussian process example to use length scale with a proper prior (e.g., half-t(4)).

[bob-carpenter]

• add horseshoe prior example using Allen's implementation:

https://ariddell.org/horseshoe-prior-with-stan.html

• explore what happens with L1 priors in Stan's MLE to see if we get truly 0 values or close; I'm not sure what L-BFGS will do

[bob-carpenter]

• add section or new chapter about arithmetic precision
• include following example I sent to stan-users:

For example, this model is problematic arithmetically:

fit <- stan(model_code = "parameters { real y; } model { y ~ normal(1,1e-20); }")

and produces

          mean  se_mean       sd      2.5%       25% 50% 75% 97.5% n_eff         Rhat
y            1      0.0        0         1         1   1   1     1  4000          Inf
lp__ -15407440 422002.9 26689805 -61629758 -15407440   0   0     0  4000 1.500106e+14

There's really nothing we can do about such models numerically. What you have to do from a user perspective is use a parameterizations z = y - 1, which can be handled stably because you can represent numbers close to zero well.

parameters { real z; } model { z ~ normal(0,1e-20); }

The problem is that there's still no way to represent y = z + 1 as a real value using double-precision floating-point arithmetic, but you can get close to log(y) by using

  log_y <- log1p(z);

[bob-carpenter]

Provide examples of log-mix for the manual mixture section.

[bob-carpenter]

From Jeffrey Arnold:

Regarding the Kalman filter examples, I have examples and fully implemented Kalman filtering (of several flavors: with / without missing values, as batch / sequentially) with backward sampling in the generated quantities block here: https://github.com/jrnold/ssmodels-in-stan. I could probably write up that section.

[bob-carpenter]

From Andrew on stan-dev, an example of a model for the problematic posteriors section:

Next I decided, just for fun, to fit a nonlinear model: y = a_exp(-b_x) + c + error, with N(0,sigma) errors. I simulated some data (x taking on the values between 0 and 23, setting a=.5, b=1/10, c=1, sigma=.2, and I fit it constraining all the parameters to be positive. (In retrospect I’m thinking I should’ve just constrained b and sigma to be positive; constraining a and c made sense in the context of my hypothetical example, but in general it’s poor practice to put in hard constraints unless absolutely necessary.)

I simulated data and fit the model and all was fine. But then I repeated and the fit blew up, I was getting some estimates of b that were either really close to 0 (1e-10, etc) or really high (1e10, etc), I can’t remember which. The point is that there can be instability at the extremes, and I realized the model need a bit of hand-holding. As a kluge I constrained b to be between 1/50 and 1 (which made sense in the context of the example). A hard constraint was ugly but was easy enough to understand. I didn’t want to confuse the class by introducing a normal prior here.

So this second example was a bit of a mess. Which suggests to me that in the manual/book we need to throw in an example to explicitly address instability and constraints. Since, ultimately, this is an inherent difficulty with modeling, it’s not anything to do with Stan in particular.

[bob-carpenter]

• add discussion of K vs. (K-1) coefficient vector parameterizations of multi-logit

Here's what the model should look like, but there's a bug in Stan 2.6 preventing the compilation:

data {
  int<lower=2> K;             // num categories
  int<lower=1> D;             // num predictors (without intercept)
  int<lower=0> N;             // num observations
  vector[D] x[N];              // predictors
  int<lower=1,upper=K> y[N];  // observations
}
transformed data {
  vector[1] zero;
  zero[1] <- 0;
}
parameters {
  vector[K-1] alpha;          // intercepts
  matrix[K-1,D] beta;         // slopes
}
model {
  // priors
  alpha ~ normal(0,20);
  to_vector(beta) ~ normal(0,10);
  // likelihood
  for (n in 1:N) {
    y[n] ~ categorical_logit(append_row(zero, alpha + beta * x[n]));
  }  
}

Here's a program that will work, but isn't optimally efficient:

data {
 int<lower=2> K;             // num categories
 int<lower=1> D;             // num predictors (without intercept)
 int<lower=0> N;             // num observations
 vector[D] x[N];              // predictors
 int<lower=1,upper=K> y[N];  // observations
}
transformed data {
 vector[1] zero;
 zero[1] <- 0;
}
parameters {
 vector[K-1] alpha;          // intercepts
 matrix[K-1,D] beta;         // slopes
}
model {
 // temps
 vector[K] lambda;
 lambda[1] <- 0;

 // priors
 alpha ~ normal(0,20);
 // likelihood
 for (n in 1:N) {
   for (k in 2:K)
     lambda[k] <- alpha[k-1] + beta[k-1] * x[n];
   y[n] ~ categorical_logit(lambda);
 }
}

[bob-carpenter]

• Add a chapter on using Stan for pure Monte Carlo integration

Computing pi

/**
 * Compute the mathematical constant pi by Monte Carlo integration.
 *
 * First, sample (x,y) uniformly from a square with vertices (-1,1),
 * (-1,1), (1,-1), (1,1) by sampling x and y independently from the
 * interval (-1,1).  
 * 
 * Second, determine if (x,y) lies within the circle circumscribed in
 * the square, by testing if (x^2 + y^2 < 1). Then multiply the
 * proportion of points inside the circle by the area of the square,
 * which is 4.
 *
 * The posterior mean of the generated quantity pi will compute pi
 * to within MC error.
 *
 * Run with fixed parameter algorithm and zero warmup iterations.
 */
model {
}
generated quantities {
  real<lower=-1,upper=1> x;
  real<lower=-1,upper=1> y;
  real<lower=0,upper=4> pi;
  x <- uniform_rng(-1,1);
  y <- uniform_rng(-1,1);
  pi <- 4 * (x^2 + y^2 < 1);
}

Generalization to volume of hypersphere

/**
 * Compute the volume of a hypersphere of a given radius in a
 * given number of dimensions by Monte Carlo integration.
 *
 * The volume will be the mean of the sapled sphere_volume values.
 */
data {
  int<lower=1> N;
  real<lower=0> r;
}
transformed data {
  real cube_volume;
  cube_volume <- (2 * r)^N;
}
model {
}
generated quantities {
  vector<lower=-1, upper=1>[N] x;
  real<lower=0,upper=cube_volume> sphere_volume;

  for (n in 1:N)
    x[n] <- uniform_rng(-1,1);

  sphere_volume <- cube_volume * (dot_self(x) < r);
}

[bob-carpenter]

• add example of WAIC calculation

[bob-carpenter]

• add example of posterior predictive checks

[bob-carpenter]

• add example of pure cross-validation
• discuss leave-one-out cross-validation and relate to WAIC (see Aki and Andrew's paper)

[bob-carpenter]

• include an example of a spline model from BDA 3 (I still don't understand the description in the book, though, so hopefully someone else can code an example)
• add Kalman filter examples (again, this is something I'm not that familiar with, but should be doable in Stan)

[bob-carpenter]

• review all the recommendations w.r.t. Cholesky factorization to bring them up to current best practice in terms of computational and statistical efficiency