3
$\begingroup$

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

Example of the data I have and the fitted portion

Improve this question
$\endgroup$
7
  • $\begingroup$ I have a hunch that the multivariate delta method could be useful here, applied to the two estimators $\hat{a}$ and $\hat{b}$ you have for true coefficients $a$ and $b$. The tricky part is going to be the fact that $\hat{a}$ and $\hat{b}$ are correlated in a specific way because they follow from linear regression. However, if you have a lot of data points, you could just ignore that correlation. $\endgroup$
    – StijnDeVuyst
    Commented Oct 29, 2019 at 17:41
  • 1
    $\begingroup$ How did you do the fitting? Is your procedure capable of reporting the full covariance matrix of the three estimated parameters? (It's likely those estimates are strongly correlated, so this matters.) $\endgroup$
    – whuber ♦
    Commented Oct 29, 2019 at 17:43
  • $\begingroup$ First thanks for the responses. For the fitting I use a Non-Linear Least-Squares Fitting algorithm in python called lmfit. But the fit in itself is not problematic, what I need is an error in the position of the vertex assuming that the errors in the parameters $a$ and $b$ are ok. Indeed both are correlated and regarding the data points I'm fitting several hundreds (or thousands) in that red window. could you please link me something regarding multivariate delta method? I don't know anything about it... $\endgroup$
    – Emberck
    Commented Oct 29, 2019 at 22:43
  • $\begingroup$ The multivariate delta method requires either you assume the estimates are uncorrelated or that you have the full covariance matrix. Whether or not the fit is problematic and regardless of how much data you have, it's likely the correlations of your estimates are large, so for many applications it would be wise to pursue a solution that estimates and exploits those covariances. $\endgroup$
    – whuber ♦
    Commented Oct 30, 2019 at 14:44
  • 2
    $\begingroup$ My remark above that "if you have a lot of data points, you could just ignore that correlation" is wrong. Sorry about that. $\endgroup$
    – StijnDeVuyst
    Commented Oct 31, 2019 at 18:18

1 Answer 1

Reset to default
0
$\begingroup$

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())
Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.