How to condition on a parameter

Questions
1

I am currently following https://am207.github.io/2017/wiki/Lab7_bioassay.html#sampling to explore PyMC3. Roughly, they build a logistic regression model that they use with a binomial likelihood.

with pm.Model() as bioassay_model2:
    alpha = pm.Normal('alpha', 0, sd=100)
    beta = pm.Normal('beta', 0, sd=100)
    theta = invlogit(alpha + beta * log_dose_shared)
    obs_deaths = pm.Binomial('obs_deaths', n=n_shared, p=theta, observed=deaths)

So far so good. Now I found the example in Gelman as well, where Gelman brings up the problem for estimating an LD50. I thought this would be a cool addition to build into my notebook and to explore PyMC3 a bit more. Now, I know that I can calculate the LD50 by - alpha / beta. It is a bit more complicated though, because, as Gelman writes

We summarize the inference on the LD50 scale by reporting two results: (1) the posterior probability that β > 0—that is, that the drug is harmful—and (2) the posterior distribution for the LD50 conditional on β > 0.

Now I can see a path for obtaining (1), because it is kind of taking Pr(alpha, β > 0 | DATA) and integrating over alpha. For (2) I am a bit more lost, because I wonder how I can condition the posterior distribution for the LD50 on β > 0.

So far I added a deterministic distribution for the ld50 into my model:

 ld50 = pm.Deterministic('ld50', -alpha / beta)

This should give me a posterior distribution for ld50, but how can I condition for β > 0?

2

If I understand, you should be able to replace beta in general for something like this with

BoundedNormal = pm.Bound(pm.Normal, lower=0.0)
beta = BoundedNormal('beta',  mu=0, sd = 100)

in fact, this is a special case, and you could (and should) just use beta = pm.HalfNormal('beta', sd=100), which is included as a special case.

3

If I bounded Beta, I would put the assumption of beta being greater than zero into the model. In fact, I would like however, to obtain that probability in a posteriori. I now think that I can actually use the trace of beta to achieve this (ergodic averaging over (trace[‘beta’] < 0).

1 Like