How to model very complicated likelihood automatically？

Questions v5
2

Trivial stuff PyMC CustomDist can do it for you. Borrowing your example:

import pymc as pm

def dist(mu, sigma, size):
  return (mu + sigma * pm.Normal.dist(size=size)) ** 2

mu = 2
sigma = 1.5
x1 = pm.CustomDist.dist(mu, sigma, dist=dist)

value = 0.7
print(pm.logp(x1, value).eval())  # -1.9362354720843369

If your complicated likelihood function is still relatively simple (mostly chained invertible functions), PyMC can figure it for you automatically. Otherwise you’ll have to step in.

But to take a step back, many times you don’t really need it. You can define an arbitrarily complex function for the mean / standard deviation and then assume a normal likelihood (or some other simple form) around this. Whether this is appropriate depends on your problem obviously.