Puzzling parsing error with parenthesis

Modeling
1

I am trying to compile a Stan program but get an error that I don’t understand:

SYNTAX ERROR, MESSAGE(S) FROM PARSER:
 error in 'modelbc551814461b_stan_bc5560116d23' at line 157, column 24
  -------------------------------------------------
   155:     qtilde3 ~ weighted_hat_lpdf(knots3, weights3, N3);
   156:  
   157:     target += gpoispoint(y | mu, xi, sigma,  u, noy) ;
                               ^
   158:  
  -------------------------------------------------

PARSER EXPECTED: "("
Error in stanc(file = file, model_code = model_code, model_name = model_name,  : 
  failed to parse Stan model 'stan-bc5560116d23' due to the above error.

The parser seems not to identify the “(” that is clearly present. What is going on please?

I give below the complete Stan program:

functions {

    real gpoispoint_lpdf(vector y, real mu, real k, real sigma, real u, real noy) {
        // generalised Pareto log pdf
        int N = rows(y);
        real tester;
        if (sigma<=0){
            reject("sigma<=0; found sigma =", sigma);
        }
        if (fabs(k) > 1e-10) {
            tester = min(1 + (y - mu) * (k / sigma));
            if(tester < 0){
              reject("parameters out of range", mu, sigma, k);
              return(0);
            }
            else{
              return -(1+inv(k))*sum(log1p((y-mu) * (k/sigma))) -N*log(sigma) -
            noy * (1 + k * (u - mu) / sigma) ^ -inv(k);
            }
        }
        else {
            return -sum(y-mu)/sigma -N*log(sigma) -
            (noy * exp(u - mu) / sigma); // limit k->0
        }
    }

     
    real felt_p1(real x){
        return (1 - cos(pi() * (x + 1))) / 2;
    }

    real felt_p2(real x) {
        return (1 + cos(pi() * x)) / 2;
    }

    real felt_c1(real x) {
        return (x + sin(pi() * x) / pi()) / 2;
    }

    real felt_c2(real x) {
        return (x + sin(pi() * x) / pi()) / 2;
    }

    real cf(real x, real lo, real mid, real hi){
        if (x < lo)
            return 0;
        else if (x <= mid)
            return (mid - lo) * (felt_c1((x - mid)/(mid - lo)) - felt_c1(-1)) ;
        else if (x  < hi)
            return (hi - mid) * (felt_c2((x - mid)/(hi - mid))) + (mid - lo) / 2;
        else
            return (hi - lo) / 2;
    }

    real pf(real x, real lo, real mid, real hi){
        real x1 = (x - mid)/(mid - lo);
        real x2 = (x - mid) / (hi - mid);
        if ((x < lo) || (x > hi))
            return 0;
        else if (x <= mid)
            return felt_p1(x1);
        else
            return felt_p2(x2);
    }


    real wted_cdf(real x, int j, vector knots, vector weights, int N){
        real lo;
        real mid = knots[j];
        real hi;
        if(j == 1)
            lo = 0;
        else
            lo = knots[j - 1];
        if(j == N)
            hi = 1;
        else
            hi = knots[j + 1];
        return weights[j] * cf(x , lo, mid, hi);
    }

    real wted_pdf(real x, int j, vector knots, vector weights, int N){
        real lo;
        real mid = knots[j];
        real hi;
        if (j == 1)
            lo = 0;
        else
            lo = knots[j - 1];
        if (j == N)
            hi = 1;
        else
            hi = knots[j + 1];
        return weights[j] * pf(x , lo, mid, hi);
    }

    real weighted_hat_lpdf(real y, vector knots, vector weights, int N){
        vector[N] pfn;
        real z = y / (y + 1);
        for (n in 1:N){
            pfn[n] = wted_pdf(z, n, knots, weights, N);
        }
        return log(sum(pfn));
    }
}

data{
    // parameters of distribution 1
    int<lower = 0> N1;
    vector[N1] knots1;
    vector[N1] weights1;
    // parameters of distribution 1
    int<lower = 0> N2;
    vector[N2] knots2;
    vector[N2] weights2;
    // parameters of distribution 1
    int<lower = 0> N3;
    vector[N3] knots3;
    vector[N3] weights3;

    int<lower=0> N; // number of observations
    vector[N] y; // observations
    real noy;

    real u; // high threshold

}
parameters {
    real<lower = 0> qtilde1;
    real<lower = 0>  qtilde2;
    real<lower = 0>  qtilde3;
}
transformed parameters{
    real mu;
    real xi;
    real sigma;
    if(fabs(qtilde3 - qtilde2) / qtilde3 > 1e-10){ 
     real nu = qtilde3 / qtilde2;
     xi = log10(nu);
     sigma  = qtilde2 * xi / (nu * (nu - 1));
     mu = qtilde1 - qtilde2 / nu;
    } 
    else {
      xi = 0;
      mu = qtilde1 - qtilde2;
      sigma = qtilde2 / log(10);
    }
}
model {
    // priors
    qtilde1 ~ weighted_hat_lpdf(knots1, weights1, N1);
  
    qtilde2 ~ weighted_hat_lpdf(knots2, weights2, N2);
 
    qtilde3 ~ weighted_hat_lpdf(knots3, weights3, N3);
 
    target += gpoispoint(y | mu, xi, sigma,  u, noy) ;
 

}
2

While this might not be the main issue, something that jumps out is you don’t need the _lpdf suffix when using the ~ statement, but you do need it when using target +=. (So the opposite of what you’re currently using).

Also, the two forms will be equivalent for user-defined lpdfs

3

Thanks. I’ll see if that helps.
But I don’t understand your final sentence. Could you put it another way?

4

The two approaches (~ vs target +=) are used to specify whether constants should be included in the calculation of the log-probability (see this section of the manual for more information).

This is handled internally in c++ by detecting whether the inputs to a given distribution (e.g. normal) are parameters or just data. In contrast, a user-defined lpdf will always return the same value regardless of which sampling notation is used.

5

Thanks again.
Unfortunately my program still does not run.

6

If I make the suggested change, the model parses without issue for me. Reprex below:

modcode = '
functions {

    real gpoispoint_lpdf(vector y, real mu, real k, real sigma, real u, real noy) {
        // generalised Pareto log pdf
        int N = rows(y);
        real tester;
        if (sigma<=0){
            reject("sigma<=0; found sigma =", sigma);
        }
        if (fabs(k) > 1e-10) {
            tester = min(1 + (y - mu) * (k / sigma));
            if(tester < 0){
              reject("parameters out of range", mu, sigma, k);
              return(0);
            }
            else{
              return -(1+inv(k))*sum(log1p((y-mu) * (k/sigma))) -N*log(sigma) -
            noy * (1 + k * (u - mu) / sigma) ^ -inv(k);
            }
        }
        else {
            return -sum(y-mu)/sigma -N*log(sigma) -
            (noy * exp(u - mu) / sigma); // limit k->0
        }
    }

     
    real felt_p1(real x){
        return (1 - cos(pi() * (x + 1))) / 2;
    }

    real felt_p2(real x) {
        return (1 + cos(pi() * x)) / 2;
    }

    real felt_c1(real x) {
        return (x + sin(pi() * x) / pi()) / 2;
    }

    real felt_c2(real x) {
        return (x + sin(pi() * x) / pi()) / 2;
    }

    real cf(real x, real lo, real mid, real hi){
        if (x < lo)
            return 0;
        else if (x <= mid)
            return (mid - lo) * (felt_c1((x - mid)/(mid - lo)) - felt_c1(-1)) ;
        else if (x  < hi)
            return (hi - mid) * (felt_c2((x - mid)/(hi - mid))) + (mid - lo) / 2;
        else
            return (hi - lo) / 2;
    }

    real pf(real x, real lo, real mid, real hi){
        real x1 = (x - mid)/(mid - lo);
        real x2 = (x - mid) / (hi - mid);
        if ((x < lo) || (x > hi))
            return 0;
        else if (x <= mid)
            return felt_p1(x1);
        else
            return felt_p2(x2);
    }

    real wted_cdf(real x, int j, vector knots, vector weights, int N){
        real lo;
        real mid = knots[j];
        real hi;
        if(j == 1)
            lo = 0;
        else
            lo = knots[j - 1];
        if(j == N)
            hi = 1;
        else
            hi = knots[j + 1];
        return weights[j] * cf(x , lo, mid, hi);
    }

    real wted_pdf(real x, int j, vector knots, vector weights, int N){
        real lo;
        real mid = knots[j];
        real hi;
        if (j == 1)
            lo = 0;
        else
            lo = knots[j - 1];
        if (j == N)
            hi = 1;
        else
            hi = knots[j + 1];
        return weights[j] * pf(x , lo, mid, hi);
    }

    real weighted_hat_lpdf(real y, vector knots, vector weights, int N){
        vector[N] pfn;
        real z = y / (y + 1);
        for (n in 1:N){
            pfn[n] = wted_pdf(z, n, knots, weights, N);
        }
        return log(sum(pfn));
    }
}

data{
    // parameters of distribution 1
    int<lower = 0> N1;
    vector[N1] knots1;
    vector[N1] weights1;
    // parameters of distribution 1
    int<lower = 0> N2;
    vector[N2] knots2;
    vector[N2] weights2;
    // parameters of distribution 1
    int<lower = 0> N3;
    vector[N3] knots3;
    vector[N3] weights3;

    int<lower=0> N; // number of observations
    vector[N] y; // observations
    real noy;

    real u; // high threshold

}
parameters {
    real<lower = 0> qtilde1;
    real<lower = 0>  qtilde2;
    real<lower = 0>  qtilde3;
}
transformed parameters{
    real mu;
    real xi;
    real sigma;
    if(fabs(qtilde3 - qtilde2) / qtilde3 > 1e-10){ 
     real nu = qtilde3 / qtilde2;
     xi = log10(nu);
     sigma  = qtilde2 * xi / (nu * (nu - 1));
     mu = qtilde1 - qtilde2 / nu;
    } 
    else {
      xi = 0;
      mu = qtilde1 - qtilde2;
      sigma = qtilde2 / log(10);
    }
}
model {
    // priors
    qtilde1 ~ weighted_hat(knots1, weights1, N1);
  
    qtilde2 ~ weighted_hat(knots2, weights2, N2);
 
    qtilde3 ~ weighted_hat(knots3, weights3, N3);
 
    target += gpoispoint_lpdf(y | mu, xi, sigma,  u, noy) ;
 

}
'

mod_transp = rstan::stanc(model_code = modcode)
mod_transp$status
#> [1] TRUE

Created on 2021-10-15 by the reprex package (v2.0.1)