using hierarchical models with random slope resutling in sampling problem

Hi, I am currently working with the HSSM package (v 0.2.3) and have encountered an issue while handling my own data, including the "cavanagh_theta" data provided by the package. Specifically, when using hierarchical models with random slope (e.g., stim|participant_id), I am facing the following problem: the sampling process terminates quickly, but the obtained samples are mostly invalid, resulting in all NaN values for r_hat. However, if remove the random slope, the sampling process will be more normal and the results of the sampling will also be acceptable. Could you please investigate this issue? Can you check this issue? Currently, there are still some problems with HSSM when dealing with hierarchical models, so it is still difficult to use it in formal research.

hierarchical_ddm_model_reg_vazt = hssm.HSSM(
        model = "ddm",
        data = cav_data,
        prior_settings="safe",
        loglik_kind = "approx_differentiable", #"analytical"      
        include = [
            {
                "name": "v",
                "formula": "v ~ 0 + stim + (0 + stim|participant_id)",      
                "link": "identity",
            },
            {
                "name": "a",
                "formula": "a ~ 0 + stim + (0 + stim|participant_id)",      
                "link": "identity",
            },
            {
                "name": "z",
                "formula": "z ~ 0 + stim + (0 + stim|participant_id)",      
                "link": "identity",
            },      
            {
                "name": "t",
                "formula": "t ~ 0 + stim + (0 + stim|participant_id)",      
                "link": "identity",
            },                                                 
                     
        ],
    )
    
    print(hierarchical_ddm_model_reg_vazt)

or the following code

# Define a basic hierarchical model with trial-level covariates
hierarchical_ddm_model_reg_vazt = hssm.HSSM(
    model = "ddm",
    data = cav_data,
    loglik_kind = "approx_differentiable", #"analytical"      
    include = [
        {
            "name": "v",
            "formula": "v ~ 1 + stim + (1 + stim|participant_id)",
            "link": "identity",
        },               
        {
            "name": "a",
            "formula": "a ~ 1 + stim + (1 + stim|participant_id)",
            #"formula": "v ~ (1|subj_idx) + theta",
            "link": "identity",
        },  
        {
            "name": "z",
            "formula": "z ~ 1 + stim + (1 + stim|participant_id)",
            "link": "identity",
        }, 
        {
            "name": "t",
            "formula": "t ~ 1 + (1|participant_id)",
            #"formula": "v ~ (1|subj_idx) + theta",
            "link": "identity",
        },                                         
    ],
)

print(hierarchical_ddm_model_reg_vazt)

hierarchical_ddm_model_reg_vazt.sample(
    sampler="nuts_numpyro", 
    chains=2,  # how many chains to run
    cores=4,  # how many cores to use
    draws=1000,  # number of draws from the markov chain
    tune=500,  # number of burn-in samples
    #target_accept = 0.90,
    )
    

Start with the HSSM.HSSM construction using the hierarchical formulas in the issue, especially random slopes such as (0 + stim|participant_id), and reproduce the behavior with cavanagh_theta or cav_data. Run sample() with sampler="nuts_numpyro" and inspect why random-slope models produce invalid samples and NaN r_hat values while models without random slopes do not. Done means hierarchical random-slope sampling returns valid diagnostics.