# Computing pdf values of a distribution in STAN

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

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.

``````  for (i in 1:nt)