You are right that you need to model it with a Bounded Normal. However, in PyMC3 you can not have bounded observed variable (yet).
In general, there are two ways to deal with it:
1, transformed the data so that the observer is now (-∞, ∞), an inverse logistic transform is one of them.
2, model it with a Beta distribution:
sd = pm.Uniform('sd', 0, tt.sqrt(score_est * (1 - score_est)))
# Parameterized the Beta with mu and sd
y_like = pm.Beta('y_like', mu=score_est, sd=eps, observed=audits.FinalScore)
Also, if you know where the percentage comes from (the n and the k), then ideally you should use a Binomial.