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

Questions
5

Alright I managed to make the Beta posterior work. I am a software engineer and still learning bayesian stats for fun so… I saw this post and transformed my code like his and used a pretty small standard deviation. Made sure my observed variable did not contain values exactly 0 or 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.HalfNormal('eps', sd=0.1)

# Expected value
u = pm.math.invlogit(a[storeIDX] + (b[storeIDX] * audits['Shift1Score'].values))

# Data likelihood    
y_like = pm.Beta('y_like', mu=u, sd=eps, observed=audits.AdjustedScore)

Thank you!