I have a bunch of data points $(x, y)$, and I know that they fit well to a model of the form $y = a + bx + c x^2$, with $a \approx 0.01, \ b \approx 1\ \textrm{and}\ c \lesssim 0.1$. I'd like to fit the actual values of the parameters $a, b, c$. This far, reasonably standard linear regression. However, the twist is that the values $y$ are known exactly (essentially by definition, see the longer explanation in the end), while the $x$ -values have errors, which are not equal between points (heteroscedasticity is the word, I believe?). The errors originate from analog measurement noise, so the simple model is that they are independent and normally distributed. They're not equal because they're actually a ratio of the raw measurements, and by error propagation the error in the ratio depends on the values of the points themselves. The error magnitudes vary by quite a lot, so I'd like to do a weighted fit.

First, let me point out that I'm a theoretical physicist by education, which means that my last encounter with serious statistics was back in undergraduate courses, so I may be missing something obvious here.

So my questions are:

  1. is there some well known method to do the fit? Googling seems to lead to total least squares and generalized least squares, which deal with the case of errors in both dependent and independent variables, and I'd imagine my case would warrant some simplification vs. the fully general case
  2. my current idea to do this easily is to just use the fact that I know the approximate parameters, so I can trade the error in $x$ for an error in $y$ via error propagation/implicit function theorem (i.e. there is some $y$ which corresponds exactly to the actual unknown $x$, and I can estimate the difference from the known $y$ via the implicit function theorem). Since $b \approx 1$, this actually just makes the $y$ error approximately the same as the $x$ error. I could also iterate, i.e. use the initial fit to get a better estimate of the errors for a second round, and so on. However, the fact that I haven't seen this mentioned anywhere makes me suspect I'd be doing something horribly wrong?

Explanation of the underlying problem, in case there's an XY-problem here:

I'm measuring the gain of an analog circuit as a function of a control voltage applied, more specifically the digital control value used to set it via a DAC. I measure the actual gain $g$ as a function of the DAC value $n$ for a number of values $n$. So $n$ is exact, and $g$ has an error that I can estimate, since I know the noise floor of the system. Using that, I want to build a calibration function $n(g)$, which gives me the $n$ I need to set in order to achieve some desired gain $g$. Both from observation and theory, I know that the fit is almost linear, but with a small $g^2$ -term (due to Early voltage in a transistor, and the exact function is a actually a rational function, but within measurement accuracy I can't tell it apart from the $O(g^2)$ approximation).

So the natural thing, at least to me, is to tabulate the points as $(g, n)$ and then fit the polynomial, but this leads to the situation that the errors are in the independent variable. If this were a linear fit, I could just fit $(n, g)$ and then invert the function, but since I'm now fitting a polynomial, the inverse function will not itself be a polynomial, so I'd have to do a non-linear fit. Or I could invert approximately as a series, but that's piling one approximation on top of another, which I'd rather avoid.

A similar question that does not solve my problem:

Here's a similar question, but as the recommendation there is to use Deming regression, which however seems to give the same result as no errors at all when applied to the case where there is no error in the dependent variable.

  • $\begingroup$ I think Deming regression is what I would do if I had a good sense of what $\delta$ should be. It would be very unusual for the results from Deming regression to be the same as the case. I didn't try to prove it, but it seems like it can only happen if $\delta=0$ in the formula here. Can you post the $\delta$ and the other summary statistics for your problem? That formula is for the case of only linear term, yours has $x^2$ too. There are a few packages in R that fit general EIV model with unknown $\delta$, you could try eivreg. $\endgroup$
    – John L
    Feb 23 at 14:50
  • $\begingroup$ What I mean is that formula is for the simpler case with no quadratic term. I was asking if you can try that to verify you are doing it correctly and getting a different answer than the ordinary regression. Then, you can try the Deming regression with the quadratic polynomial you want. $\endgroup$
    – John L
    Feb 23 at 14:56
  • $\begingroup$ @JohnL sorry I was being unclear:$\delta$ is zero, since $\sigma_\epsilon=0$. It seems to me that then the error in the dependent variable disappears entirely from the formulas? $\endgroup$
    – Timo
    Feb 23 at 14:58
  • $\begingroup$ Yes, that means there is no error. $\endgroup$
    – John L
    Feb 23 at 14:59

You could try Heteroskedastic GP Regression methods which have an R implementation. In python you might find what you need using GPflow or BoTorch's HeteroskedasticSingleTaskGP.

I'll post an example if I have the time.

  • $\begingroup$ Cool, I'll check those out! The eventual implementation needs to be in, or interfaced to, C++. $\endgroup$
    – Timo
    Feb 23 at 9:49
  • $\begingroup$ You might learn a thing or two by looking at Binois' HetGP code here is the source code. $\endgroup$
    – ArnoV
    Feb 23 at 9:53
  • $\begingroup$ After a quick look, that looks like quite heavy artillery for what (perhaps naively) seems like still a reasonably simple problem to me. Do you have an opinion on my point 2., i.e. swapping the error to the $y$ -values using error propagation (the values of the parameters shouldn't vary by more than 10% or so from their ideal values) and then using weighted least squares? $\endgroup$
    – Timo
    Feb 24 at 7:50
  • $\begingroup$ I was perhaps too quick to answer. It seems your situation is not exactly heteroskedasticity as you want to model $y\approx f(x)$ but instead of having an error on the $y$ values $y = f(x) + \varepsilon$ you instead have an error $y = f(x + \varepsilon)$ on the inputs. What is the noise variance $\mathrm{Var}[\varepsilon]$ ? If it's small enough you should be able to use error propagation with a first of second order approximation. $\endgroup$
    – ArnoV
    Feb 24 at 8:04
  • $\begingroup$ Your problem reminds me of Post Non-Linear Modelling (PNL). This models a distortion $y = g(f(x) + \varepsilon)$, so you could use PNL related code to solve your problem by setting $f(x) = x$. $\endgroup$
    – ArnoV
    Feb 24 at 8:06

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.