Find Error in the vertex of a fitted parabola I have a set of data (a Cross-Correlation function) to which I've adjusted a parabola (using lmfit python package ), from the fit I got the values of the parameters and their error:


*

*(Model): $f(x; a, b, c) = a x^2 + b x + c $
$ a \pm \delta a$
$ b \pm \delta b$
$ c \pm \delta c$
the only thing I'm interested with the fit is the position of the vertex of the parabola, but I need to have a measure of the error of such position... and I'm unsure of how to obtain it.
Initialy I thought that since the vertex position is defined by $v=\frac{-b}{2a}$ I just need to use the standard error propagation formula $\delta v = |\frac{\partial v}{\partial a}|\delta a +|\frac{\partial v}{\partial b}|\delta b $ but I'm unsure the at is the correct approach... 
Any help will be most appreciated! Thanks

 A: I'd suggest you to use Stan (https://mc-stan.org/), you can build your statistical model and infer the parameters by sampling: i.e. stan can sample from the posterior distribution of the parameters. Having samples, you can calculate estimates of the mean and standard deviation of your model parameters.
It can be a little trickier to get Stan to work. Below I present my code using python and pystan.
For the Stan model, it is crucial to have prior distribution. Therefore, by any preferred method (e.g. maximum likelihood) I get prior guesses of the parameters a_prior, b_prior, p_prior. Then I assume that prior distribution is gaussian around these means with e.g. 30% std (error_precentage=0.3). The model that I'm fitting is:
a ~ normal(a_prior, square(a_prior*error_precentage));
b ~ normal(b_prior, square(b_prior*error_precentage));
p ~ normal(p_prior, square(p_prior*error_precentage));
    
y ~ normal(a-b*square(x-p), square(y_std));

Full code:
import pystan

stan_code = """
data {
    int N;
    vector[N] y;
    vector[N] x;
    real y_std;
    real a_prior;
    real b_prior;
    real p_prior;
    real error_precentage;
}
transformed data {
    
}
parameters {
    real a;
    real b;
    real p;
}
transformed parameters {
}
model {
    a ~ normal(a_prior, square(a_prior*error_precentage));
    b ~ normal(b_prior, square(b_prior*error_precentage));
    p ~ normal(p_prior, square(p_prior*error_precentage));
    
    y ~ normal(a-b*square(x-p), square(y_std));
}
generated quantities {
}
"""

fit = sm.sampling(data={"N":len(temp), "x":x, "y":y, "y_std":y.std(), 
                        "a_prior":a_prior, "b_prior":b_prior, 
                        "p_prior":p_prior, "error_precentage":error_precentage},
                        iter=5000, chains=8, 
                        )
lb = fit.extract(pars=["a", "b", "p"],permuted=False, inc_warmup=False);

# the statistics for the vertex of parabola
np.mean(lb["p"].flatten()), np.std(lb["p"].flatten())

