# Plotting survival curve with a Bayesian inference?

**URL:** <https://discourse.pymc.io/t/plotting-survival-curve-with-a-bayesian-inference/12770>\
**Category:** v5\
**Created:** [August 21, 2023, 12:59pm UTC](https://discourse.pymc.io/t/plotting-survival-curve-with-a-bayesian-inference/12770 "2023-08-21T12:59:43Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![Vincent](https://avatars.discourse-cdn.com/v4/letter/v/e9bcb4/32.png) [@Vincent](https://discourse.pymc.io/u/Vincent)\
**Post date:** [August 21, 2023, 12:59pm UTC](https://discourse.pymc.io/t/plotting-survival-curve-with-a-bayesian-inference/12770/1 "2023-08-21T12:59:43Z")

</div>

```python
import pymc as pm
import numpy as np
import pytensor
import pytensor.tensor as pt
import arviz as az

# Set values for the prior parameters
gamma_mode = 0.5
gamma_variance = 1.0
t_star_mode = 15.0
t_star_variance = 100.0
s_star = 0.88
lambda_ = -pt.log(s_star)

# Definition of prior for gamma
aux1 = 2.0 + (gamma_mode**2.0) / gamma_variance
a = 0.5 * (aux1 + np.sqrt(aux1**2.0 - 4.0))
b = gamma_mode / (a - 1.0)
b1 = 1.0 / b

# Definition of prior for t_star
aux2 = 2.0 + (t_star_mode**2.0) / t_star_variance
ap = 0.5 * (aux2 + np.sqrt(aux2**2.0 - 4.0))
bp = t_star_mode / (ap - 1.0)
b2 = 1.0 / bp

# Define the model
model = pm.Model()

# Observed data (censored)
t = np.array([25.62191781, 25.35616438, 25.34246575, 25.29589041, 24.61643836, 24.34794521, 24.30410959,
      24.06027397, 23.81643836, 23.50410959, 23.2109589, 23.14520548, 22.97808219, 22.70410959,
      21.90684932, 21.81917808, 21.36986301, 20, 19.5260274])

# Function to calculate log complementary cdf of Weibull distribution
def weibull_lccdf(x, t_star, gamma):
    return -(lambda_ * (x / t_star) ** gamma)

with model:
    # Definition of priors for gamma and t_star
    gamma = pm.Gamma('gamma', alpha=a, beta=1/b)
    t_star = pm.Gamma('t_star', alpha=ap, beta=1/bp)

    # Bus survivor calculation
    surv = pm.Deterministic('surv', pt.exp(-lambda_ * (pt.arange(0, 26) / t_star) ** gamma))

    # Define likelihood
    y_cens = pm.Potential("y_cens", weibull_lccdf(t, t_star, gamma))

    # Sample from the model
    trace = pm.sample(draws=10000, tune=1000)

summary = pm.summary(trace)
print(summary)

```

I wanted to plot the survival fonction “surv” in my code but it is not working. Does anyone can help on this task ?

---

<div class="post-metadata">

**Author:** ![cluhmann](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/cluhmann/32/3083_2.png) [@cluhmann](https://discourse.pymc.io/u/cluhmann)\
**Post date:** [August 21, 2023, 4:02pm UTC](https://discourse.pymc.io/t/plotting-survival-curve-with-a-bayesian-inference/12770/2 "2023-08-21T16:02:31Z")

</div>

I tried fixing the formatting of your code, but please let me know if I didn’t get it quite right.

I would suggest checking out the [survival analysis notebooks](https://www.pymc.io/projects/examples/en/latest/gallery.html#survival-analysis) for examples of who to visualize what is going on. If you have tried plotting and it is not working, you can share that code as well.

---

<div class="post-metadata">

**Author:** ![Vincent](https://avatars.discourse-cdn.com/v4/letter/v/e9bcb4/32.png) [@Vincent](https://discourse.pymc.io/u/Vincent)\
**Post date:** [August 22, 2023, 8:09am UTC](https://discourse.pymc.io/t/plotting-survival-curve-with-a-bayesian-inference/12770/3 "2023-08-22T08:09:52Z")

</div>

```python
# Assuming you have already calculated the "surv" function from your first code
# Replace the following lines with your calculated "surv" function and interval bounds
# For example, assuming you have `surv` and `interval_bounds` from your first code
surv_function = surv # Replace with your actual "surv" function
interval_bounds = np.arange(0, 26) # Replace with your actual interval bounds

# Plot survival function
fig, surv_ax = plt.subplots(figsize=(8, 6))

# Plot credible intervals of survival function for your data
az.plot_hdi(interval_bounds[:-1], surv_function, ax=surv_ax, smooth=False, color="C0")

# Plot mean survival function
surv_ax.plot(interval_bounds[:-1], get_mean(surv_function), color="darkblue")

# Customize the plot for survival function
surv_ax.set_xlim(0, interval_bounds.max())
surv_ax.set_xlabel("Months since mastectomy")
surv_ax.set_ylabel("Survival function $S(t)$")

# Add a title to the figure
fig.suptitle("Survival Function Visualization")

```

Thanks, so I add the following code above but it still didn’t work:

```python
UFuncTypeError Traceback (most recent call last)
<ipython-input-9-d5ec0b9e0bd3> in <module>
     68 
     69 # Plot credible intervals of survival function for your data
---> 70 az.plot_hdi(interval_bounds[:-1], surv_function, ax=surv_ax, smooth=False, color="C0")
     71 
     72 # Plot mean survival function

~\AppData\Roaming\Python\Python38\site-packages\arviz\plots\hdiplot.py in plot_hdi(x, y, hdi_prob, hdi_data, color, circular, smooth, smooth_kwargs, figsize, fill_kwargs, plot_kwargs, hdi_kwargs, ax, backend, backend_kwargs, show)
    154 elif not 1 >= hdi_prob > 0:
    155 raise ValueError("The value of hdi_prob should be in the interval (0, 1]")
--> 156 hdi_data = hdi(y, hdi_prob=hdi_prob, circular=circular, multimodal=False, **hdi_kwargs)
    157 
    158 hdi_shape = hdi_data.shape

~\AppData\Roaming\Python\Python38\site-packages\arviz\stats\stats.py in hdi(ary, hdi_prob, circular, multimodal, skipna, group, var_names, filter_vars, coords, max_modes, dask_kwargs, **kwargs)
    588 if isarray and ary.ndim <= 1:
    589 func_kwargs.pop("out_shape")
--> 590 hdi_data = func(ary, **func_kwargs) # pylint: disable=unexpected-keyword-arg
    591 return hdi_data[~np.isnan(hdi_data).all(axis=1), :] if multimodal else hdi_data
    592 

~\AppData\Roaming\Python\Python38\site-packages\arviz\stats\stats.py in _hdi(ary, hdi_prob, circular, skipna)
    631 interval_idx_inc = int(np.floor(hdi_prob * n))
    632 n_intervals = n - interval_idx_inc
--> 633 interval_width = np.subtract(ary[interval_idx_inc:], ary[:n_intervals], dtype=np.float_)
    634 
    635 if len(interval_width) == 0:

UFuncTypeError: Cannot cast ufunc 'subtract' input 0 from dtype('O') to dtype('float64') with casting rule 'same_kind'

```