# Help with simple survival data simulation and exponential model fitting

**URL:** <https://discourse.mc-stan.org/t/help-with-simple-survival-data-simulation-and-exponential-model-fitting/40193>\
**Category:** Modeling\
**Tags:** fitting-issues\
**Created:** [August 18, 2025, 7:10am UTC](https://discourse.mc-stan.org/t/help-with-simple-survival-data-simulation-and-exponential-model-fitting/40193 "2025-08-18T07:10:26Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![ruben-s](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/ruben-s/32/20878_2.png) [@ruben-s](https://discourse.mc-stan.org/u/ruben-s)\
**Post date:** [August 18, 2025, 7:10am UTC](https://discourse.mc-stan.org/t/help-with-simple-survival-data-simulation-and-exponential-model-fitting/40193/1 "2025-08-18T07:10:26Z")

</div>

I’m doing a really basic hazard model fitting using simulated data. Here is the Stan model:

```stan
data {
    int N_subjects;
    vector[N_subjects] time;
    vector[N_subjects] event;
    vector[N_subjects] group;
}
parameters {
    real<lower=0> h_0; // baseline mean survival time
    real B_h1;
}
transformed parameters {
    vector[N_subjects] lambda;

    lambda = (1.0 / h_0) * exp(group * B_h1);
}
model {
    h_0 ~ lognormal(log(180), 0.2);
    B_h1 ~ normal(0, 10);

    target += event .* exponential_lpdf(time | lambda);
    target += (1 - event) .* exponential_lccdf(time | lambda);
}

```

I’m the simulating some data using the following code (this is LLM-generated code but seems reasonable enough for the simulation data):

```no-highlight
N_SUBJECTS = 2000
H0_TRUE = 180.0 # True baseline mean survival time
B_H1_TRUE = -0.6 # The true coefficient we want to recover
CENSOR_RATE = 0.3 # Percentage of subjects to be censored

# --- 2. Simulate Data ---
np.random.seed(42)

# Create the covariate (group)
group = np.random.randint(0, 2, size=N_SUBJECTS)

# Calculate the rate for each subject based on the PH model
# rate = (1 / baseline_mean_time) * hazard_ratio
lambda_ind = (1.0 / H0_TRUE) * np.exp(group * B_H1_TRUE)

# Generate true event times from an exponential distribution
# T = -log(U) / lambda, where U ~ Uniform(0,1)
u = np.random.uniform(0, 1, size=N_SUBJECTS)
true_event_time = -np.log(u) / lambda_ind

# Generate censoring times from a separate exponential distribution
# This creates random right-censoring
censor_time = np.random.exponential(scale=H0_TRUE / CENSOR_RATE, size=N_SUBJECTS)

# Observe the final time and event status
time = np.minimum(true_event_time, censor_time)
event = (true_event_time <= censor_time).astype(int)

# Create a DataFrame
df = pd.DataFrame({"time": time, "event": event, "group": group})

```

For some reason, I’m able to recover the `B_h1` relative risk more closely using sksurv’s `CoxPHSurvivalAnalysis`, while the `B_h1` estimate produced by Stan is the following:

```no-highlight
            Mean MCSE ... ESS_tail R_hat
h_0 200.653000 0.004460 ... 1536.30 1.00487
B_h1 -0.377075 0.000034 ... 1569.86 1.00420

```

(For reference, the `CoxPHSurvivalAnalysis` model produces an estimate of `-0.54`)

I’m completely stumped here - what could possibly be going wrong? I’m new to Stan and just trying to figure out how to get through the basics, but somehow got caught up with a very simple model.

Thank you for the help!

---

<div class="post-metadata">

**Author:** ![amas0](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/amas0/32/9754_2.png) [@amas0](https://discourse.mc-stan.org/u/amas0)\
**Post date:** [August 18, 2025, 3:49pm UTC](https://discourse.mc-stan.org/t/help-with-simple-survival-data-simulation-and-exponential-model-fitting/40193/2 "2025-08-18T15:49:14Z")

</div>

Hi @ruben-s , welcome to the Stan community!

There is an error with how the `target += ` statements are written, where you aren’t actually fitting the model you intended to fit. If you write it more carefully in an iterative, you should get the results you are looking for:

```stan
for (n in 1:N_subjects) {
    if (event[n] == 1)
        target += exponential_lpdf(time[n] | lambda[n]);
    else
        target += exponential_lccdf(time[n] | lambda[n]);
}

```

When I run this updated model using your data generation script, I get a posterior mean for `B_h1` of `-0.561` (and a 90% interval of `(-0.65, -0.47)`).

The core disconnect here lies in how Stan’s vectorization works. When you write `exponential_lpdf(time | lambda)` where `time` and `lambda` are vectors, this does not produce a vector of log probabilities, it produces the sum of those log probabilities. So, when you write `event * exponential_lpdf(time | lambda)`, this produces a vector that is 0 wherever event is 0 and \sum\_n \log p(t\_n | \lambda\_n) wherever event is 1.

(As a side note; the Stan code you provided is not valid, `event .* exponential_lpdf(time | lambda)` should not compile with the left-side being a vector and the right side being a real number. I’m assuming you meant `*`, but it’s unclear.)

Some other thoughts on the code provided. You code all your input data as vectors, which can be fine, but given that `event` and `group` are binary flags, it may be more reasonable to code them as `array[N_subjects] int<lower=0,upper=1>`. This will perform data validation on inputs, which can avoid some pernicious bugs.

Also, you construct `lambda` within the `transformed parameters` block. Since there are only two values for `group`, you are performing unnecessary calculations (and returning the `N_subjects` length array as outputs). Instead, you can drop the `transformed parameter` block and instead compute `lambda` within the `model` block like:

```stan
for (n in 1:N_subjects) {
    real lambda = (1.0 / h_0) * exp(group[n] * B_h1);
    if (event[n] == 1)
        target += exponential_lpdf(time[n] | lambda);
    else
        target += exponential_lccdf(time[n] | lambda);
}

```

---

<div class="post-metadata">

**Author:** ![ruben-s](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/ruben-s/32/20878_2.png) [@ruben-s](https://discourse.mc-stan.org/u/ruben-s)\
**Post date:** [August 18, 2025, 4:31pm UTC](https://discourse.mc-stan.org/t/help-with-simple-survival-data-simulation-and-exponential-model-fitting/40193/3 "2025-08-18T16:31:46Z")

</div>

Wow, thank you so much for the awesome explanation. I really appreciate it! I can’t believe I got so much wrong in such a small amount of code…haha