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

In Binomial, the n could be a vector of trial Ns, so in you case, after
score_est = pm.math.sigmoid(a[storeIDX] + (b[storeIDX] * audits.Shift1Score.values))
you can do something like:

    #likelihood
    y_like = pm.Binomial('observed', n=audits.trial_n, p=score_est, observed=audits.observed_n)