Using a “black box” likelihood function with multilevel modeling

Questions
1

Based on the article Using a “black box” likelihood function 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 :slightly_smiling_face:

2

Reinstalling the environment solved the problem :slightly_smiling_face: