How to use maximum function in my model?

Hi, everyone. Recently I have met a problem in my project and I will show it in short.

Suppose y=b_1*\max(x-\gamma,0)+b_2*\min(x-\gamma,0)+\epsilon with \epsilon \sim Normal(0,\sigma^2), and now i have dataset (x_i,y_i). I want to estimate 4 parameters b_1,b_2,\gamma and \sigma^2. How should I use maximum or minimum functions in my Stan code?

data {
    int<lower=0> N;
    vector[N] y;
    vector[N] x;
}
parameters{ 
    real<lower=0> sigma;
    real a1;
    real a2;
    real gamma;
}
model{
    for (n in 1:N)
        y[n] ~ normal(a1*max(x[n]-gamma,0)+a2*min(x[n]-gamma,0),sigma);
}

or the last line can be changed by

        y[n] ~ normal(a1* (x[n]-gamma+abs(T x[n]-gamma))/2
                            +a2* (x[n]-gamma- abs(T x[n]-gamma))/2
                             ,sigma);

The maximum and minimum functions are documented here: https://mc-stan.org/docs/2_20/functions-reference/step-functions.html

However, I suspect that what you are seeking to do will not work very well in Stan because when fmin and fmax are applied to parameters, they can introduce discontinuities in the gradients used by Stan’s algorithms.

2 Likes