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