Thank you for your detailed reply, @jlindbloom. It helped with my understanding.
It is not straightforward to me on how to compute gradients using an existing trained ML model (from sklearn/tensorflow/pytorch) that doesn’t have Theano as the backend. Converting the trained model to Theano may be the way to go. I will leave this to expert.
Using a Bayesian neural network to obtain probability distributions for a set of parameters \theta^* was my first intuition in solving the problem. After much deliberation, it is easier to explain the posterior probability of \theta^* obtained using the mapping from \theta^* \rightarrow \theta (this mapping is described by a semi-analytical physics model) than to explain using the mapping from \theta \rightarrow \theta^* from learned weights and bias distributions of BNN.