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

Yes, you can extract the mu from the trace: mu = trace['mu'], do the logistic transform: p = logistic(mu), and add it back to the trace trace.add_values(p), then you can visualize it using traceplot etc.

You Beta model is not working initially because you fix the sd to sd = tt.sqrt(score_est * (1 - score_est)), but for the mu and sd parametrization to work, sd must within the bound of (0, sqrt(mu*(1-mu)))