Writing a user defined function for inverse-Gaussian distribution

Modeling
1

Hi,

I am trying to fit the inverse-Gaussian distribution to my failure data. I wrote a function to define the inverse-Gaussian distribution and tried to get posterior samples. I have no experience in defining new distributions in RStan. My codes produced the following error.

SYNTAX ERROR, MESSAGE(S) FROM PARSER:

 No matches for: 

  real[] ~ IG(real, real)

 Available argument signatures for IG:

   vector ~ IG(real, real)


 ERROR at line 31

 30:    //Likelihood of the data
 31:    y ~ IG(mu, lambda);
                       ^
 32:    }

Error in stanc(file = file, model_code = model_code, model_name = model_name,  : 
  failed to parse Stan model 'a49daf541cb95cd29896bb541eb39aa1' due to the above error.

I found an implementation of inverse-Gaussian distribution here, which worked finely with my data. But I’d like to fix the errors in my Stan program.

The codes I used to define the distribution are,

failure_data <- list (J = 88,
                  y = rep (c (8, 16, 32, 40, 56, 60, 64, 72, 80, 96, 104, 108, 112, 114, 
                              120, 128, 136, 152, 156, 160, 168, 176, 184, 194, 208, 216, 
                              224, 232, 240, 246, 256, 264, 272, 280, 288, 304, 308, 328,
                              340, 352, 358, 360, 384, 392, 400, 424, 438, 448, 464, 480, 
                              536, 552, 576, 608, 656, 716), times = c (1, 4, 2, 4, 3, 1, 1, 
                                                                        5, 4, 2, 1, 1, 2, 1, 1, 1, 1, 3, 1, 1, 5, 1, 3, 1, 2, 1, 4, 1,
                                                                        1, 1, 1, 2, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
                                                                        1, 1, 1, 1, 1, 1, 1)) 
)
#pdf of inverse Gaussian distribution
#f(x; mu, shape) = (shape/(2*pi*x^3))^0.5*exp(-shape*(x - mu)^2/(2*mu^2*x))

model_IG <- "functions {

  real IG_log (vector x, real mu, real lambda){
vector [num_elements (x)] prob;
real lprob;
for (i in 1:num_elements(x)) {
  prob[i] <- (lambda/(2*pi()*(x[i]^3)))^0.5*exp(-lambda*(x[i] - mu)^2/(2*mu^2*x[i]));
}
lprob <- sum (log(prob)); 
return lprob;

}
}

data {
int<lower=0> J; // number of observations
real y[J]; // failure in hours

}
parameters {
real<lower=0> mu; // mean
real<lower=0> lambda; // shape

}

model {
//define priors of mu and lambda


//Likelihood of the data
y ~ IG(mu, lambda);
}
"

mod_IG <- stan_model(model_code = model_IG)

fit_IG <- sampling (mod_IG, data = failure_data, chains = 4, iter = 10^4,
                   warmup = 5000, thin = 20)
fit_IG

Thank you in advance for any help and suggestions to improve my coding.

2

It looks like the function is defined with a vector as input data (“vector x”), but you’re passing it an array of reals. You could either change the function code, or change your data statement from “real y[J]” to “vector[J] y”.

3

Thank you. Now it works well!!