I am trying to fit a quadratic curve across a scatter plot of two variables. Since both variables are noisy I cannot use an ordinary least square regression (OLS) and I would like to have a symmetrical answer regardless which one I choose to be the independent variable.
Total (orthogonal) least squares regression seem to be the answer.
Unfortunately it seems that the two implementations that I have tried to use so far give me very bad fits.
In the figure below the blue line is the OLS fit, which obviously could be improved. The red is the TLS fit using the matlab/octave code below which seems to be the standard approach using single-value decompositions (SVD). I copied it from the corresponding Wikipedia article here.
I can see how it is trying to fit the upper-right denser cloud of points to the (very large) expense of the long tail. I have tested the code with simulated data and it works so is does not seem to be a bug.
I have also tried this other matlab package that implements a different iterative approach but unless I give it very tight bounds to the model parameters also gives me "crazy" fits.
I want to be able to perform TLS between many variables with many bootstrapping repeats to obtain confidence intervals and I cannot afford to come out with bounds for all those repeats eyeballing each graph; that would also tamper with the confidence interval estimation I guess.
I reckon the problem is that the data is a bit too noisy but I can do little about it. You can find the data here.
Any ideas how to go around this issue? Perhaps some know a more robust TLS method with an easy to use implementation available? There seem to be many academic articles but little runnable code.
Many thanks in advance! V.
function B = tls(xdata,ydata) % xdata = nx3 matrix [x.^2 x ones(size(x))] % ydata = nx1 matrix m = length(ydata); %number of x,y data pairs X = [xdata]; Y = ydata; n = size(X,2); % n is the width of X (X is m by n) Z = [X Y]; % Z is X augmented with Y. [U S V] = svd(Z,0); % find the SVD of Z. VXY = V(1:n,1+n:end); % Take the block of V consisting of the first n rows and the n+1 to last column VYY = V(1+n:end,1+n:end); % Take the bottom-right block of V. B = -VXY/VYY; end