Having problems defining random function for custom distribution with mixture of observed and random variables

Questions
6

Since your logp of the observed is:

def survival(event, logalpha, beta, time):
    log_h = logalpha + beta * time
    H = (np.exp(logalpha)/beta) * (np.exp(beta * time) - 1)
    return ((event * log_h) - H).sum()

what I would do (not sure if it is correct or not!) is get a temporal function of the latent alpha and beta:

def survival2_random(logalpha, beta, time):
    log_h = logalpha + beta * time
    H = (np.exp(logalpha)/beta) * (np.exp(beta * time) - 1)
    logp_temporal = np.stack([(0. * log_h) - H, (1. * log_h) - H]) <== observed being 0 or 1
    p_temporal = np.exp(logp_temporal)
    p_norm = p_temporal/p_temporal.sum(axis=1, keepdim=True)
    return np.random.bernoulli(p_norm)

# and call it like:
for point in trace:
    randomsample = survival2_random(point['logalpha'], point['beta'], time)

I did not run the code so there might be mistake :wink: