Fitting Weibull CDF to data in rstan

Modeling
1

I have a rather complicated model, however I think the main issue is stemming from the fit of a Weibull CDF to data. My issue can be reproduced in this way:
generate data in R:
y ← pweibull(1:24, shape=2, scale=8)
N ← length(y)

The Stan model in a .stan file is :
data {
int<lower=0> N;
vector[N] y;
}

parameters {
real<lower=0> shape;
real<lower=0> scale;
}

model {
shape ~ normal(3, 1);
scale ~ normal(7, 2);
target += weibull_lcdf(y | shape, scale);
}

The problem is this: After running stan from R using the above simulated data, I get very strange results. the parameter scale is set to 0.02 with a sd of 0.02, so statistically zero, and the shape is 2.87 with sd of 1.07
What is wrong? Why I cannot fit such an easy model? I am new to Stan, so any help is appreciated.

2

Welcome to the community. Is there a reason why you are using the CDF rather than the PDF? I would try passing 1:24 as y and using weibull_lpdfrather than weibull_lcdf.

3

Is there a particular reason why you cannot model y directly? If you have access to the data, simply doing:

y <- rweibull(24, shape=2, scale=8)
N <- length(y)

data {
int<lower=0> N;
vector[N] y;
}

parameters {
real<lower=0> shape;
real<lower=0> scale;
}

model {
shape ~ normal(3, 1);
scale ~ normal(7, 2);
y ~ weibull(shape, scale);
}

will work just fine (note y ~ weibull(shape, scale); is more or less equivalent to target +-= weibull_lpdf)

1 Like
4

Thank you for responding so quickly.
Maybe this can explain my problem better.
I have measured the percentage of a drug absorbed at time t. So, for example I have 0.1 at t=1, 0.3 at t=2 and so on. The plot of the absorption over time looks like a Weibull CDF. My observations are the percentages over time (1:24).
How can I estimate the parameters of that Weibull CDF knowing these percentages and incorporating the times?

1 Like
5

Ah, I see more clearly what you are trying to do now. A generative model for you situation might be the concentration evolves over time with a deterministic mean trend defined by the Weibull CDF curve plus some observation noise.

t <- 1:24
N <- length(t)
shape <- 2
scale <- 8
sigma <- 0.01
y <- pweibull(t, shape, scale) + rnorm(N, 0 , sigma)

data {
int<lower=0> N;
vector[N] y;
vector[N] t;
}

parameters {
real<lower=0> shape;
real<lower=0> scale;
real<lower=0> sigma;
}

model {
shape ~ normal(3, 1);
scale ~ normal(7, 2);
sigma ~ normal(0,1);
vector[N] f;

for (i in 1:N){
  f[i] = weibull_cdf(t[i] | shape, scale); //webiull curve as deterministic trend
}

y ~ normal(f, sigma); //observational model

}

Of course, this will allow for possible values less than 0 or above 1, but it might be a good start. For refinement, perhaps consider a beta observational model (or a similar distribution with values constrained to be between 0 and 1).

2 Likes