since numpyro supports enumerating discrete latent variables, imputing missing values for discrete covariates should be a possibility (which makes numpyro suitable for many more applied projects!)

Since array shapes will be altered when using parallel enumeration it is not directly evident how to adapt the continuous imputation example to discrete covariates, an example may be helpful

right now we use "twice_nll" as a fit objective and in the test statistic a simple diffence

twice_nll_constrfit - twice_nll_globalfit

but rather we should just to a NLL fit and in the test stat do

2*(nll_constrfit - nll_globalfit)

this will require updating some test reference numbers in the tests

Because / at end of question is missing

Description
Update README and add experiments.

Solution
Create experiments or examples package to work