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.