Mixed distribution: coding custom lpdf

Modeling
1

I am trying to write a custom lpdf, but am having trouble working out how to code it in Stan.

In words, the value of the log pdf f(x) depends on which of three cases applies, as shown, where x, y and z are all scalars:
(i) x < y < z: f(x) = a(x),
(ii) y < x < z: f(x) = b(x), or
(iii) y < z < x: f(x) = c(x)
where the functions a, b, c are valid log pdf ( normalised if necessary so the the resultant f(x) will be a valid log pdf).

I can write code for this in R, with which I am very familiar, but whatever I try in Stan raises errors that I do not know how to deal with.

For example

functions{
  real a(real x){
    return 0.1;
  }
  real b(real x){
    return 0.1;
  }
  real c(real x){
    return 0.1;
  }
  real my_lpdf(real[] x, real y, real z){
  
    int N = size(x);
    vector[N] f; 
    for (i in 1:N){
      if (x[i] <= y) {
        f[i] = a(x[i]);
      } else if (x[i] <= z){
        f[i] = b(x[i]);
      } else {
        f[i] = c(x[i]);
      }
    }
  return f;  
  }
}

fails, with the message:


Ill-typed arguments supplied to infix operator <=. Available signatures: 
(int, int) => int
(int, real) => int
(real, int) => int
(real, real) => int
Instead supplied arguments of incompatible type: real, array[] real.

Obviously my functions a, b,c are wrong but were needed to get code that could be checked. That code reflects how I would code it in R.
The underlying problem is that as a very occasional user of Stan I am not used to its much tighter discipline over types etc. than R.

How should I code this function in Stan?

2

It looks like this is the line throwing the error. x[I] is a real while x is an array of reals.

3

Oops! You are right that it is an error but it was my carelessness in putting this topic together. I have corrected that mistake with an edit.

The original error is still in there somewhere.

4

Could you share your full model code? That might reveal any issues.

I also think you need to change your return type to a vector rather than real, since f is a vector.

5

I have made some progress but have now run into another error. Here is my code of a user-defined lpdf function which uses a second user-defined function.

functions{
  real minpdf( real z, real[,] v){
    real result;
    for (j in 1:4) {
       if (z <= v[j + 1, 1]){ 
          result = v[j, 2] + 
            (v[j + 1,2]-v[j, 2]) * (z-v[j, 1]) / (v[j + 1, 1]-v[j, 1]); 
          break;
       }  
     } 
     return(result);
  }

  real minknot_lpdf(vector y, real a, real b, real t1, real t2, real w){
    int N = rows(y);
    array[5,2] real v;
    real A;
    real B;
    real T1;
    real T2;
    real W;
    real C;
    real D;
    real E; 
    vector[N] f;
    int k;
    v[1, 1] = a;
    v[5, 1] = b;
    v[1, 2] = 0;
    v[5, 2] = 0;
    if (w < t1) {
      A = t1 - a;
      B = t1 - w;
      C = t2 - t1;
      D = b - t2; 
      T1 = (2.0/3.0) * (1/C - 1/D);
      T2 = (2.0/3.0) * (1/D);
      W = (2.0/3.0) * (1 - (B/C) + (B/D)) * (1/A);
      
      v[2, 1] = w;
      v[3, 1] = t1;
      v[4, 1] = t2;
      v[2, 2] = W;
      v[3, 2] = T1;
      v[4, 2] = T2;
      
  } else if (w < t2) {
    A = t2 - t1;
    B = w - t1;
    C = t2 - w;
    D = b - t2;
    E = t1 - a;
    
    T1 = (2.0/3.0) * (1/E);
    T2 = (2.0/3.0) * (1/D);
    W = (2.0/3.0) * (1 - (B/E) - (C/D)) * (1/A);
    
    v[2, 1] = t1;
    v[3, 1] = w;
    v[4, 1] = t2;
    v[2, 2] = T1;
    v[3, 2] = W;
    v[4, 2] = T2;
  } else {
    A = b - t2;
    B = w - t2;
    C = t2 - t1;
    D = t1 - a;
    
    T2 = (2.0/3.0) * (1/C - 1/D);
    T1 = (2.0/3.0) * (1/D);
    W = (2.0/3.0) * (1 - (B/C) + (B/D)) * (1/A);
    
    v[2, 1] = t1;
    v[3, 1] = t2;
    v[4, 1] = w;
    v[2, 2] = T1;
    v[3, 2] = T2;
    v[4, 2] = W;
     
  }
  for (j in 1:N){
    f[j] = minpdf(y[j], v)
  }
      
  return log(f);  
  
   }
}

When I check it I get:

rstan:::rstudio_stanc("minknot1.stan")
Error in stanc(filename, allow_undefined = TRUE) : 0
Syntax error in 'string', line 83, column 2 to column 3, parsing error:
   -------------------------------------------------
    81:    for (j in 1:N){
    82:      f[j] = minpdf(y[j], v)
    83:    }
           ^
    84:  
    85:    return log(f);
   -------------------------------------------------

Ill-formed phrase. Found L-value "=" expression. There are many ways in which this can be completed to a valid phrase.

If you cannot help, do I need to post this as a new query?

6

You’re missing a semicolon after the minpdf call

3 Likes
7

Thanks. It’s the influence of R again.

1 Like