Hello, I have met same error as yours.
I tried to construct my own likelihood function with gradient refered to [Using a “black box” likelihood function]. But it failed in the “model.compile_dlogp()(ip)”, occuring same NotImplementedError.Then, I tried to run the example codes in colab. It turns out the codes in example doesn’t work 
After a lot of investigations, I figure it out what happen !!!
First, if you delete the sentence in grad():" if m.type.dim!=0 ……", you would meet an error: “LokLIkeWithGrad.grad returned a term with 0 dimensions, but 1 are required”, and
the above sentence: “if hasattr(var, “ndim”) and term.ndim != var.ndim: raise ValueError(……）”
These errors mean that the first term returned from grad was a scalar, but the variable from the node is a vector with shape(1,).
The problem is clarified! And the solution is quite simple:
In make_node() of LogLikeWithGrad(Op), after converting inputs to tensor variables “m=pt.as_tensor(m)”, we could add a sentence “m = m[0]”, that could convert a vector with shape(1,) into a scalar with same dimensioin as the output of grad(). Besides, we could add some dimension judgment statement (such as "if m.type.ndim!=0 ") to avoid potential error.
Hope my experience could help you.