# Recreating convoys in pymc

**URL:** https://discourse.pymc.io/t/recreating-convoys-in-pymc/7923
**Category:** Questions
**Created:** [August 22, 2021, 5:49pm UTC](https://discourse.pymc.io/t/recreating-convoys-in-pymc/7923 "2021-08-22T17:49:21Z")
**Posts on this page:** 2
**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: [August 22, 2021, 5:49pm UTC](https://discourse.pymc.io/t/recreating-convoys-in-pymc/7923/1 "2021-08-22T17:49:21Z")

</div>

I’m interested in recreating the algorithm from [convoys](https://github.com/better/convoys/) for generalized gamma survival analysis (that converges to some conversion rate, c). The loss function is posted [here](https://github.com/better/convoys/blob/99e832d2170ba9670e5e16bb3f632ac9055291f3/convoys/regression.py#L75)

Does anyone have any pointers on how to get started for this?  
I understand how the `generalized_gamma_loss` function calculcates the log probability, but I have no idea how to make this compatible with pymc3

---

<div class="post-metadata">

### Author: ![KyleJCaron](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/kylejcaron/32/2715_2.png) [@KyleJCaron](https://discourse.pymc.io/u/KyleJCaron)
#### Post date: [August 23, 2021, 4:12pm UTC](https://discourse.pymc.io/t/recreating-convoys-in-pymc/7923/2 "2021-08-23T16:12:00Z")

</div>

I found an [old PR from the convoys package](https://github.com/better/convoys/blob/b40c25e8f5ecd43f49c61640ec0c26fb50464015/convoys/bayesian.py) that actually attempted using a Weibull distribution with pymc3 and adapted it slightly to have a parameterization that more closely mimics pymc3’s. Here’s the code. It shouldnt be too difficult to adapt this to the generalized gamma

And if anyone has suggestions for priors to reduce some of the divergences and improve regularization, I’d be very happy to hear them.

```auto
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import scipy.stats as stats
import pymc3 as pm
import arviz as az
from scipy.special import expit
from pymc3.math import dot, sigmoid, log, exp
logistic = lambda x: np.exp(x)/ (1+np.exp(x))

N = 10000
days = 100
t = np.arange(0, days,1)
T_true = pm.Weibull.dist(2.5, 70).random(size=N) # true time to event (independent of conversion itself)

X = np.random.normal(size=N)[:,None]
b = -0.2
CVR = logistic( b*X + -0.5 )

# only call it a conversion if they wouldve converted by now AND they convert
p = np.random.binomial(1, CVR)
B = ((T_true < t.max()) & (p.ravel() == 1))*1

# adjust T for non-converts
ttc = (t >= np.where(B==1, T_true, 10000)[:,None]*1)
T = np.where((p.ravel()==0) | (T_true > t.max()), t.max(),T_true)

plt.plot(ttc.mean(axis=0))
plt.title("Time to Conversion ")
plt.ylabel("CVR")
plt.xlabel("Time")

n, ncoefs = X.shape

with pm.Model() as m:
    alpha_c = pm.Normal('alpha_c', -5, 1) # intercept for conversion rate
    
    beta_sd = pm.Exponential('beta_sd', 1.0) # Weak prior for the regression coefficients
    beta = pm.Normal('beta', mu=0, sd=beta_sd, shape=(ncoefs,)) # Regression coefficients
    
    c = sigmoid(dot(X, beta) + alpha_c) # Conversion rates for each example
    k = pm.Lognormal('k', mu=0.5, sd=1.0) # Weak prior around k=1
    lambd = pm.Exponential('lambd', 0.01) # Weak prior

    # PDF of Weibull: k / lambda * (x / lambda)^(k-1) * exp(-(t / lambda)^k)
    LL_observed = log(c) + log(k) - log(lambd) + (k-1)*(log(T) - log(lambd)) - (T/lambd)**k
    # CDF of Weibull: 1 - exp(-(t / lambda)^k)
    LL_censored = log((1-c) + c * exp(-(T/lambd)**k))

    # We need to implement the likelihood using pm.Potential (custom likelihood)
    logp = B * LL_observed + (1 - B) * LL_censored
    logpvar = pm.Potential('logpvar', logp.sum())

    trace = pm.sample(init="advi+adapt_diag", return_inferencedata=True)

```