Estimate Shape and Scale for Weibull distribution from the mean and standard deviation

General
,
1

Hello,
I am trying to write a function in Stan that will let me convert a mean and standard deviation into the shape and scale parameters expected by the weibull sampling statement and associated functions.

I can do this in R using the mixdist package like:

mu = 1/runif(2, min=0, max=2)
#0.6216184 0.6359826
weibullpar(mu, sigma=2)
#   shape     scale loc
#1 0.3942742 0.1796633   0
#2 0.3997025 0.1909753   0

I’ve tried to translate the R code into a stan function like:

functions {
  real weibullshape(real capmu, real capsig){
    real cv;
    real shape;
    real scale;
    cv = capsig/capmu;
    real aa;
    real gb;
    real nu;
              aa = log(cv^2 + 1);
                nu = 2 * cv/(1 + cv);
               while (abs(gb) > 1e-12){
                  gb = (lgamma(1 + 2 * nu) - 2 * lgamma(1 + 
                    nu) - aa)/(2 * (digamma(1 + 2 * nu) - digamma(1 + 
                    nu)));
                  nu = nu - gb;
                }
                shape = 1/nu;
                scale = exp(log(capmu) - lgamma(1 + nu));

  return shape;
  }
real weibullscale(real capmu, real capsig){
   real cv;
    real shape;
    real scale;
    cv = capsig/capmu;
    real aa;
    real gb;
    real nu;
              aa = log(cv^2 + 1);
                nu = 2 * cv/(1 + cv);
                while (abs(gb) > 1e-12){
                  gb = (lgamma(1 + 2 * nu) - 2 * lgamma(1 + 
                    nu) - aa)/(2 * (digamma(1 + 2 * nu) - digamma(1 + 
                    nu)));
                  nu = nu - gb;
                shape = 1/nu;
                scale = exp(log(capmu) - lgamma(1 + nu));

  return scale;
  }

I’m just trying to test this function, so the rest of the stan model is:

data {
  int<lower=0> N;
  vector[N] capmu;
  real capsig;
}
generated quantities {
  vector[N] wshape; 
  vector[N] wscale;
  for (i in 1:N){
    wshape[i] = weibullshape(capmu[i], capsig);
  }
  for (i in 1:N){
    wscale[i] = weibullscale(capmu[i], capsig);
  }
}

When I run this with algorithm="Fixed_param", I get the following:

d <- extract(mod, "wshape")$wshape[1,]
#[1] 0.6554046 0.6589957

d <- extract(mod, "wscale")$wscale[1,]
#0.4591410 0.4725465

Which isn’t really close to the values I get from R. I’m trying to determine whether:
a) the differences in calculations are due to an inherent difference in the lgamma and digamma implementations in R and Stan, or
b) there is a mismatch in the translation between the R function and the stan syntax. Based on my tests it appears that some part of the while loop is behaving differently (the values of gb produced by my stan function differ from those I calculate in R), or
c) there is a much better way to do this…

Maybe @Bob_Carpenter or @bgoodri has an idea?

Incidentally, I’d appreciate any good resources for folks just learning to write custom functions in Stan (i.e., efficient testing, etc)
Thanks!!
M

1 Like
2

Also, please let me know if this is an ill-formed question and I’m happy to try and add more details.

3

I don’t know if there is a better, more numerically stable way to do this calculation but I think the error is that you don’t initialize gb before using it in the loop condition: while (abs(gb) >= 1e-12). Try setting gb beforehand to a value > 1e-12, to make sure the loop executes at least once as in the R code.

1 Like
4

That fixed it!! Thanks for catching that - I didn’t realize that they were being initialized differently.