Newton Method with Random Variables

I’ll try to briefly explain my problem: I have a Bayesian inference pipeline over hydrodynamic models.
I test my model by coupling ‘sub’ models present in the bibliography to evaluate each part of them, at a certain moment I have to find a variable as a function of the observed variables, e.g I have x and y observed, but I want to find z.
Generally I can manipulate the expression and write my variable of interest as a function of the others, but not this time. For a simple multi-objective optimization I can use Newton’s method or any other numerical method to find the root of functions and it works fine. But using PyMC3 doesn’t work as the methods cant find a solution by passing parameters as random variables (theano structures) to the optimizer.

Just to illustrate, having lhs as my known information and params my random variables from PyMC, with an optimization (where there are no RVs) I was able to simply do the following and it worked:

def f(s, p):
   return p - self.krw_parameterized(s, params)
se = np.array([opt.fsolve(f, 0.5, args=(l))[0] for l in lhs])

Do you guys have any suggestions on how to deal with this issue?

But using PyMC3 doesn’t work as the methods cant find a solution by passing parameters as random variables (theano structures) to the optimizer.

Since the optimizer expects functions defined over Python numeric types or Numpy arrays, this isn’t surprising. It sounds like you want to run a PyMC3 model, get the results from the inference, and then pass it to the optimizer. If that’s the case, can you just use sampled values from pm.sample()? Do you need to perform optimization for every MCMC iteration?

If you do indeed need to do optimization within the main inference loop, I would also like to suggest writing out your objective function as a pm.Potential. If it’s differentiable and not too nasty to deal with, you may be able to directly optimize it by using the default No U Turn sampler. It’s important to remember that, in some sense, NUTS is performing a variant of gradient ascent itself, so it is perfectly capable of rendering solutions which optimize a target objective function which depends on other random variables.

Yea, this isn’t surprising. It’s also weird thinking mathematically about finding the root of a function of random variables, so in that case I’d have to sample points just like in the optimizer through a search algorithm.
And that’s exactly why I’ve already explained the reasons for the problem in the topic, looking for some suggestion to deal with it.

And yes again, I need an optimization for every MCMC iteration. It’s an inverse problem of uncertainty quantification, where I can usually isolate the necessary variables algebraically, but not this time.

What did I do in that time interval? Well I found this decorator method theano.compile.ops.as_op and passed to my function that propagates the randomness

@as_op(itypes=[tt.dscalar,tt.dscalar,tt.dscalar,tt.dscalar,tt.dscalar,tt.dscalar,tt.dscalar], otypes=[tt.dvector]) 

The result is strange but seems to have worked around the problem somehow. I don’t know if it is the best solution.