Caching results of repeated evaluations of integrate_1d in custom distribution?

Modeling
1

Hi,

I am defining a custom distribution (i.e. a user defined *_lpdf function) which makes use of integrate_1d to numerically solve for the normalization constant. This is a 1D integration over y for a distribution whose shape (sigma), location (y_0), etc. parameters are parameters of the Stan model. While debugging, I added some print statements and noticed that the integrate_1d function seems to be called many times for the same parameter values. Presumably this happens because the *_lpdf function is called once per data sample, where the normalization will be constant because the parameters are not changing within the leapfrog step.

In a situation like this, will Stan automatically cache the integrator output so that it can skip repeated integrate_1d evaluations? If not, is it possible to manually do that in some way?

In particular, I wondered if I could extract the integration out of the *_lpdf function so that I could manually enforce that it be run just once per step, but I wasn’t sure how to do that while still making use of sampling statements (with the ~ operator) associated with a custom distribution.

Thanks in advance for any guidance!

PS it looks like the link to this “performance” category in the welcome page is broken.

2

There’s no such caching that I’m aware of. Presumably you’re sure you need the normalization constant?

Can you post your code? If your lpdf function doesn’t literally have it inside a loop over data points, and it’s not being called inside such a loop in the model block, it really shouldn’t be computed repeatedly on a given leapfrog step as you report.

3

Hi @mike-lawrence, thanks for following up. Here are my work-in-progress custom functions,

    real unnorm_weibull_logistic_pdf(real y, // Function argument
                                     real xc, // Complement of function argument on the domain (not explicitly defined in model code)
                                     real[] theta, // parameters
                                     data real[] x_r, // data (real)
                                     data int[] x_i // data (integer)
                                     ) {
        // This is an un-normalized logistic-modified Weibull distribution. It is used as an integrand
        // to determine the normalization.
        real alpha = theta[1];
        real sigma = theta[2];
        real kappa = theta[3];
        real y_0 = theta[4];
        real b = x_r[1];
        int y_min = x_i[1];
        real log_wb;
        if (y >= y_min) {
            log_wb = weibull_lpdf(y | alpha, sigma);
        }
        else {
            log_wb = negative_infinity();
        };
        real log_lg = log(b) - log(1 + exp(-kappa * (y - y_0)));
        return exp(log_wb + log_lg);
    }

    real weibull_logistic_lpdf(real y, real alpha, real sigma, data real b, real kappa, real y_0, data int y_min) {
        // This is a logistic-modified, left-truncated Weibull distribution.  The normalization is
        // derived by numeric integration.
        real wb = weibull_lpdf(y | alpha, sigma);
        // Adjust the Weibull to account for truncation
        if (y >= y_min) {
            wb += -weibull_lccdf(y_min | alpha, sigma);
        } else {
            wb += negative_infinity();
        }
        real lg = log(b) - log(1 + exp(-kappa * (y - y_0)));
//         print(" alpha ", alpha, " sigma ", sigma, " kappa ", kappa, " y_0 ", y_0);
        real norm = integrate_1d(unnorm_weibull_logistic_pdf,
                            y_min * 1.0,
                            positive_infinity(),
                            {alpha, sigma, kappa, y_0}, {b}, {y_min},
                            1e-4);
        return -log(norm) + wb + lg;
    }

This is invoked in the model block like this:

    for (n in 1:N) {
        k = K_id[n];
        y_real[n] ~ weibull_logistic(alpha[k], sigma[k], 1.0, kappa[k], y_0[k], y_min);
    }

Where alpha, sigma, kappa, and y_0 are Stan model parameters.

4

Ah, then you are looping over the data in the model block. Put the loop over data in your function instead and compute the normalization constant outside the loop.

5

Hi @mike-lawrence would that approach entail directly adding to target in the model block rather than using a sampling statement (~)? I can understand the reason for doing that, if so, but want to make sure I go about this in the right way.

6

Nope, just define your lpdf function as taking vector inputs rather than reals

1 Like
7

Hi Mike,

I gave this a try last night and it seems to work well! Thanks so much for the guidance. I’ll plan to write up some examples and post some updated code on GitHub for anyone interested in what the result looks like.

~Nathan

1 Like