Option to keep constant terms in log probability via standard sampling syntax?

Developers
1

Dear Stan developers,

Question. Might it be possible as a future development to have a run-time option to have standard sampling statements, like y ~ normal(mu, sigma), sample in a way that keeps the constant aspects of the log probability (called e.g. keep_constants_in_lp), like the target += ... notation does now?

Rationale. Sampling statements like y ~ normal(mu, sigma) drop constant terms, while “target” notation like target += normal_lpdf(mu, sigma) keeps them [1]. Some applications like bridge sampling require these constants [2, 3]. So, Stan code needs to be rewritten for that purpose. The code gets ugly fast, and this increases the barriers to entry and the chance of errors. Some of those difficulties relate to the relative lack of flexibility of “target” notation. For example, if you want to sample with boundaries, e.g. vector<lower=a, upper=b> y, then you also have to adjust target using <distribution>_lcdf(upper | parameters) etc. [2]. If you try to automate this, to reduce the chance of programming errors, you run into the problem that Stan functions can’t have the same names but different signatures – my approach has been to write code [4] to write a bunch of Stan functions [5] that handle sampling with “target” notation. In contrast, ~ notation handles lots of variable types automatically with the same notation, which is clean and simple.

Would this be feasible in Stan itself? That is, if keep_constants_in_lp is false, it would works as it does now (for speed). If keep_constants_in_lp is true, it would behave like target += notation does now, and if it’s sampling a variable with explicit lower and/or upper limits, it would apply the correction(s). That would allow much simpler syntax for this sort of situation.

Thanks as ever for a fantastic piece of software!

all the best,
Rudolf.

[1] 7.4 Sampling Statements | Stan Reference Manual
[2] CRAN - Package bridgesampling
[3] https://doi.org/10.1016/j.jmp.2017.09.005
[4] e.g. https://egret.psychol.cam.ac.uk/rlib/make_commonfunc_stan.py
[5] e.g. https://egret.psychol.cam.ac.uk/rlib/commonfunc.stan

2

Hi Rudolf,

you might want to look into the _lupdf suffix that is available in Stan 2.25 and later. This might help with your use-case.

Consider this example:

functions {
    real foo_lpdf(real y, real x) {
        return normal_lupdf(y| x, 1);
    }
    real goo_lpdf(real y, real x) {
        return normal_lpdf(y| x, 1);
    }
}
parameters {
    real y;
}
model{
    target += foo_lpdf(y| x); // normal inherits propto=false
    target += foo_lupdf(y| x); // normal inherits propto=true
    target += goo_lpdf(y| x); // normal uses propto = false
    target += goo_lupdf(y| x); // normal uses propto = false
}

So you can write you entire likelihood with lupdf suffixes and then call you user-defined lpdf in the model block. If you call with lupdf, you drop the constants, and if you call it with lpdf the constants will not be dropped.

2 Likes
3

Also, I noticed that you have a comment in your collection of functions: “size() doesn’t work on a plain “vector”. Use num_elements().”.

size() works with all types as of at least Stan 2.24.

2 Likes
4

Thanks, Rok! Good to know, both re _lupdf() and the update to size().

The net syntax is still significantly more cumbersome, though; that was part of the reason for my request. For example,

parameters {
    vector<lower=LOWER, upper=UPPER>[N] y;
}
model {
    y ~ normal(MU, SIGMA);  // drop constants
}

versus

parameters {
    vector<lower=LOWER, upper=UPPER>[N] y;
}
model {
    target +=
        normal_lpdf(y | MU, SIGMA) -  // don't drop constants
        N * log_diff_exp(  // adjust each for truncation
            normal_lcdf(UPPER | MU, SIGMA),
            normal_lcdf(LOWER | MU, SIGMA)
        );
}

With a dozen or so of those, the potential for error increases (and the user-friendliness decreases).

1 Like
5

Oh, in case of truncation, the target += syntax is much more cumbersome, you are right.
I see that rarely and forgot about that.

1 Like