# Metropolis with Scipy.fsolve

**URL:** <https://discourse.pymc.io/t/metropolis-with-scipy-fsolve/7300>\
**Category:** Questions\
**Created:** [April 21, 2021, 7:42pm UTC](https://discourse.pymc.io/t/metropolis-with-scipy-fsolve/7300 "2021-04-21T19:42:27Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![RandomVar](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/randomvar/32/4623_2.png) [@RandomVar](https://discourse.pymc.io/u/RandomVar)\
**Post date:** [April 21, 2021, 7:42pm UTC](https://discourse.pymc.io/t/metropolis-with-scipy-fsolve/7300/1 "2021-04-21T19:42:27Z")

</div>

I am trying to use Pymc3 on a model that uses scipy’s fsolve and integrate.quad. I understand that due to not being able to give the sampler gradients, my model cannot use the HMC/NUTS samplers but Metropolis samplers work.

I wrote a small test script to see if I could get Metropolis work with fsolve and I got the following issues when I try sampling.

**Libraries**

```
import numpy as np
import matplotlib.pyplot as plt
import pymc3 as pm
from theano.compile.ops import as_op
import scipy
import seaborn as sn

```

**Functions**

```
def basic_func(y, p):
    return y**2 - p

def test_basic_solver(p):
    y_guess = 1
    sol = scipy.optimize.fsolve(basic_func, y_guess, args=(p), full_output=0)
    return sol[0]

@as_op(itypes=[tt.dscalar], otypes=[tt.dscalar])
def test_basic_solver_OP(p):
   y_guess = 1
    sol = scipy.optimize.fsolve(basic_func, y_guess, args=(p), full_output=0)
    return sol[0]

```

**Synthetic Data**

```
p_real = 4
p_noise = scipy.stats.norm.rvs(loc=0, scale=0.01, size=100)
p_data = p_real + p_noise

p_mean = np.mean(p_data)
p_std = np.std(p_data)

y_data = [test_basic_solver(p) for p in p_data]

```

**Model**

```
 if __name__ == ' __main__':
    with pm.Model() as model_test:
        p_guess = pm.Normal('p', mu=p_mean, sigma=p_std)
        error = pm.HalfCauchy('sigma', 0.1)
    
        obs = pm.Normal('y', mu=test_basic_solver_OP(p_guess), sigma=error, observed = y_data)
        step = pm.Metropolis([p_guess, error])

        trace_trial = pm.sample(draws=200, cores=2, step=step, progressbar=True)

```

**NotImplementedError** : input nd

During handling of the above exception, another exception occurred:

**NotImplementedError** : input nd  
Apply node that caused the error: InplaceDimShuffle{x}(FromFunctionOp{test\_basic\_solver\_OP}.0)

 ![image](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/2X/8/83f6826048e26045adc86c9771381adc2000321c.png)

---

<div class="post-metadata">

**Author:** ![RandomVar](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/randomvar/32/4623_2.png) [@RandomVar](https://discourse.pymc.io/u/RandomVar)\
**Post date:** [April 26, 2021, 11:34pm UTC](https://discourse.pymc.io/t/metropolis-with-scipy-fsolve/7300/2 "2021-04-26T23:34:57Z")

</div>

Nobody has had this issue or tried to sample from a model using Scipy methods?

---

<div class="post-metadata">

**Author:** ![BioGoertz](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/biogoertz/32/4489_2.png) [@BioGoertz](https://discourse.pymc.io/u/BioGoertz)\
**Post date:** [May 13, 2021, 2:37pm UTC](https://discourse.pymc.io/t/metropolis-with-scipy-fsolve/7300/3 "2021-05-13T14:37:15Z")

</div>

I’m not sure if this can do what you need, but if you can reformulate your problem in JAX (e.g. through the [`minimize`](https://jax.readthedocs.io/en/latest/_autosummary/jax.scipy.optimize.minimize.html) or [`custom_root`](https://jax.readthedocs.io/en/latest/_autosummary/jax.lax.custom_root.html), then you can use [this Op](https://discourse.pymc.io/t/theano-op-using-jax-for-lightning-fast-ode-inference/7244). I formulated it for solving ODEs, but it should work with any JAX-friendly function. You still can’t use NUTS for some unkown reason (you _should_ be able to, but there’s a bug somewhere), but you can use ADVI, and probably Metropolis.