Taking these backwards:
(5)
It’s hard to proffer advice on prior selection or alternative parameterizations without seeing the form of the likelihood function. Given that they’re not particularly identifiable, it means that either the likelihood is largely invariant to them (in which case, use a strong prior for sampling), or varies (largely) with some scalar function g(a,b) - in which case you might consider sampling from g and the preimage {(a,b) : g(a,b)=g} (apologies for using g as a value and a function here).
(4) Is PyMC4 dead or just resting? - #5 by bridgeland
(3) 45 minutes seems like a long time for 4 parameters. Can the likelihood function implementation be optimized?
(2) The values you have plotted are samples from the posterior. You can use the model to generate additional samples from the posterior.
(1) This is a perennial problem that has been asked many times before, and it’s deeply related to your observation that you don’t have a CDF. You have samples from P_{\mathrm{post}}[\theta] = P[\mathcal{D}|\theta]P[\theta]; but if that distribution isn’t closed-form, the only way to sample from P_\mathrm{post}[\alpha] = P[\mathcal{D}'|\alpha, \theta]P_{\mathrm{post}}[\theta]P[\alpha] is via
P_\mathrm{post}[\alpha] = P[\mathcal{D}'|\alpha, \theta]P[\mathcal{D}|\theta]P[\theta]P[\alpha]
i.e., by including both datasets. Solutions posed here, so far, have been to approximate the posterior CDF with some (typically non-parametric) distribution Q_{\mathrm{post}}[\theta] \approx P[\mathcal{D}|\theta]P[\theta].
One possibility is taking a close look at your likelihood function and seeing if it will factor into something of the form
P[\mathcal{D} | \theta] = \lambda(\mathcal{D})f(\mu(\mathcal{D}), \theta)
where \lambda and f are arbitrary functions, but \mu computes sufficient statistics for \theta. Then using P[\theta] \propto f(\mu_0, \theta) is conjugate. This would give you a closed-form posterior an also a mechanism of speeding up your sampling via a custom step.