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)
```