@Xiang_Ye has worked on circular distributions initially for INLA, but as he was visiting at Aalto, he started making three Stan case studies that are now available and I think they are great and worth checking out by anyone working with circular valued covariates or targets. Case studies include also links to @Xiang_Ye’s papers on these topics
Thanks, @Xiang_Ye (and @avehtari for posting)! This is a great start.
We could use more on this in the User’s Guide, too, @Xiang_Ye, if you want to write something in that form. The User’s Guide is actually simpler in that we just talk about how to code the models and don’t engage in any simulation, fitting, etc. As it stands now, we don’t have much about hyper spherical statistics.
I’m looking at the second case study on regression and am wondering why it’s not using the unit vector type? Isn’t this going to run into the problem that 0 and 2 * pi are the same, but highly separated in the 1D representation?
I’m also wondering about the numerics which implies a stop gradient for extreme inputs in the second pass of conditionals:
if (d_sq < 1e-12) d_sq = 1e-12;
if (numerator < 1e-12) numerator = 1e-12;
If these ever trigger, they introduce a stop gradient and cut’s Stan’s ability to feed information from the data back to the parameters. It will devolve to rejection sampling.
I’m also not clear on why this conditional isn’t rolled into the definition of lavm_lpdf
There are also numerics that can be cleaned up. Whenever you see log(1 - alpha), it should be turned into log1m(alpha) in Stan and log1p(-alpha) in non-Stan systems.
I’d also recommend dropping the redundant documentation, like “# 1. Extract draws directly into a data frame” and “# 2. Build individual density plots” and " // 2. Custom Random Number Generator for Posterior Predictive Checks" and "// Calculate securely in log space to prevent overflow
" (also “securely” isn’t the adverb you want here—I think you want to just drop the adjective, and it also tends to prevent underflow as much as overflow). Also, the numbers aren’t giving you anything in comments as you can see the layout of the code order. Other comments are wrong, like “// The latent Gaussian process” on a variable (a variable isn’t a GP!). Also, you usually don’t want to document the language, but something like “// Local scope to prevent saving intermediate variables to MCMC draws” can be OK in a tutorial.
But the real win’s taking things like where the documentation is “PC Prior Implementation for kappa” and turning this into a named function. Then you get documentation and encapsulation for free and there’s no chance of the doc being out of synch with the code.
I don’t care so much about things like “the slope is clearly separated from zero.” I think this is leading people to think like a frequentist about whether an effect is “significant” or not.
Is there anything being done for the problem of discontinuity where 0 = 2 \pi? I didn’t see it on a quick glance, but I’ve heard this can be a serious problem for these models in practice.
Rather than " // Soft constraints to preserve identifiability", you want to use our new sum-to-zero-vector types. they’re super efficient compared to soft identifiability like this and @spinkney’s shown how to use them to back off to what you would’ve gotten without the sum-to-zero constraint.
You a usually do better than wide normal priors on a RW prior. But for this,
for (i in 3:N) {
w[i] ~ normal(2.0 * w[i-1] - w[i-2], sigma_wi);
}
you want to vectorize as we describe in the User’s Guide, to:
w[3:N] ~ normal(2 * w[2:(N - 1)] - w[1:(N - 2)], sigma_wi);
I’d just use sigma_w for the name here if you want to follow the verbose Gelman-style conventions. But given there are no other sigma variables here, I’d just drop the suffix.
One of the main things you can do to accelerate Stan code is avoid duplicating computation. So let’s look at this loop:
for (g in 1:G) {
real x = start_val + (g - 1) * step_size;
real tan_half_x = tan(x / 2.0);
real A = cos(2.0 * atan(tan_half_x - eta));
real D = 1.0 + square(eta) - eta * sin(x) - square(eta) * square(sin(x / 2.0));
// Calculate securely in log space to prevent overflow
log_probs[g] = (kappa * A) - log(D);
}
This can all be vectorized to much more efficient code:
vector[G] x = start_val + linspaced_vector(0, G - 1, 1) * step_size;
vector[G] tan_half_x = tan(0.5 * x);
vector[G] A = cos(2 * atan(tan_half_x - eta));
real square_eta = square_eta;
vector[G] D_minus_1 = square_eta - eta * sin(x) - square_eta * square(sin(0.5 * x));
vector[G] log_probs = kappa * A - log1p(D);
You can simplify
categorical_rng(softmax(log_probs));
to
categorical_logit_rng(log_probs);
(We made the same mistake as the ML world in calling the argument to softmax “logits” when they are correctly log probabilities as @Xiang_Ye coded it.
When you have five calls to ggplot with the same arguments, it’s time to pull them out into a function. The general rule in software engineering is do it once, copy and paste once, then generalize into a function rather than copy and paste a third time.
I’d say that whenever you’re tempted to do something like this:
vector[N] w_unc;
w_unc[1] = z_w[1] * 5.0;
w_unc[2] = z_w[2] * 5.0;
it’s better to use an appropriate prior or likelihood on z_w to get the same effect. It can also be coded more compactly (though not really more efficiently) as
w_unc[1:2] = 5 * z_w[1:2];
(You don’t need all the .0 markings in Stan—the compiler in C++ figures out whether it needs to be promoted to floating point at compile time, so there’s no efficiency overhead.)
I’d also add a conclusion of some kind.
Anyway, the only reason I added all these notes is that I’m excited to see this revised and maybe something like it added to the User’s Guide. So thanks again. We can also publish on the Stan web site under case studies (I think @avehtari probably has enough privileges for that, but if not, you can ask on the forums here and we’ll find someone to post).
Hi Professor Carpenter,
Thank you so much for taking the time to go through the tutorials and provide so many detailed comments. I really appreciate it.
I don’t usually work with MCMC or Stan directly, so I’m not as experienced with many of the Stan-specific implementation and numerical practices. Your comments are therefore extremely helpful. I’ll go through them carefully, check the issues you pointed out, and revise the tutorials accordingly.
I’ll reach out again after I’ve had a chance to read through everything carefully and make the revisions. Thanks again!
Let’s hope you won’t be able to say that for much longer.
Aki should also be able to help with all of the MCMC and Stan parts of this. Also, I don’t really know much at all about hyper spherical stats other than that you calculate means and variances differently.
It’d be great to get a version of this suitable for the User’s Guide, but I’d recommend finishing the case studies first so it’s clear what to put in the guide. I’m happy to provide feedback on revisions if you’d like. But it’s also useful to just get stuff out there—anything can be made better.
P.S. Everyone’s on a first name basis here, so please just “Bob”! I haven’t been a professor since 1996!
Thanks, Bob! I’ve now made several revisions based on your suggestions. They’ve all been very helpful. I agree that refining the case studies first makes sense.
Regarding the unit-vector and boundary questions, our aim is to incorporate circular variables as responses or covariates into a general regression framework with real-valued predictors and latent effects. LAvM uses a local real-valued representation of the circular variable to form the residual
g^{-1}(y_i)-\eta_i
=
\tan(y_i/2)-\eta_i,
\qquad
\eta_i=\mathbf{x}_i^\top\boldsymbol{\beta}+w_i,
before mapping it back into the circular density.
The inverse link maps (-\pi,\pi) to \mathbb{R}, with the endpoint limits corresponding to -\infty and +\infty. The circular point -\pi(or \pi) is therefore singular for this representation. This is a trade-off of this particular construction: we obtain the desired linear-scale regression structure, but restrict its intended use to reasonably concentrated data, rather than approximately uniform data, with a suitable angular origin that keeps the relevant observations and fitted directions away from that point. Concentration alone is insufficient if it straddles the chosen cut.
I see the advantage of a unit vector for a freely estimated circular location, particularly when its posterior crosses the cut. Here, however, the regression is parameterized through real-valued coefficients and latent effects. A unit-vector representation would not by itself remove the inverse-link singularity while preserving this construction, so I would be interested in your thoughts on its role here.
For a simple circular location model, ordinary von Mises already covers the entire circle and may suffice without the LAvM adjustment, while the PC prior remains applicable. LAvM is instead motivated by the potential ambiguity in the usual linked-location regression formulation discussed in the tutorial and paper, rather than a need to replace von Mises for circular location estimation in general.
I hope this explanation is clear and that I have understood your concern correctly.
Hi @Xiang_Ye, really cool stuff! You might be interested in this alternative parameterization of the unit vector A better-er unit vector where you could put regression coefficients and then run the transform over the regression vector to output the mean unit-vector for the Von Mises distribution.
Yes, and thanks for responding. I just saw this now.
Also, I’d check out anything @spinkney is doing—he’s been massively improving the geometry of our transforms.