How can I specify a non-theano based likelihood?

Question

I saw a post from a few days ago by someone else: pymc3 likelihood math with non-theano function. Even though I think the problem at its core is the same, I thought I would ask with a simpler example:

Inside logp_wrap, I put some made up definition of a likelihood function. It depends on the rv and an observation. In this case I could do this with theano operations, but let's say that I want this function to be more complex and so I cannot use theano.

The problem comes when I try to define the likelihood both in terms of an RV and observations. From what I have seen, this format would work if I was specifying everything in 'logp_wrap' as theano operations.

I have searched around for a solution to this, but haven't found anything where this problem is fully addressed.

The problem in my attempt to do this is actually that the logp_ function is correctly decorated, but the logp_wrap function is only correctly decorated for its input, and not for its output, so I get the error TypeError: 'TensorVariable' object is not callable.

Would be great if someone had a solution - don't think I am the only one with this problem.

The theano version of this that works (and uses the same function within a function definition) without the @as_op code is here: https://pymc-devs.github.io/pymc3/notebooks/lda-advi-aevb.html?highlight=densitydist (Specifically the sections: "Log-likelihood of documents for LDA" and "LDA model section")

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
"""
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import pymc3 as pm
from theano import as_op
import theano.tensor as T
from scipy.stats import norm

#Some data that we observed
g_observed = [0.0, 1.0, 2.0, 3.0]

#Define a function to calculate the logp without using theano
#This as_op is where the problem is - the input is an rv but the output is a
#function. 
@as_op(itypes=[T.dscalar],otypes=[T.dscalar])
def logp_wrap(rv):
    #We are not using theano so we wrap the function.
    @as_op(itypes=[T.dvector],otypes=[T.dscalar])
    def logp_(ob):
        #Some made up likelihood -
        #The key here is that lp depends on the rv input and the observations        
        lp = np.log(norm.pdf(rv + ob))
        return lp
    return logp_

hb1_model = pm.Model()
with hb1_model:
    I_mean = pm.Normal('I_mean', mu=0.1, sd=0.05)
    xs = pm.DensityDist('x', logp_wrap(I_mean),observed = g_observed)

with hb1_model:
    step = pm.Metropolis()
    trace = pm.sample(1000, step)