Nonlinear Weighted Total Least Squares

I'm having some trouble with algorithms for total nonlinear least squares weighted by errors in both x and y.

There is a solution by (Malengo and Pennecchi, 2013), that I am trying to recreate, but it doesn't seem to be correct. The example they fit is Pearson's data, which is a linear data set and so can be solved using York's method for comparison.

The dataset and current solution is:

d = {{0.0, 5.9}, {0.9, 5.4}, {1.8, 4.4}, {2.6, 4.6}, {3.3, 3.5}, {4.4,
3.7}, {5.2, 2.8}, {6.1, 2.8}, {6.5, 2.4}, {7.4, 1.5}};(*data*)
e = {{1000., 1.}, {1000., 1.8}, {500., 4.}, {800., 8.}, {200.,
20.}, {80., 20.}, {60., 70.}, {20., 70.}, {1.8, 100.}, {1.0,
500.}} // Sqrt[1/#] &;(*errors*)
Ux = DiagonalMatrix[e[[All, 1]]^2];(*x Covariance Matrix*)
Uy = DiagonalMatrix[e[[All, 2]]^2];(*y Covariance Matrix*)
f[r_] := a*r + b;(*Function*)
ps = Variables[Level[f[Nothing], {-1}]];(*parameter list*)
xs = Array[x, {Length[d]}];(*list of "x" values, located {xs,f[xs]}*)
dx = d[[All, 1]] - xs;(*x-minimization potential*)
dy = d[[All, 2]] - f[xs];(*y-minimization potential*)
fx = NMinimize[dx.Inverse[Ux].dx + dy.Inverse[Uy].dy,
Join[xs, ps]][[2, -2 ;;]](*Total potential*)

(*or equivalently*)

fx = NMinimize[
Total[(d[[All, 1]] - xs)^2/e[[All, 1]]^2 + (d[[All, 2]] - f[xs])^2/
e[[All, 2]]^2], Join[xs, ps]][[2, -2 ;;]]


Which results in:

{a -> 0.248787, b -> 1.63261}


Far from the actual values calculated via York potential:

NMinimize[
Total[(f[d[[All, 1]]] - d[[All, 2]])^2/(a^2 e[[All, 1]]^2 +
e[[All, 2]]^2)], ps,
Method -> "RandomSearch"][[2]](*York*)


Which has a result of:

{a -> -0.480533, b -> 5.47991}


In their paper it works fine, and I'm just not sure how my code is different from theirs. Can anyone spot where I'm going wrong?

• data is undefined. – JimB Jan 5 '18 at 17:55
• It seems to want a constraint in order to avoid a local min. fx = FindMinimum[{Total[(d[[All, 1]] - xs)^2/ e[[All, 1]]^2 + (d[[All, 2]] - f[xs])^2/e[[All, 2]]^2], b >= 4}, Join[xs, ps]] agrees with the Yorke computation. From what I can tell, those values for {a,b} will not lead to a minimization of the original objective though. – Daniel Lichtblau Jan 5 '18 at 19:44
• I can't access the referenced article without paying. Is the value of the objective function available from the article? Also, having 10 data points and 12 parameters seems a bit suspicious. Can their code be supplied? And maybe a reference to "York potential" would be helpful. I say that because solutions can be provided by folks unfamiliar with the jargon and subject matter. – JimB Jan 9 '18 at 17:12

So I managed to solve this issue some time ago, which was actually a simple case of incorrect starting values during optimization. A condensed form for the Nonlinear Weighted Total Least Squares (NWTLS) algorithm is below.

d = {{0.0, 5.9}, {0.9, 5.4}, {1.8, 4.4}, {2.6, 4.6}, {3.3, 3.5}, {4.4,3.7},
{5.2, 2.8}, {6.1, 2.8}, {6.5, 2.4}, {7.4, 1.5}};(*data*)

e = {{1000., 1.}, {1000., 1.8}, {500., 4.}, {800., 8.}, {200., 20.}, {80., 20.},
{60., 70.}, {20., 70.}, {1.8, 100.}, {1.0, 500.}} // Sqrt[1/#] &;(*errors*)

f[r_] := a*r + b;(*Function*)

NWTLS[{d_, e_, f_}] := Module[{ps, xs, pi, fx},
ps = Variables[Level[f[Nothing], {-1}]];
xs = Array[x, {Length[d]}];
pi = FindArgMin[Total[(f[d[[All, 1]]] - d[[All, 2]])^2], ps];
fx = FindArgMin[
Total[(d[[All, 1]] - xs)^2/e[[All, 1]]^2 + (d[[All, 2]] - f[xs])^2/e[[All, 2]]^2],
Join[{xs, d[[All, 1]]}\[Transpose], {ps, pi}\[Transpose]]];
{fx[[;; Length[xs]]], f[fx[[;; Length[xs]]]] /.
Rule @@@ ({ps, fx[[-Length[ps] ;;]]}\[Transpose])}\[Transpose];(*Optimized {x,y}*)
fx[[-Length[ps] ;;]]]

NWTLS[{d, e, f}]

Out= {-0.480533,5.47991}


The output now matches the York linear fitting method, which indicates that it is functioning properly.

• +1 Cool. What about the errors of the fitted parameters? – corey979 May 2 '18 at 9:58
• @ChaSta ditto for corey979's comment – Sterling Mar 26 at 21:56
• I use Bootstrapping to get parameter distributions, as nonlinear terms are typically not well approximated by a normal distribution, so are generally better described using distribution quantiles. I haven't looked if there is an analytical way to get normal error terms, but imagine that there is. – ChaSta Mar 31 at 8:55