Hi there,
I’m looking for help defining an appropriate likelihood function, which accepts multiple inputs. I’m fairly new to pymc3
and Bayesian modelling in general, so please excuse any oversights on my end. Also, excuse some of the nomenclature; most of my exposure to Bayesian methods in practice has been through geophysical inversions literature.
I have a mass balance model, in my case a simple representation of the surface input and outputs for a glacier system, for which I need to invert for a set of unknown model parameters. The final result of the mass balance model Y is the sum of three components ( i.e. Y = X_1 + X_2 - X_3).
I initially attempted my inversion solely inverting for Y from observations Y_{\rm obs}. While I was able to reproduce Y_{\rm obs}, my posterior distributions were nonphysical for the climatic region I’m working in. I know this could be rectified by defining tighter priors, but many of these parameters have well define distributions from previous literature which are the most scientifically defensible.
As an alternative approach, in hopes of more physical results, I am now attempting the inversion for the individual components of the model X_1, X_2, X_3 based on observations X_{1_{\rm obs}}, X_{2_{\rm obs}}, X_{3_{\rm obs}}. Where I need help defining an appropriate likelihood function which accepts multiple inputs.
As an initial approach I’ve started with:
# Define Forward model (wrapped through theano)
with model:
R, A, M = PDD_forward.forward(z_obs, f_s_prior, C_prior, f_r_prior, grad_a, A_m_prior)
# net balance [m i.e. / yr ]
B = A + R - M
# Define likelihood
with model:
# Individual likelihood functions for each component
R_est = pm.Normal("R_est", mu=R, sigma=R_sigma, observed=R_obs)
A_est = pm.Normal("A_est", mu=A, sigma=A_sigma, observed=A_obs)
M_est = pm.Normal("M_est", mu=M, sigma=M_sigma, observed=M_obs)
potential = pm.Potential("obs", R_est.sum()*A_est.sum()*M_est.sum())
but I’m not sure how appropriate this is. I settled on this as my initial approach by modifying (and simplifying) the example shown here.
A little unclear whether this is the proper way to go about this. While looking through the examples on the pymc3
website, I noticed the Non-linear Differential Equations tutorial from the pymc3.ode
example notebook is at least analogous to my problem, in so far as the forward model has multiple outputs. In that case there is a single likelihood function which takes in multidimensional inputs. Since the code is using the pymc3.ode
API I haven’t been able to 1:1 translate over to my example. I’ve looked at the source code for pymc3.ode.DifferentialEquation
to understand a bit more of what’s returned and see how I could try mimicking it. This seems somewhat complicated, so I’d thought it be best to reach out for help before diving to deep down the rabbit hole.
Thank you for taking to time to look at this. If there’s anything left out and can further clarify please let me know.
Best,
Andrew