Find fit function with linear and quadratic component

Questions
2

That data looks like it could be pretty well represented with a standard polynomial regression. Have you tried that?
You probably need a cubic fit here due to that dip around x=1.25.

b0 = pm.Normal('b0', 0, sd=10**4)
b1 = pm.Normal('b1', 0, sd=10**4)
b2 = pm.Normal('b2', 0, sd=10**4)
model_error = pm.HalfNormal('error', 10**4)
y_obs ~ pm.Normal('y', mu=b0+b1*np.pow(x,2)+b2*np.pow(x,3), sd=model_error, observed=y)

I just chose relativly flat priors on b0-2, you might want something more informative