Thanks for getting back to me on this @cluhmann. I tried a skewed-normal distribution, in the following way:
xx = physDist.copy()
with pm.Model() as model_skewNormal:
γ = pm.Normal('γ', mu=1.5, sigma = 0.5) ## analagous to initial spread of variance (natural log of)
δ = pm.Normal('δ', mu=-0.0006, sigma=0.005) ## slope, rate of tightening of variance
κ = pm.Normal('κ', mu=200, sigma = 10) ## same old k
α = pm.Normal('α', mu=-1.0, sigma = 0.5)
μ = pm.Deterministic('μ', xx/(xx+κ))
ε = pm.Deterministic('ε', np.power(np.e, -γ + δ*xx))
y_pred = pm.SkewNormal('y_pred', mu=μ, sd=ε, alpha=α, observed=bcDist)
trace_sn = pm.sample(init='adapt_diag')
This works, but the results are really far shifted from my priors, and the mean function \mu of the posterior sits pretty far out from the rest of the distribution. I guess this is possible because of the skew? Using the skewed normal distribution pushed the asymptote into an almost right-angle pattern. I think because the low X values \{x \lt 500\} in the steep area of the asymptote are difficult to fit with this skewed distribution. Looks like this :
Orange areas representing HPDs of 50 and 94 credible intervals. The general fit of the skewed \epsilon model is pretty poor, R2= 0.13 .
This can be compared to when I assign a non-skewed normal distribution to \epsilon, which doesn’t respect the theoretical limit of my data ( y\leq 1 ) but that predicts a lot more of the variance around the mean ( R2=0.53):
To keep it simple I think that a non-skewed noising will do just fine for my purposes, maybe best to leave good-enough alone here. So I think you answered the “how-to-do-it” of my question, for anyone else who might end up here, and I figured out that my idea wasn’t great for my data 

