Self-defined joint likelihood function


I’m trying to code using package “rstan” in R my own joint likelihood function for the MCMC sampling to run on. Currently I have:

real newexp_lpdf(vector x, real r, real mu){
vector[num_elements(x)] prob;
real lprob;
for (i in 1:num_elements(x)){

prob[i] = -r/mu*(gamma_p(r+1,(x[i]+1)*mu) ;


lprob = sum(log(prob));
return lprob;

data {
int N;
vector[N] Y;
parameters {
real <lower=0> mu;
real <lower=0> r;
model {
real alpha;
real beta;
alpha = 1.0;
beta = 1.0;
mu ~ gamma(alpha, beta);
r ~ gamma(alpha, beta);
Y ~ newexp(r,mu);

My question is, if my observed data Y contains multiple data points per observation. That is, instead of observing (Y1, Y2,…Yn), observe (A1, B1, C1, D1), (A2, B2, C2, D2), …(An, Bn, Cn, Dn). Is there a way I can code the observed “vectors” as:

vector~ newexp(r, mu)

Instead of Y~ newexp(r,mu)?

Thanks for your help.


Y ~ newexp(r, mu) isn’t going to work because your function defines its first argument as a vector instead of an array of vectors. Looping over observations and doing Y[i] ~ newexp(r, mu) could work but I haven’t verified the details of this model.