Using a built-in distributions in customised function

Modeling
1

Hi! is my first comment on the group. I want to congratulate everyone for a such a remarkable work!
This is my question:
I want to build a customised likelihood following this paper:
https://link.springer.com/article/10.3758/s13428-016-0746-9
I want to use the built-in gamma distribution in my function but I don’t know if it is legal, the PARSER shows the
error : expression not assignable to the left side. You will forgive me for not giving the full script but I attach the function in particular!

  real ptcm_log(vector yobs,  vector theta, matrix X,  real alpha, real sigma){
    vector[num_elements(yobs)] lpdf;
    vector[num_elements(yobs)] pdf;
    vector[num_elements(yobs)] cdf;
    real term1;
    real term2;
    vector[num_elements(yobs)] prob;
    vector[num_elements(yobs)] p;
    real lprob;
    
    p = exp(-exp(X * theta));
    lpdf = gammal_lpdf(yobs | alpha, beta);
    pdf = exp(lpdf);
    cdf = gamma_cdf(yobs , alpha, beta);
    
    for(i in 1:num_elements(yobs)) {
        term1 = (-log(p[i])*pdf[i]);
        term2 = exp(log(p[i])*cdf[i]);
        prob[i] = term1 * term2;
    }
    lprob = sum(log(prob));
    return lprob;
  }
2

Mind pasting the parser error? It makes it much easier for someone to read and help if the parser error is already provided.

P.S. I edited your original post by adding code blocks so the function is easier to read.

3

I think I found a way out which is using a not vectorised form of the built-in function in the following way:

for(i in 1:num_elements(yobs)) {
    lpdf = gammal_lpdf(yobs[i] | alpha, beta);
    pdf = exp(lpdf);
    cdf = gamma_cdf(yobs[i] , alpha, beta);
    term1 = (-log(p[i])*pdf[i]);
    term2 = exp(log(p[i])*cdf[i]);
    prob[i] = term1 * term2;
}
lprob = sum(log(prob));
return lprob;

}

I don’t know why it does not work with the vectorised version.
Thanks very much,

4

Great, that’s the way you need to do currently.

It’s a known design issue. When that function was designed, the designer was thinking the case where it’s used to increment log density and then it’s handy that it returns sum instead of a vector. There has been discussions how to change the design to have a different name or some option to allow vector output, too.

5

Hi,
Alright thank you very much for the clarification.