Pymc3 approach

foo.py

%matplotlib inline

import pymc3 as pm

with pm.Model() as model:
    _a = pm.Uniform('a', lower=-10, upper=0)
    _h = pm.Uniform('h', lower=4, upper=12)
    _k = pm.Uniform('k', lower=18000, upper=30000)
    _sigma = pm.Uniform('sigma', lower=5, upper=25)

    predicted_y = _a * (x_data - _h)**2 + _k
    prediction_error = observed_y - predicted_y
    
    y_probability = pm.Normal('y_probability', mu=0, observed=prediction_error, sd=_sigma)
    
    start = pm.find_MAP() # Find starting value by optimization
    step = pm.NUTS(scaling=start) # Instantiate MCMC sampling algorithm
    trace = pm.sample(25000, step, start=start)
    
pm.traceplot(trace)