Custom distribution by pm.DensityDist not coincide with built-in distribution in PyMC4

Questions v5
1

I’m new to pymc. Thank you for your great efforts for this amazing library.
Now I’m learning the behavior of custom distribution defined by pm.DensityDist in v4.1.3. For this, I made two equivalent simple gaussian models as follows.

Both models coverge well, but the latter custom model shows larger value of estimation for sigma, probably due to my misunderstanding of the usage of pm.DensityDist. Could you specify any reasons?

Any suggestions are welcome. Thank you in advance.

import arviz as az
import numpy as np
import pymc as pm
import aesara.tensor as at

RANDOM_SEED = 8927
rng = np.random.default_rng(RANDOM_SEED)

#create data from true distribution
data = rng.normal(loc=5,scale=2, size=1000)
print( data.mean(), data.std() )

#1st model using pm.Normal()
model_builtin = pm.Model()
with model_builtin:
    mu = pm.Normal("mu", mu=0, sigma=5)
    sigma = pm.HalfNormal("sigma", sigma=5)
    obs = pm.Normal("obs",mu=mu,sigma=sigma, observed=data)
    idata_builtin = pm.sample()

#2nd model using pm.DensityDist
def _logp(obs, mu, sigma):
    z = (obs-mu)/sigma
    return -at.sum(z**2) / 2

model_custom = pm.Model()
with model_custom:
    mu = pm.Normal("mu", mu=0, sigma=5)
    sigma = pm.HalfNormal("sigma", sigma=5)
#    obs = pm.Normal("obs",mu=mu,sigma=sigma, observed=data)
    obs = pm.DensityDist("obs", mu,sigma, logp=_logp, observed=data) #only changed here
    idata_custom = pm.sample()

#check posteriors
display(pm.summary(idata_builtin))
display(pm.summary(idata_custom))
2

I think your normal logp is wrong. This is the one used by PyMC aeppl/logprob.py at 8f65562d00025b955e3d7d6f52609c54b736afc4 · aesara-devs/aeppl · GitHub

Ignoring the optional normalization constant, I see the -log(sigma) term is missing from your definition.

3

Thank you so much for your quick reply.
Before the 1st post, I tried

def _logp(obs, mu, sigma):
    z = (obs-mu)/sigma
    return -at.sum(z**2) / 2 - at.log(2*np.pi*sigma**2) /2

and had a similarly biased posterior. But I should have included the normalization constant inside the at.sum() … silly mistake.

Now, both

    return at.sum(-z**2 / 2 - at.log(sigma) )

and

    return -z**2 / 2 - at.log(sigma)

work fine. Thank you!