1
$\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

$\endgroup$
6
  • $\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$ Oct 29 '19 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
    Oct 29 '19 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
    Oct 29 '19 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
    Oct 30 '19 at 14:44
  • 1
    $\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$ Oct 31 '19 at 18:18
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())
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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