Hi guys. I want to do bayesian model averaging for different modflow models using PyMC. According to this script defining the models, we will predict output basing on inputs and target observed value. However in modflow the predicted value is already calculated using the modflow software. So, I found a problem to construct my models in PyMC and to complete the BMA process. Anyone could help me please ? Thank you in advance

```
with pm.Model() as model_0:
alpha = pm.Normal("alpha", mu=0, sigma=10)
beta = pm.Normal("beta", mu=0, sigma=10)
sigma = pm.HalfNormal("sigma", 10)
mu = alpha + beta * d["neocortex"]
kcal = pm.Normal("kcal", mu=mu, sigma=sigma, observed=d["kcal.per.g"])
trace_0 = pm.sample(2000, return_inferencedata=True)
```

I have read this question multiple times and I am still not sure if I completely get it, but here is an answer to what I understood:

I assume that you have some prior (deterministic) models that outputs certain predictions about a parameter (I think this is d[“neocortex”] correct me if I am wrong). So each model produces a parameter that it calculates. You want to now build up a Bayesian model which uses these predictions. If what I understood is correct one possibility could be as follows. You get the distribution of values for d[“neocortex”] and try to fit some probability distribution (with some continuous distribution available in pymc) that represents this distribution of values. For demonstration purposes lets assume your parameters fits a normal distribution with mean = 5 and sd=1. Then to your model you can add an extra prior of the form

```
p = pm.Normal("p", 5, 1)
```

and replace d[“neocortex”] with that. This would represent your prior scientific knowledge about d[“neocortex”] as a probability distribution. The model would then go on to update this parameter (and the others) further to fit the data best. However while updating this parameter, it would need to “respect” the constraint sets by the parameter prior. So how far it could steer of from your pre-calculated parameters would be an interplay between how strict your prior is and the likelihood of observing the data given your model.

Note that this is not exactly Bayesian model averaging though. The only applications of Bayesian model averaging I have previously seen (and I am not an expert) are when averaged models themselves are Bayesian models and you for instance incorporate an extra “index” or “selection probability” prior that allows you to hop between these models while making predictions. When your models are not Bayesian, but their outputs are deterministic (as I assumed in this case), this was one way that came to my mind.

If this is not what you were looking for, you may want to explain further:

1- What are the “predicted values” coming out of the modflow software (my assumption was these are the d[“neocortex”] values)?

2- Are these models Bayesian (As far as I can understand it involves solving PDEs (without any noise) so I assumed they were producing deterministic outputs) ?

3- What is the relation between these models and the Bayesian model you have presented here?

1 Like