Below is a mwe with some suggested input-data. The under-/overflow isn’t an issue here as the cdf is not included. However it doesn’t work properly for the mwe either. Hopefully something is obvious to you. Thanks!

```
OS.event.times <- c(1,1,1,3,3,4,6,7,7,8)
OS.event <- c(0,0,0,0,1,0,1,0,1,0)
trial.data <- data.frame(OS.months = OS.event.times, OS.event = OS.event)
trial.data
```

#### Folder for STAN models

```
Stanfolder <- "./R/STAN/STAN_models/" # What folder to create for JAGS-model-file
if(!dir.exists(Stanfolder)) {dir.create(Stanfolder)} # checks if dir exists, creates it if not
```

# Setting up values for STAN

```
trial.N <- sum(trial.data$OS.event) # Number of datapoints
trial.events.N <- sum(trial.data$OS.event==1) # Number of events
trial.censored.N <- sum(trial.data$OS.event==0) # Number of censored
Trial.censored.vect <- rep(0, trial.censored.N) # Vector of censored events
trial.event.times <- trial.data$OS.months[trial.data$OS.event==1] # Trial event times
trial.censored.times <- trial.data$OS.months[trial.data$OS.event==0] # Trial censored times
```

```
# Stan code
cat("functions {
// shape b, scale eta
real gompertz_cdf (real y, real b, real eta) {
return 1 - exp(-eta*(exp(b*y) - 1));
}
real gompertz_lpdf (real y, real b, real eta) {
return log(b*eta) + eta + b*y - eta*exp(b*y);
}
real gompertz_lcdf (real y, real b, real eta) {
return log1m_exp(-eta*expm1(b*y));
}
real gompertz_lccdf (real y, real b, real eta){
return -eta*(exp(b*y) - 1);
}
}
data {
int N;
int<lower=0> Nobs;
int<lower=0> Ncen;
vector[Nobs] yobs;
vector[Ncen] ycen;
}
parameters {
real <lower=0, upper=0.1> b; // param 1
real <lower=0, upper=0.1> eta; // param 2
}
model {
// priors
b ~ gamma(1, 20);
eta ~ gamma(1, 20);
for (i in 1:Nobs) {
target += gompertz_lpdf(yobs[i] | b, eta);
}
for (j in 1:Ncen) {
target += gompertz_lccdf(ycen[j] | b, eta);
}
}
```

```
", file=paste0("./R/STAN/STAN_models/Stan_Model_gompertz.stan"))
fit <- stan("./R/STAN/STAN_models/Stan_Model_gompertz.stan", data=list(N = trial.N, Nobs = trial.events.N,
Ncen = trial.censored.N, yobs = trial.event.times, ycen = trial.censored.times,
trial_cen_vector = Trial.censored.vect), iter=5000, chains=4)
fit
```

[edit: escaped code—it requires back-ticks, not apostrophes (```)]