4
$\begingroup$

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?

$\endgroup$
  • 2
    $\begingroup$ data is undefined. $\endgroup$ – JimB Jan 5 '18 at 17:55
  • $\begingroup$ 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. $\endgroup$ – Daniel Lichtblau Jan 5 '18 at 19:44
  • $\begingroup$ 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. $\endgroup$ – JimB Jan 9 '18 at 17:12
3
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ +1 Cool. What about the errors of the fitted parameters? $\endgroup$ – corey979 May 2 '18 at 9:58

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.