Computing pdf values of a distribution in STAN

Modeling
1

The distribution of a variable y over time looks similar to a gamma distribution (additionally, including a baseline and an exponential decay). So, I would like to estimate the functional form of this variable, where I only consider a specific time span x. In order to do so, I would like to directly compute the pdf values of the gamma and the epxonential distribution in STAN for the respective values x. Thus, I would actually need the equivalent to the R-function dgamma() and dexp() in STAN (respective part of the code is in bold).

Can anyone tell me a way to do this? Thanks!

data {
int<lower=1> nt;
int<lower=1> n;
vector[nt] x;
vector[nt] y;
real<lower=0> meany;
real<lower=0> sdy;
}
parameters {
real<lower=0> alpha;
real<lower=0> beta;
real<lower=0> gamma;
real a;
real<lower=0> b;
real<lower=0> c;
real<lower=0> hyperalpha;
real<lower=0> hyperbeta;
real<lower=0> hypergamma;
real<lower=0> hyperb;
real<lower=0> hyperc;
}
model {
vector[nt] y_hat;
hyperalpha ~ exponential(0.01);
hyperbeta ~ exponential(0.01);
hypergamma ~ exponential(0.01);
hyperb ~ exponential(0.01);
hyperc ~ exponential(0.01);
alpha ~ gamma(5,5/hyperalpha);
beta ~ gamma(5,5/hyperbeta);
gamma ~ gamma(5,5/hypergamma);
a ~ normal(meany1, sdy1);
b ~ gamma(5,5/hyperb);
c ~ gamma(5,5/hyperc);
y_hat = a + b * dgamma(x, alpha, beta) + c * x * dexp(x, gamma);
y ~ normal(y_hat,10);
}

2

Those would be exp() of gamma_lpdf and exponential_lpdf, which could have been found by evaluating rstan::lookup("dgamma"). However, you have to loop over the observations and evaluate them with scalars.

3

Thank you very much for your reply!

However, I am not sure how I can “evaluate them with scalars” - would the following line of code be enough?

y_hat = a + b * for(i in 1:nt){exp(gamma_lpdf(x[i] | alpha, beta))} + c * x * for(i in 1:nt){exp(exponential_lpdf(x[i] | gamma));

4

Not quite. Stan evaluates like any language based on C, so you need

  for (i in 1:nt)
    y_hat[i] = a + b * exp(gamma_lpdf(x[i] | alpha, beta)) + 
               c * x[i] * exp(exponential_lpdf(x[i] | gamma));
5

Thank you so much - this helps me a lot!