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. $\endgroup$ – JimB Jan 5 '18 at 17:55fx = 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