'Returns a scipy distribution object corresponding to the
distribution of uncensored endog (duration) values for each
case.
Parameters
params : array-like
The proportional hazards model parameters.
Returns
A list of objects of type scipy.stats.distributions.rv_discrete
Notes
The distributions are obtained from a simple discrete estimate
of the survivor function that puts all mass on the observed
failure times within a stratum.'

def get_distribution(self, params):

    from scipy.stats.distributions import rv_discrete
    surv = self.surv
    bhaz = self.baseline_cumulative_hazard(params)
    (pk, xk) = ([], [])
    for stx in range(self.surv.nstrat):
        exog_s = surv.exog_s[stx]
        linpred = np.dot(exog_s, params)
        if (surv.offset_s is not None):
            linpred += surv.offset_s[stx]
        e_linpred = np.exp(linpred)
        pts = bhaz[stx][0]
        ichaz = np.outer(e_linpred, bhaz[stx][1])
        usurv = np.exp((- ichaz))
        z = np.zeros((usurv.shape[0], 1))
        usurv = np.concatenate((usurv, z), axis=1)
        probs = (- np.diff(usurv, 1))
        pk.append(probs)
        xk.append(np.outer(np.ones(probs.shape[0]), pts))
    mxc = max([x.shape[1] for x in xk])
    for k in range(self.surv.nstrat):
        if (xk[k].shape[1] < mxc):
            xk1 = np.zeros((xk[k].shape[0], mxc))
            pk1 = np.zeros((pk[k].shape[0], mxc))
            xk1[:, 0:xk[k].shape[1]] = xk[k]
            pk1[:, 0:pk[k].shape[1]] = pk[k]
            (xk[k], pk[k]) = (xk1, pk1)
    xka = (np.nan * np.ones((len(self.endog), mxc)))
    pka = (np.ones((len(self.endog), mxc), dtype=np.float64) / mxc)
    for stx in range(self.surv.nstrat):
        ix = self.surv.stratum_rows[stx]
        xka[ix, :] = xk[stx]
        pka[ix, :] = pk[stx]
    dist = rv_discrete_float(xka, pka)
    return dist
