Unidentified discrete variable

1

I am implementing an intermediate complexity model on the analysis of asthma characteristics taken from Chapter 6 of the book by Winn & Bishop, and I wish to run it in pymc version 4 under numpyro. I am getting a message that NUTS cannot run with discrete variables. I supposedly took care of that by replacing the output of a table by a Bernoulli mixture when the variables were non-observed. However, I still have a problem when the variables are observed. Would this occur if there were missing values that were imputed automatically?

Specifically, in the line below:

            ige_test[ti][al]    = pm.Bernoulli("ige_test"+st, table_ige_[ti][al], observed=data_people[ti]['IgE'][al])

there are imputed value. Would this be sufficient to prevent NUTS from running? Thanks.

Here is my model:

def create_mwe_clean(data, data_people, nb_people, nb_allergens):
    # data (Pandas DataFrame)
    # data_people (observables)
    # Hanlde 4 time frames, single allergy
    skin_test = {}
    ige_test = {}
    sens = {}
    sens[1] = {}
    table_skin_ = {}
    table_ige_ = {}
    table_retain_gain_ = {}
    p_retain = {}
    p_gain = {}
    nb_aller = nb_allergens
    shape = (nb_people, nb_allergens)
    shape = nb_allergens
    
    with pm.Model() as asthma_model:
        p_skin_if_sens = pm.Beta("p_skin_if_sens", 2, 1)
        p_skin_ifnot_sens = pm.Beta("p_skin_ifnot_sens", 1, 2)
        p_ige_if_sens = pm.Beta("p_ige_if_sens", 2, 1)
        p_ige_ifnot_sens = pm.Beta("p_ige_ifnot_sens", 1, 2)
        p_sens_prior = {}
        p_sens_prior[1] = {}
        w = pm.Dirichlet('w', a=np.array([1, 1]))  # 2 mixture weights
        for al in allergens:
            p_sens_prior[1][al] = pm.Beta("p_sens1_"+str(al)+"_prior", 1, 1, shape=nb_people)  
            sens[1][al] = pm.Bernoulli("sens1_"+str(al), p=p_sens_prior[1][al], shape=nb_people)  ### NOT OBSERVED
        
        for ix, ti in enumerate(times):                             
            table_skin_[ti] = {}
            table_ige_[ti] = {}
            skin_test[ti] = {}
            ige_test[ti] = {}
            
            for al in allergens:
                st = "_" + al + "_" + str(ti)

            if ti != 1:
                p_retain[ti] = {}
                p_gain[ti] = {}
                table_retain_gain_[ti] = {}
                sens[ti] = {}
                for al in allergens:
                    st = "_" + al + "_" + str(ti)
                    p_retain[ti][al] = pm.Beta("p_retain"+st, 1, 1) 
                    p_gain[ti][al]   = pm.Beta("p_gain"+st, 1, 1)   
                    (sens_true, sens_false) = table_retain_gain(sens[times[ix-1]][al], p_retain[ti][al], p_gain[ti][al], nb_people, nb_aller)
                    components  = [
                        pm.Bernoulli.dist(p=sens_true, shape=nb_people), 
                        pm.Bernoulli.dist(p=sens_false, shape=nb_people)
                    ]
                    sens[ti][al] = pm.Mixture("sens"+st, w=w, comp_dists=components)

            table_skin_[ti][al] = table_skin(sens[ti][al], p_skin_if_sens, p_skin_ifnot_sens, nb_people, nb_aller) 
            table_ige_[ti][al]  = table_ige(sens[ti][al], p_ige_if_sens, p_ige_ifnot_sens, nb_people, nb_aller) 
                        
            skin_test[ti][al]   = pm.Bernoulli("skin_test"+st, table_skin_[ti][al], observed=data_people[ti]['Skin'][al])
            ige_test[ti][al]    = pm.Bernoulli("ige_test"+st, table_ige_[ti][al], observed=data_people[ti]['IgE'][al])  

    return asthma_model
2

Yes

3

Thanks. I will close this thread and start another asking questions about Mixture models.