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);

}