# Simple Exponential Survival Function

**URL:** <https://discourse.pymc.io/t/simple-exponential-survival-function/4961>\
**Category:** Questions\
**Created:** [April 28, 2020, 3:50pm UTC](https://discourse.pymc.io/t/simple-exponential-survival-function/4961 "2020-04-28T15:50:53Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![john\_c](https://avatars.discourse-cdn.com/v4/letter/j/ee59a6/32.png) [@john\_c](https://discourse.pymc.io/u/john_c)\
**Post date:** [April 28, 2020, 3:50pm UTC](https://discourse.pymc.io/t/simple-exponential-survival-function/4961/1 "2020-04-28T15:50:53Z")

</div>

I’m having trouble replicating survival analysis in PyMC3. I’ve been following along with the Statistical Rethinking lectures, and there’s an example about survival analysis on cat adoptions that isn’t in the published book. I’m wondering how to convert it into a PyMC3 model, particularly this sort of logic:

![Screen Shot 2020-04-28 at 11.46.29 AM](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/2X/6/6038181d6601bc42634e9b098cefeba299451dac.png)

If there is an instance of adoption it isn’t censored data so I can use a regular exponential distribution, and if it isn’t labeled as an adoption it is censored data because the cat might still end up being adopted so I’d use the exponential-CCDF. How do I implement this sort of logic in PyMC3?

Here’s stan code for the model as well:

```auto
model{
vector[22356] lambda; 
b ~ normal( 0 , 1 ); 
a ~ normal( 0 , 1 ); 
for ( i in 1:22356 ) {
        lambda[i] = a + b * black[i];
        lambda[i] = exp(lambda[i]);
    }

for ( i in 1:22356 )
    if ( adopted[i] == 0 ) target += exponential_lccdf(days_to_event[i] | lambda[i]);
for ( i in 1:22356 )
    if ( adopted[i] == 1 ) days_to_event[i] ~ exponential( lambda[i] );
}

```

---

<div class="post-metadata">

**Author:** ![junpenglao](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/junpenglao/32/8_2.png) [@junpenglao](https://discourse.pymc.io/u/junpenglao)\
**Post date:** [April 29, 2020, 7:43am UTC](https://discourse.pymc.io/t/simple-exponential-survival-function/4961/2 "2020-04-29T07:43:36Z")

</div>

You can use a potential for that. Note that we have logcdf implemented for `pm.Exponential` but not the survival function (i.e., lccdf), so the best way I can think of is:

```auto
with pm.Model():
    a = pm.Normal('a', 0, 1)
    b = pm.Normal('b', 0, 1)
    log_rate = a + b * black
    y = pm.Exponential('observed', tt.exp(log_rate), observed=days_to_event)
    pm.Potential('survival', -log_rate[adopted == 0]) ==> exp_lccdf = exp_log_prob(lambda, rate) - log_rate

```

---

<div class="post-metadata">

**Author:** ![john\_c](https://avatars.discourse-cdn.com/v4/letter/j/ee59a6/32.png) [@john\_c](https://discourse.pymc.io/u/john_c)\
**Post date:** [April 29, 2020, 1:51pm UTC](https://discourse.pymc.io/t/simple-exponential-survival-function/4961/3 "2020-04-29T13:51:39Z")

</div>

> [@junpenglao](#):
>
> ```auto
> pm.Potential('survival', -log_rate[adopted == 0]) ==> exp_lccdf = exp_log_prob(lambda, rate) - log_rate
> 
> ```

I couldn’t end up following your solution as well as I had hoped, part of it is being new to PyMC3 and part of it is needing to work on my theoretical knowledge more. I did however find this exponential survival function from PyMC3s website (which I just managed to derive myself after and figure out how they got it) ![Screen Shot 2020-04-29 at 9.37.46 AM](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/2X/a/a8b5ec905cf59f8668cc8b983555e43757b1c6f3.png)  
and a corresponding function

```
def exponential_logp(censored, value):
    return (censored * tt.log(lam) - value*lam).sum()

```

So the model ends up being

```
with pm.Model() as m1:
    a_c = pm.Normal('a_color', 0, 1, shape=len(df.color_id.unique()))
    mu = tt.exp(a_c[cid])
    lam = 1/mu
    exp_surv = pm.DensityDist('exp_surv', exponential_logp, observed={'censored':censored, 'value':days_to_event})
    
    m1_trace = pm.sample(draws=1000)

```

And this model is returning the exact same inference as the Stan model. However, when I try to to run  
`pm.sample_posterior_predictive(trace_m1)`  
I get the following error:  
`ValueError: Distribution was not passed any random method Define a custom random method and pass it as kwarg random`

Are there any good resources on how to fix this issue and create a custom distribution from my exponential\_logp function? I attempted to troubleshoot the issue but I haven’t been getting very far

edit: with a little more effort I ended up following along with your code, which unlike with my approach, allows for sampling. Thanks!

---

<div class="post-metadata">

**Author:** ![roesta07](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/roesta07/32/6060_2.png) [@roesta07](https://discourse.pymc.io/u/roesta07)\
**Post date:** [June 18, 2021, 1:04pm UTC](https://discourse.pymc.io/t/simple-exponential-survival-function/4961/4 "2021-06-18T13:04:10Z")

</div>

@john_c I have been trying to solve the same problem. Did you end up using Potential? If so, How did you incorporated your model? Seems like you are missing priors for exponential\_logp !! Thanks

---

<div class="post-metadata">

**Author:** ![broesler](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/broesler/32/5953_2.png) [@broesler](https://discourse.pymc.io/u/broesler)\
**Post date:** [January 5, 2024, 10:32pm UTC](https://discourse.pymc.io/t/simple-exponential-survival-function/4961/5 "2024-01-05T22:32:14Z")

</div>

The error you’re getting with `DensityDist` is telling you that the model does not know how to generate a random sample from your custom distribution (see [Custom Distribution](https://www.pymc.io/projects/docs/en/latest/api/distributions/generated/pymc.CustomDist.html)).

You need to define a sampling function `random`, then pass `DensityDist(..., random=random)`. In the cats example, the function would be:

```auto
from scipy import stats

def random(*dist_params, rng=None, size=None):
    """The function from which to generate random samples.

    Parameters
    ----------
    dist_params : list
        The same parameters as given to `exponential_logp`. (t, lam) in this case.
    rng : Generator
        The random number generator used internally.
    size : tuple
        The desired size of the random draw.

    Returns
    -------
    result : (N,) ndarray
        A random sample of `size` from the desired distribution.
    """
    t, lam = dist_params
    return np.where(
        np.array(adopted == 1),
        stats.expon.rvs(scale=1/lam, size=size), # not censored = Exp(lam)
        stats.expon.isf(stats.uniform.rvs(size=size), scale=1/lam) # censored = Exp(lam).isf()
    )

# NOTE DensityDist -> CustomDist in new API
exp_surv = pm.CustomDist('exp_surv', days_to_event, lam,
                         logp=logp, observed=adopted,
                         random=random)

```

This addition will allow you to sample from the posterior predictive. Note that I’ve only done a basic check on the mean of the posterior predictive, so I’m not sure if this function is entirely accurate.

@john_c Can you post your code for using `pm.Potential` to solve the issue? I cannot seem to get it to work properly either.

---

<div class="post-metadata">

**Author:** ![broesler](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/broesler/32/5953_2.png) [@broesler](https://discourse.pymc.io/u/broesler)\
**Post date:** [January 5, 2024, 11:01pm UTC](https://discourse.pymc.io/t/simple-exponential-survival-function/4961/6 "2024-01-05T23:01:43Z")

</div>

Update: I figured out the `pm.Potential` solution. In @junpenglao 's terms, the line is

```auto
pm.Potential('survival', pm.math.switch(adopted, 0, -log_rate))

```

In @john_c 's terms, the line is

```auto
pm.Potential('survival', pm.math.switch(adopted, 0, a_c[cid]))

```

For reasons discussed in McElreath’s book, use of the index variable `cid` in the latter is preferred over the indicator variable `black` from the former.

---

<div class="post-metadata">

**Author:** ![broesler](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/broesler/32/5953_2.png) [@broesler](https://discourse.pymc.io/u/broesler)\
**Post date:** [January 5, 2024, 11:21pm UTC](https://discourse.pymc.io/t/simple-exponential-survival-function/4961/7 "2024-01-05T23:21:32Z")

</div>

I get the following error when I try to run `pm.sample_posterior_predictive` on your model with `pm.Potential`:

```auto
UserWarning: The effect of Potentials on other parameters is ignored during posterior predictive sampling. This is likely to lead to invalid or biased predictive samples.

```

Is there any further information on this error? Since the Potential is only used to estimate the main parameter in this case, does this warning even matter?

---

<div class="post-metadata">

**Author:** ![ricardoV94](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ricardov94/32/5775_2.png) [@ricardoV94](https://discourse.pymc.io/u/ricardoV94)\
**Post date:** [January 6, 2024, 8:45am UTC](https://discourse.pymc.io/t/simple-exponential-survival-function/4961/8 "2024-01-06T08:45:00Z")

</div>

If you’re using a Potential instead of CystomDist for the likelihood PyMC won’t know how to perform posterior predictive sampling of censored variables (hence the warning). If you don’t care about censoring in posterior predictive you can ignore the warning.

Note that we now have a Censored class that can be used in many of these survival models [Censored — PyMC dev documentation](https://www.pymc.io/projects/docs/en/latest/api/distributions/censored.html)

Which alleviates the need for CustomDist or Potentials. It doesn’t work based on indicator variables but based on specifying the censoring bounds. If you specify `Censored(..., upper=[5, np.inf], observed=[5, 6])` would correspond to a case where the first observation had a value of 5 and was censored and the second had a value of 6 that was not censored.

---

<div class="post-metadata">

**Author:** ![swpease](https://avatars.discourse-cdn.com/v4/letter/s/6de8d8/32.png) [@swpease](https://discourse.pymc.io/u/swpease)\
**Post date:** [October 21, 2024, 11:16pm UTC](https://discourse.pymc.io/t/simple-exponential-survival-function/4961/9 "2024-10-21T23:16:34Z")

</div>

In case anyone else goes the `Censored` class route, blindly copies over the example, and gets an error about “non-finite logp”, that’s just due to the prior on the `a`’s being Normal(0,1). `a` is supposed to be the log average wait time in days, so setting it to something more reasonable, like just Normal(2,1), solves the issue. The wait times of ~1000 days do not compute (ha!) with that N(0,1) prior.