How to set up a custom likelihood function for two variables

pm.DensityDist needs to evaluate on an observed value, you should pass a function to it instead of a tensor. For example, in this case you should pass likelihood instead of the output of likelihood(theta_1, theta_2):

exp_surv = pm.DensityDist('like', likelihood, observed={'theta_1':theta_1, 'theta_2':theta_2})

Otherwise, you can also add the custom logp directly to the logp as a potential:

like = pm.Potential('like', likelihood(theta_1,theta_2))