PyMC for Monte Carlo Simulation

I am new to Bayesian methods and PyMC, so forgive me if this question might sound odd:

For the sake of simplicity, suppose a linear model. However, the actual model might consist of several parameters.

y = pm.Normal(beta*X, sigma**2) 

beta = [pm.Normal(..., name=intercept), pm.Normal(..., name=slope)] 

X being a controllable physical parameter e.g. concentration, pressure etc. If I have allowed limits on y (e.g. toxicological limits), I would like to know which ranges of X would yield acceptable limits of y. Once I fitted the model parameters (beta), intuitively I would perform a Monte Carlo simulation if the problem cannot be solved analytically. Since PyMC already has efficient samplers built in, I was wondering if the following task or steps could be achieved using PyMC and if such an approach would make sense for you or if you would suggest a different approach?

1.) Fit the model and find the unknown parameters
2.) Convert the input X to a probability distribution to sample from
3.) Use PyMC.sample to draw samples and calculate y
4.) Define an acceptable range for X which yields y values within the acceptable range (e.g. below toxic concentration)

when you do pm.sample(...) the resulting fit represent your unknown parameter as MCMC samples

you just need to more X as well in your model, this is an old blogpost but the idea stands Motif of the Mind | Junpeng Lao, PhD

for calculating y I suspect you are looking for pm.sample_posterior_preditive, which give you prediction conditioned on your observation and inference result

You can also do that use posterior predictive sample, by finding all the y within the acceptable range, and index to X, and compute the range of X[index, :]

1 Like