# Using a “black box” likelihood function with multilevel modeling

**URL:** <https://discourse.pymc.io/t/using-a-black-box-likelihood-function-with-multilevel-modeling/4762>\
**Category:** Questions\
**Tags:** theano\
**Created:** [March 29, 2020, 11:08am UTC](https://discourse.pymc.io/t/using-a-black-box-likelihood-function-with-multilevel-modeling/4762 "2020-03-29T11:08:03Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![Philipp](https://avatars.discourse-cdn.com/v4/letter/p/46a35a/32.png) [@Philipp](https://discourse.pymc.io/u/Philipp)\
**Post date:** [March 29, 2020, 11:08am UTC](https://discourse.pymc.io/t/using-a-black-box-likelihood-function-with-multilevel-modeling/4762/1 "2020-03-29T11:08:03Z")

</div>

Based on the article [Using a “black box” likelihood function](https://docs.pymc.io/notebooks/blackbox_external_likelihood.html) I would like to use a blackbox model including a likelihood function. In addition I want to use a hierarical approach for my unknown distribution.

Enclosed a simple example to illustrate the structure.

**Generate the data**

```
def my_model(theta):
    m = theta
    return m

M = 100 # number of samples
sigmatrue = 0.1 # standard deviation of noise

mmeantrue = 10.0 # true mean
msigmatrue = 1.0 # true sigma

m = mmeantrue + np.random.randn(M) * msigmatrue
truemodel = my_model(m)

data = sigmatrue * np.random.randn(M) + truemodel

```

**Likelihood function and Theano Op**

```
def my_loglike(theta, data):
    m = theta
    sigma = 0.1

    model = my_model(m)
    logl = -(0.5 / sigma **2) * np.sum((data - model)** 2)

    return logl

class LogLikeGrad(tt.Op):
    itypes = [tt.dvector]
    otypes = [tt.dvector]

    def __init__ (self, loglike, data):
        self.likelihood = loglike
        self.data = data

    def perform(self, node, inputs, outputs):
        theta, = inputs

        def lnlike(values):
            return self.likelihood(values, self.data)

        eps = np.sqrt(np.finfo(float).eps)
        grads = scipy.optimize.approx_fprime(theta, lnlike, eps * np.ones(len(theta)))

        outputs[0][0] = grads

class LogLike(tt.Op):
    itypes = [tt.dvector]
    otypes = [tt.dscalar]

    def __init__ (self, loglike, data):
        self.likelihood = loglike
        self.data = data
        self.logpgrad = LogLikeGrad(self.likelihood, self.data)

    def perform(self, node, inputs, outputs):
        theta, = inputs
        logl = self.likelihood(theta, self.data)

        outputs[0][0] = np.array(logl)

    def grad(self, inputs, g):
        theta, = inputs
        return [g[0] * self.logpgrad(theta)]

```

**Pymc3 model**

```
# create Op
logl = LogLike(my_loglike, data)

# use PyMC3 to sampler from log-likelihood
with pm.Model():
    # uniform priors on m and c
    mmean = pm.Uniform('mmean', lower=5.0, upper=20.0)
    msigma = pm.Uniform('msigma', lower=0.5, upper=2.0)
    theta = pm.Normal('m', mu=mmean, sd=msigma, shape=M)

    pm.DensityDist('likelihood', lambda v: logl(v), observed={'v': theta})

    trace = pm.sample(3000, chains=1, tune=1000, step=pm.NUTS(), discard_tuned_samples=True)

```

When I want to run the script I get the following error:  
`Exception: ('Compilation failed (return status=1): C:\\Users\\*\\AppData\\Local\\Temp\\4\\ccptC33H.s: Assembler messages:\r. C:\\Users\\*\\AppData\\Local\\Temp\\4\\ccptC33H.s:282: Error: invalid register for .seh_savexmm\r. ', '[Elemwise{EQ}(<TensorType(float64, vector)>, <TensorType(int8, (True,))>)]')`

Any idea why the compilation fails? I’m happy for every tip 🙂

---

<div class="post-metadata">

**Author:** ![Philipp](https://avatars.discourse-cdn.com/v4/letter/p/46a35a/32.png) [@Philipp](https://discourse.pymc.io/u/Philipp)\
**Post date:** [April 17, 2020, 7:21am UTC](https://discourse.pymc.io/t/using-a-black-box-likelihood-function-with-multilevel-modeling/4762/2 "2020-04-17T07:21:05Z")

</div>

Reinstalling the environment solved the problem 🙂