Difficulties in using .dist and pm.Potential for inference

@jessegrabowski, many thanks for checking! I really like your solution.

However, I tried to run a similar code, but I’ve got many divergences 20221016-pymc-question.ipynb · GitHub

Just to confirm, is it not about priors or something like that?

with pm.Model() as model:    
    # Data
    T = tt.as_tensor_variable(counts.shape[0], dtype="int64") 
    D = tt.as_tensor_variable(counts.shape[1] - 1, dtype="int64") 
    D_grid = at.arange(0, D + 1)
    
    counts_matrix = tt.as_tensor_variable(counts, dtype="int64")
    
    # Priors
    μ = pm.HalfNormal('μ', 5.0, initval=2.0)
    σ = pm.HalfNormal('σ', 5.0, initval=2.0)
    
    # Likelihood
    ### Parameters of the lognormal distribution
    σlog = tt.sqrt(tt.log((σ / μ)**2 + 1))
    μlog = tt.log(μ) - σlog**2 / 2
    
    log_norm = pm.Lognormal('log_norm', mu = μlog, sigma = σlog)
    def cdf(x):
        return tt.exp(pm.logcdf(log_norm, x))
    π, _ = ae.scan(lambda t: cdf(t + 0.5) - cdf(t - 0.5), sequences = [D_grid])
    
    pm.Deterministic('p', π)

    log_like = pm.Multinomial('obs', 
                              n=counts_matrix.sum(axis=1), 
                              p = π / π.sum(),
                              observed=counts_matrix)
    
    trace = pm.sample(tune=3000, chains=4, init='advi')