Issue with custom likelihood

Questions
4

Hi,

Thanks for the suggestion. However, for the scale of the data that prior is reasonably informative (data values range from 1e2 to 1e6).

I found that the divergences actually came about because, while trying out different things to get the likelihood working, I’d switched the prior on the amplitude to a Normal (as opposed to the half Normal). So ruling out non-physical values probably stopped the leap-frog steps wandering in to a region the lead to divergences.

I am now encountering some odd behaviour when using the Arviz HDI routines. Apologise for conflating the two topics in this thread!!

I am using the same model as above (although somewhat modified), but I want a deterministic quantity in there. I will restate so nothing is ambiguous.

def pow_law_mod(x,p1,p2,p3):
    return p1*tt.power(x,-p2)+p3

def likelihood(model, theta1, theta2, theta3, nu):
    '''
    model:    - function
                functional form of curve
    theta1-3: - float
                parameters
    nu:       - int 
                constant
    '''
    def logp_(value):
        mod_ = model(value[0],theta1,theta2, theta3) 
        return -nu*((value[1]/mod_)+tt.log(mod_)+(2/nu-1)*tt.log(value[1]))
    return logp_


with pm.Model() as pow_law_v2:
    
    #prior
    slope = pm.Normal('slope', 1, 1 )
    amp = pm.HalfNormal('amp', 1 )
    const = pm.Normal('const', 1e3, 100)
    
    #likelihood
    like = pm.DensityDist('like', likelihood(pow_law_mod, amp, slope, const,2), observed=[x,test_data])
    mu = pm.Deterministic('mu', pow_law_mod(x, amp, slope, const))

    trace = pm.sample(1000, tune = 2500, target_accept=0.85)
    res = az.from_pymc3(trace)

The model samples fine, and posterior samples give good estimates for actual parameter values.

I am then using the Arviz hdi and plot_hdi functions to generate a HDI for mu. However, when plotting:

az.plot_hdi(x,hdi_data=hdi_data['mu'], ax=ax)

gives some odd curve:

example_curve
example_curve1800×1200 129 KB

Basically it’s the strange dip around 0.002. This curve is the same when using the trace of mu in plot_hdi instead of pre-computing the HDI data.

I checked to see whether it was something to do with the samples from the posterior, by re-creating the posterior curves using the model parameters (i.e. amp, slope and const); which give the grey lines. Also the curves for mu from the trace give the same results.

Finally (and I just checked this before hitting submit), if I plot the upper and lower HDI curves from az.hdi manually:

plt.plot(x,hdi_data['mu'][:][:,0])
plt.plot(x,hdi_data['mu'][:][:,1])

They are also fine.

Am I missing something with the plot_hdi function?