# Bayesian Survival Function
**URL:** https://discourse.pymc.io/t/bayesian-survival-function/3126
**Category:** Questions
**Created:** [April 18, 2019, 1:25am UTC](https://discourse.pymc.io/t/bayesian-survival-function/3126 "2019-04-18T01:25:03Z")
**Posts on this page:** 5
**Page:** 1
### Author:  [@sachinruk](https://discourse.pymc.io/u/sachinruk)
#### Post date: [April 18, 2019, 1:25am UTC](https://discourse.pymc.io/t/bayesian-survival-function/3126/1 "2019-04-18T01:25:03Z")
This is a crosspost: [https://stats.stackexchange.com/questions/403714/bayesian-survival-function](https://stats.stackexchange.com/questions/403714/bayesian-survival-function)
I have managed to estimate the posterior of the latent variables of my model which can be stated as follows (adapted from [https://docs.pymc.io/notebooks/bayes\_param\_survival\_pymc3.html](https://docs.pymc.io/notebooks/bayes_param_survival_pymc3.html)):
y\_i = \beta^T x\_i + \xi\_i
where y\_i\equiv \log T\_i, is the log duration and,
\xi\_i \sim Gumbel(0, s) \\ s \sim HalfNormal(5)
Note that this can be equivalently stated that y\_i\sim Gumbel(\beta^T x\_i, s).
Therefore, the likelihood of the dataset for both censored and uncensored data is:
\log p\left(\cup\_{i=1}^N y\_i|\beta, s\right) = \sum\_{i \in uncensored} \log p\left(y\_i=\log T\_i|\beta, s\right) + \sum\_{i \in censored} \log p\left(y\_i\>\log T\_i|\beta, s\right)
**The question is regarding getting a Bayesian estimate of the survival function**. The author of the above notebook, samples from the posterior predictive distribution of y using:
```auto
X_.set_value(np.array(x_test)[None,:])
cens_.set_value(np.array([False]))
Y_.set_value(np.zeros(1))
with weibull_model:
post_trace = pm.sample_posterior_predictive(trace, vars=[y_obs])
```
and uses the samples of y (converted back to time domain) to generate the survival function with the following logic (for time steps 1 to 200):
```auto
for t in range(1, 201):
frac = sum(t_samples > t) / len(t_samples)
survival_func[t] = frac
```
It makes logical sense to state that the fraction of samples greater than a particular time is the survival rate. However, this is a point estimate of the survival function.
There’s two ways I can think of getting a ‘Bayesian’ estimate of the Survival Function:
1. Sample N number of traces (with say 1000 points) of the posterior predictive distribution and do:
```auto
for n in range(N):
for t in range(1, 201):
frac = sum(t_samples[n] > t) / len(t_samples[n])
survival_func[N, t] = frac
```
1. We know the CDF of the gumbel distribution _given the latent variables_ to be, \exp\left(-\exp\left(\left(\log T\_i - X\_i^T\beta\right)/s\right)\right). Hence we can let the j-th sample for the survival function for a given value of T to be 1 - \exp\left(-\exp\left(\left(\log T\_i - X\_i^T\beta^{(j)}\right)/s^{(j)}\right)\right) where \beta^{(j)}, s^{(j)} are the samples from the posterior of the latent variables.
Are the two methods equivalent? If not which method is the more correct way of getting a Bayesian sample of the Survival Function.
## Edit
I feel like I should point out after @Gon_F’s answer that in method 2. What I am calculating is the expectation of \exp\left(-\exp\left(\left(\log T - X^T\beta\right)/s\right)\right) under the distribution p(\beta, s| \mathcal{D}\_{train}) where \mathcal{D} is the training dataset. Whereas, the first method is looking directly at the distribution, p(T| \mathcal{D}\_{train}) where, \beta, s has **already been marginalised out**. So I’m wondering if (1) is an approximation to (2) after all. Hopefully I’m not complicating things more.
---
### Author:  [@chartl](https://discourse.pymc.io/u/chartl)
#### Post date: [April 18, 2019, 1:51pm UTC](https://discourse.pymc.io/t/bayesian-survival-function/3126/2 "2019-04-18T13:51:25Z")
I love this question.
The way it is written, (1) is a Monte Carlo approximation to (2). In order to sample from the density in (2) you would:
(i) Sample \beta^{(j)}, s^{(j)} from P\_\mathrm{post}(\beta, s)
(ii) Sample y\_i^{(j)} from \mathrm{Gum}(x\_i^T\beta^{(j)}, s)
(iii) Take T\_i^{(j)} = \exp(y\_i^{(j)})
This is precisely the content of the trace for y\_i.
---
### Author:  [@Gon\_F](https://discourse.pymc.io/u/Gon_F)
#### Post date: [April 18, 2019, 6:30pm UTC](https://discourse.pymc.io/t/bayesian-survival-function/3126/3 "2019-04-18T18:30:46Z")
If we do find an effective way to do Survival functions, it would open a path to implementing joint models, like [Stan](http://www.hammerlab.org/2017/06/26/introducing-survivalstan/) has, which allow modelling evolving Time-to-Event’s (TTE’s) over time, which I had been looking at recently.
---
### Author:  [@lucianopaz](https://discourse.pymc.io/u/lucianopaz)
#### Post date: [April 20, 2019, 8:07pm UTC](https://discourse.pymc.io/t/bayesian-survival-function/3126/4 "2019-04-20T20:07:28Z")
Maybe @AustinRochford, can comment on this
---
### Author:  [@chartl](https://discourse.pymc.io/u/chartl)
#### Post date: [April 23, 2019, 10:31pm UTC](https://discourse.pymc.io/t/bayesian-survival-function/3126/5 "2019-04-23T22:31:22Z")
> [@sachinruk](#):
>
> What I am calculating is the expectation of \exp(−\exp((\mathrm{log}T−X^TTβ)/s))\exp\left(-\exp\left(\left(\log T - X^T\beta\right)/s\right)\right) under the distribution p(\beta,s|\mathcal{D}\_{train}) where \mathcal{D} is the training dataset. Whereas, the first method is looking directly at the distribution, p(T|\mathcal{D}\_{train})
To be clear set \mathcal{D}\_{train} = (X, y); and note that y only influences the predictive posterior p(T|X, y) through its influence on the posteriors p(\beta, x | X, y):
p(T|X, y) = \int\_{\beta, s} p(T|\beta, s, X)p(\beta, s|X, y)d\beta ds
or:
p(T|D) = \int\_{\beta ,s} p(T|\beta, s, D)p(\beta, s|D)d\beta ds
So sampling from p(\beta, s|D) and then from p(T|\beta, s) and computing computing \frac{1}{K}\sum\_{i=1}^kT\_i or \frac{\sum\_{i=1}^K T\_i P(T\_i|\beta\_i,s\_iX)p(\beta\_i, s\_i | D)}{\sum p(\beta\_i, s\_i|D)} is a monte-carlo approximation to the (numerical) expectation \int T p(T|D) d T