I see, why not model the lognormal directly with the mean of the binomial: n * p or np.log(n) + np.log(p) on the logscale?
llike = pm.Lognormal("llike", mu=mu + pm.math.log(n), sigma=sd, observed=data)
Or is the data generating process really the one with a noisy binomial? In that case you can model it like that of course (with the caveats you can’t use NUTS alone)