# Inverse gaussian for degradation modelling

Hello there!
I’m trying to reproduce/translate the following Winbugs code from this [paper](// from Inverse Gaussian process models for degradation analysis: A Bayesian perspective - ScienceDirect), tbh really far from being a success,
Could anyone advise on any inconsistency you might spot?
Attached data and original winbugs model
data_stan.Rdata (915 Bytes)

``````// Adaptation from
// from https://www.sciencedirect.com/science/article/pii/S0951832014001276
//
functions {
// attempt to code inverse gaussian from http://paul-buerkner.github.io/brms/articles/brms_families.html
real IG_lpdf(real y, real alpha, real beta ){
real prob;
prob = pow( (alpha / 2*pi()*(pow(y,3))), 0.5) * exp( - alpha * pow( (y - beta ), 2) / (2* pow(beta,2)*y));
return prob;
}
}

data{
int<lower=1> nUnits;
int<lower=1> nTime;
matrix[nUnits, nTime] t;
matrix[nUnits, nTime] Y;
}

transformed data{
matrix[nUnits, nTime] deltaY;

for (i in 1:nUnits){
for (j in 2:nTime){
deltaY[i,j] = Y[i,j] - Y[i,j-1]; // diff
}

}
}

parameters{
real<lower=0> lambda; //
real<lower=0> q;
real<lower=0> mu_b;
real<lower=0> mu_a;
real<lower=0> mu_raw;
}

transformed parameters{
real<lower=0> mu[nUnits];
matrix[nUnits, nTime] deltaLambda;
real<lower=0> a;
real<lower=0> b;
//
//
for (i in 1:nUnits){
mu[i] = mu_a + mu_b*mu_raw;
for (j in 2:nTime){
deltaLambda[i,j] = pow(t[i,j],q) - pow(t[i,j-1],q);
a = mu[i] * deltaLambda[i,j];
b = pow(mu[i],3) *pow(deltaLambda[i,j],2) / lambda ;
}
}

}

model{
mu_b ~ normal(0,50);
mu_a ~ normal(0,50);
lambda ~ gamma(4.5, 0.1);
q ~ normal(0,10);

mu_raw ~ std_normal();

for (i in 1:nUnits){
for (j in 2:nTime){
deltaY[i,j] ~ IG(a, b);
}
}
}
``````