How to model observed percentages (bounded from 0 to 1)

Questions
4

Thank you for your replies! I tried transforming my FinalScore with a logit function and then fit the Normal to that. It worked but then it makes interpreting the traceplot graphs difficult. Would it be possible to transform the trace to be able to interpret it more easily?

Here is my code when I try to fit it to a Beta distribution and I will also include the exception stack trace as an attachment because it is huge…

with pm.Model() as hierarchical_model:
    # Hyperpriors
    mu_a = pm.Normal('mu_alpha', mu=0, sd=1)
    sigma_a = pm.HalfCauchy('sigma_alpha', beta=1)
    mu_b = pm.Normal('mu_beta', mu=0, sd=1)
    sigma_b = pm.HalfCauchy('sigma_beta', beta=1)
    
    # Intercept for each store, distributed around group mean mu_a
    a = pm.Normal('intercept', mu=mu_a, sd=sigma_a, shape=len(uniqueStores))
    # Intercept for each store, distributed around group mean mu_a
    b = pm.Normal('b', mu=mu_b, sd=sigma_b, shape=len(uniqueStores))
    
    # Model error
    eps = pm.HalfCauchy('eps', beta=1)
    
    # Expected value
    score_est = a[storeIDX] + (b[storeIDX] * audits['Shift1Score'].values)
    
    # Data likelihood
    sd = tt.sqrt(score_est * (1 - score_est))
    
    y_like = pm.Beta('y_like', mu=score_est, sd=sd, observed=audits.FinalScore)

    trace = pm.sample(progressbar=True, njobs=4)

exception.txt (62.8 KB)