Comparing Multinomial Models/Multinomial Distribution/LOO

1

I have two models I am applying to the same dataset. I would like to compare the fit of the two models to the data and justify why one model is a better fit than the other. Posterior predictive checking is not favoring one model over the other.

Both models predict probabilities associated with many classes and then the multinomial distribution is used to link the probabilities with the frequencies observed in the data.

The models run well but when I try to apply leave-one-out to a trace, I get error-messges.

     Estimate       SE

elpd_loo -204.77 0.00
p_loo 32.87 -

There has been a warning during the calculation. Please check the results.

Pareto k diagnostic values:
Count Pct.
(-Inf, 0.5] (good) 0 0.0%
(0.5, 0.7] (ok) 0 0.0%
(0.7, 1] (bad) 0 0.0%
(1, Inf) (very bad) 1 100.0%

UserWarning: Estimated shape parameter of Pareto distribution is greater than 0.7 for one or more samples. You should consider using a more robust model, this is because importance sampling is less likely to work well if the marginal posterior and LOO posterior are very different. This is more likely to happen with a non-robust model and highly influential observations.
warnings.warn(
\lib\site-packages\arviz\stats\stats.py:842: UserWarning: The point-wise LOO is the same with the sum LOO, please double check the Observed RV in your model to make sure it returns element-wise logp.
warnings.warn(

I have managed to replicate the problem in a simple example(bottom of this message)

I can’t say I understand the error-message, I could be missing something obvious.

Am I using LOO correctly? Is there a better way to compare the fit of the two models? Any help or advice is appreciated.

Thanks,
Wayne Hajas

-- coding: utf-8 --

“”"
Demonstrate LOO-analysis with multinomial distribution

“”"

import aesara
import aesara.tensor as at
import arviz as az
import matplotlib.pyplot as plt
import numpy as np
import pymc as pm
import pandas as pd

RANDOM_SEED = 20230424
from numpy.random import seed
seed(RANDOM_SEED)
aesara.config.exception_verbosity=‘high’
aesara.config.optimizer=‘None’
aesara.config.traceback__limit=3

def FMN(SampleFreq, RANDOM_SEED, niter=10000,n_init=10000):
nclass=len(SampleFreq)
SampleSize=sum(SampleFreq)

#Define model
MNmodel=pm.Model()
with MNmodel:

    #Variability of probabilities (logit scale)
    lnsigma=pm.Normal('lnsigma',0,sigma=1)
    sigma=at.exp(lnsigma)
    
    #Logistic of probabilities(unnormalized)
    LogProb=pm.Normal('LogProb', sigma=sigma,shape=nclass)
    
    #Un-normalized Probabilities
    UnNormProb=at.exp(LogProb)
    
    #Normalize Probabilities to sum to one
    Prob=UnNormProb/at.sum(UnNormProb)
    
    #Apply Multinomial
    Obs=pm.Multinomial('Obs',SampleSize,Prob, observed=SampleFreq)
return(MNmodel)

if name == “main”:

SampleFreq=[0,1,0,0,1,2,2,1,0,0,0,0,1,1,1,3,1,5,12,23,4,1,0,2,2,2,1,3,1,2,3,0,\
        1,1,0,3,5,3,5,4,2,6,4,11,3,2,4,3,4,3,2,5,2,4,1,3,1,1,3,4,4,0,1,4,6,\
                1,6,11,13,5,6,6,6,14,21,7,4,3,3,14,6,1,2,22,39,1,2,5,2,2,1,\
                        1,2,2,1,2,1,2,2,2,0,1,0,2,8,1]
niter,n_init=10000 ,10000
    
dummy=FMN(SampleFreq, RANDOM_SEED, niter=niter,n_init=n_init)
with dummy:
    step = pm.NUTS()
    trace=pm.sample(niter,n_init=n_init,cores=1,random_seed=(RANDOM_SEED),target_accept=0.98,idata_kwargs={"log_likelihood": True})
    
    az.plot_energy(trace)
    print(az.summary(trace))
    pooled_loo = az.loo(trace)
    print(pooled_loo)