Inverse gaussian for degradation modelling

Modeling
1

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)

image
image892×619 75.9 KB 

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