Hi, I am new to PyMC3 and need some help modeling a custom likelihood function. What I am trying to model is essentially the likelihood of getting zj failures out of a total of nj buildings:

But here we replace the p in the binomial with the function of a vulnerability curve. The log likelihood, then, is:

The goal is to try to estimate the parameters theta and beta using observations of xj, nj, and zj.

It seems to me that I probably need to define a custom likelihood function using DensityDist (i.e., hard code the entire expression for the likelihood function) for this but I am wondering if there is a more elegant way to do this using the pm.Binomial likelihood? For example, I would like to do something like:

but am not sure how I would model the input (custom) p. This is the standard normal CDF.