# Sum of two custom likelihoods

I am trying to define custom likelihoods with density dist and add them up but it complains that the RVs cannot be added using a + operation (I also tried to use theano.tensor.add and it did not work). I would appreciate it if someone could help me understand how to compute this with pymc3.

def create_model_comorbidity(events, trajectories_time, trajectories_contributions, trajectories_age, comorbidities=[], status_col='status'):
with pm.Model() as model:
BoundedNormal = pm.Bound(pm.Normal, lower=-30, upper=0)
BoundedNormal2 = pm.Bound(pm.Normal, lower=0, upper=0.5)
BoundedNormal3 = pm.Bound(pm.Normal, lower=-2, upper=2)
BoundedHalfNormal = pm.Bound(pm.HalfNormal, lower=0, upper=2)
a0 = BoundedNormal('a0', mu=0, sd=30)
a_sigma = BoundedHalfNormal('a_sigma', sd=10)
a_Sex = BoundedNormal3('a_Sex',
mu=0,
sd=a_sigma
)
b0 = BoundedNormal2('b0', sd=10)
b_sigma = BoundedHalfNormal('b_sigma', sd=10)
b_Sex = BoundedNormal3('b_Sex',
mu=0,
sd=b_sigma
)
mu_a = a0 + a_Sex*events['Sex']
mu_b = b0 + b_Sex*events['Sex']
alpha_com, beta_com = [], []
for curr_com in comorbidities:
alpha_com.append(BoundedNormal3('a_{}'.format(curr_com), mu=0, sd=a_sigma))
beta_com.append(BoundedNormal3('b_{}'.format(curr_com), mu=0, sd=b_sigma))
num_events = trajectories_contributions.shape[2]
print(num_events)
for i in range(num_events):
logalpha = mu_a + trajectories_contributions[:, 1, i] * alpha_com[0]
beta = mu_b + trajectories_contributions[:, 1, i] * beta_com[0]
logalpha += beta*trajectories_age[:, i]
if i == 0 and num_events != 1:
S = pm.DensityDist('S_{}', survival,
observed={'event': np.zeros((events.shape[0],)),
'logalpha': logalpha,
'beta': beta,
'time': trajectories_time[:, i]
}
)
elif i == 0 and i == (num_events - 1):
S = pm.DensityDist('S_{}'.format(i),
survival,
observed={'event': events[status_col],
'logalpha': logalpha,
'beta': beta,
'time': trajectories_time[:, i]
}
)
elif i == (num_events - 1):
print(1)
t = trajectories_time[:, i] - trajectories_time[:, i-1]
S1 = pm.DensityDist('S_{}'.format(i),
survival,
observed={'event': events[status_col].values,
'logalpha': logalpha,
'beta': beta,
'time': t
}
)
else:
t = trajectories_time[:, i] - trajectories_time[:, i-1]
S1 = pm.DensityDist('S_{}'.format(i),
survival,
observed={'event': events[status_col].values,
'logalpha': logalpha,
'beta': beta,
'time': t
}
)

S_tot = S + S1

return model

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()

model = create_model_comorbidity(events, trajectories_time, trajectories_contributions, trajectories_age, comorbidities=['B'], status_col='status')


---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-49-e05f6466f5f3> in <module>
----> 1 model = create_model_comorbidity(events, trajectories_time, trajectories_contributions, trajectories_age, comorbidities=['B'], status_col='status')

<ipython-input-48-2db965ebc090> in create_model_comorbidity(events, trajectories_time, trajectories_contributions, trajectories_age, comorbidities, status_col)
71                                                   )
72                          )
---> 73         S_tot = S + S1
74
75     return model

TypeError: unsupported operand type(s) for +: 'MultiObservedRV' and 'MultiObservedRV'


Anurag

What’s strange to me is that S1 may not even be defined (if num_events == 1) then it looks like only S is defined?

In general, if you have conditionally independent observations like this, where

P[D|\theta] = P[D_1|\theta]P[D_2|\theta]P[D_3|\theta]\dots P[D_k|\theta]

you can use

  for dat_idx in data:
pm.DensityDist('D%d' % dat_idx, my_lik_funcs[dat_idx], observed=my_args[dat_idx])


and the likelihoods are added up by default.

So: assuming there’s no problem with S1 not being defined in your script; just drop “S_tot = S + S1” and it will work just fine.

Thanks @charti. That was good to confirm/know.

Regarding your other question: Yes, this is survival with time dependent hazard functions (step function changes). So, they survive for some time under one hazard function and then they have a change of hazard function and may or may not have an event (i.e., they survived until final time T without undergoing the event or they have an event at t < T). In addition, there is stochasticity associated with the exact time each individual undergoes a change in the hazard function. I probablt could have probably used pm.switch as well but this was another way to code it up and I went with this.