Nonlinear Weighted Total Least Squares

Question

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?

Answer 1

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.