8
$\begingroup$

I have this set of points:

sTest={{265., -12.4317}, {269., -12.4839}, {270., -12.6383}, {274., \
-12.6947}, {278., -12.6996}, {279., -12.7755}, {300., -12.8652}, \
{309., -12.9276}, {313., -13.0571}, {314., -13.0811}, {315., \
-13.0985}, {316., -13.1501}, {322., -13.1908}, {323., -13.2129}, \
{324., -13.2311}, {325., -13.2757}, {326., -13.3849}, {332., \
-13.7497}, {335., -13.762}, {355., -13.7981}, {359., -13.9145}, \
{364., -14.0233}, {368., -14.0583}, {369., -14.1957}, {371., \
-14.2555}, {377., -14.4896}, {378., -14.5427}, {390., -14.5961}, \
{399., -14.6112}, {412., -14.7939}, {414., -14.805}, {415., \
-14.9441}, {423., -15.1411}, {424., -15.2414}, {445., -15.256}, \
{454., -15.3825}, {458., -15.4697}, {459., -15.6437}, {460., \
-15.6559}, {461., -15.7129}, {467., -15.8281}, {468., -15.9638}, \
{470., -15.987}, {489., -16.295}, {500., -16.3044}, {501., -16.3501}, \
{504., -16.3555}, {513., -16.4882}, {514., -16.7162}, {544., \
-16.7981}, {547., -16.9239}, {548., -16.9534}}

And I initially draw a line from the first to last point:

enter image description here

And notice there are some points above and below the line. I which to find the line such that all points are below the line and the sum of the vertical distance from each point to the line is a minimum. I can manually, find an approximate minimum using Manipulate:

enter image description here

which is the line $-7.52-0.017x$ giving a value of the sum of $5.77$. I would think I could use NMinimize to find the minimum line even if I have to set a limit on the values of $a$ and $b$ via:

sumF[a_, b_, buf_] := 
 Sum[(a + b buf[[i, 1]]) - buf[[i, 2]], {i, 1, Length@buf}]

checkBuf[a_, b_, buf_] := Module[{},
  VectorQ[Table[(a + b buf[[i, 1]]) - buf[[i, 2]],
    {i, 1, Length@buf}], Positive]
  ]

NMinimize[{sumF[a, b, sTest], -8 < a < -6 && -1 < b < 0 && 
   checkBuf[a, b, sTest] == True}, {a, b}]

where checkBuf returns True if all points are below the line. However, this gives the error:

NMinimize::nsol: There are no points that satisfy the constraints {False}.

I was wondering if this is not a problem NMinimize can solve or if it is, if I'm not coding it correctly and if so could someone help me do so?

Thanks for reading.

$\endgroup$
2
  • $\begingroup$ The appears to be a strong periodicity in the residuals: ListPlot[Transpose[{sTest[[All, 1]], LinearModelFit[sTest, x, x]["FitResiduals"]}], Joined -> True]. Is there a rationale for ignoring that? $\endgroup$
    – JimB
    Commented Oct 25, 2022 at 16:42
  • $\begingroup$ @JimB: I see what you mean. This data plots the log of the accuracy as a function of number of terms for a fractional power expansion $\sum_{n=1}^N c_i (z^{1/9})^n$. I assume the periodicity stems from the exponential terms. I'm interested in generating a simple linear expression $A(N)=a+bN$ which gives the minimum accuracy for $N$ terms (hense having all the points below the line. $\endgroup$
    – josh
    Commented Oct 25, 2022 at 17:35

4 Answers 4

9
+100
$\begingroup$

We may define a line: c0+ c1 x with parameters c0 and c1. To determine the parameters we request that the following sum is minimized:

Total[((c0  + c1 #[[1]]) - #[[2]]) & /@ sTest];

under the conditions that the line is above the points:

 And @@ (((c0  + c1 #[[1]]) >= #[[2]]) & /@ sTest)

With this we can determine the parameters c0 and c1:

line = c0 + c1 x /. 
  Minimize[{Total[((c0  + c1 #[[1]]) - #[[2]]) & /@ sTest], 
     And @@ (((c0  + c1 #[[1]]) >= #[[2]]) & /@ sTest)}, {c0, c1}][[2]]
Plot[line, {x, sTest[[1, 1]], sTest[[-1, 1]]}, Epilog -> Point[sTest]]

enter image description here

$\endgroup$
3
  • $\begingroup$ Outstanding Daniel! Let me check it a bit on some data sets and get back to you later. $\endgroup$
    – josh
    Commented Oct 25, 2022 at 22:03
  • $\begingroup$ Works well for my application over varied data sets. Also have some fitted to quadratics nicely! Thanks! $\endgroup$
    – josh
    Commented Oct 29, 2022 at 11:59
  • $\begingroup$ Glad to be of some help. $\endgroup$ Commented Oct 29, 2022 at 14:48
8
$\begingroup$

It seems plausible that a line that connecting two points on the convex hull is the solution.

(* Find points in convex hull *)
chr = ConvexHullRegion[sTest][[1]];

(* Find the two points on the convex hull border that satisfy the constraints and
minimize the total error *)
minError = ∞;
besta = 0;
bestb = 0;
(* Loop through all possible pairs of points on the convex hull *)
Do[Do[
  (* Slope *)
  bb = (chr[[i, 2]] - chr[[j, 2]])/(chr[[i, 1]] - chr[[j, 1]]);
  (* Intercept *)
  aa = chr[[i, 2]] - bb chr[[i, 1]];
  
  If[Min[aa + bb chr[[All, 1]] - chr[[All, 2]]] == 0,
   error = 
    aa Length[sTest] + bb Total[sTest[[All, 1]]] - 
     Total[sTest[[All, 2]]];
   If[minError == ∞, minError = error; besta = aa; 
    bestb = bb];
   If[minError > error, minError = error; besta = aa; bestb = bb]],
  
  {j, 1, i - 1}], {i, 2, Length[chr]}]

(* Show results *)
{besta, bestb}
(* {-7.8383, -0.0164702} *)

ListPlot[{sTest, {{Min[sTest[[All, 1]]], 
    besta + bestb Min[sTest[[All, 1]]]}, {Max[sTest[[All, 1]]], 
    besta + bestb Max[sTest[[All, 1]]]}}}, PlotStyle -> {Blue, Red}, 
 Joined -> {False, True}]

Data and line minimizing error that is "above" all points

Another alternative is to use NonlinearModelFit:

sTest = Rationalize[sTest, 0];
nlm = NonlinearModelFit[sTest, {a + b x, Min[a + b sTest[[All, 1]] - sTest[[All, 2]]] > 0},
  {{a, -8}, {b, -1/100}}, x, WorkingPrecision -> 20]
nlm["BestFitParameters"]
Show[ListPlot[sTest], Plot[nlm[x], {x, Min[sTest[[All, 1]]], Max[sTest[[All, 1]]]},
  PlotStyle -> Red]]

Data and line with smallest error all above the data points

$\endgroup$
1
  • $\begingroup$ Thanks Jim! I suspected a convex hull might be applicable. Will check it with some larger data sets and get back to you. $\endgroup$
    – josh
    Commented Oct 25, 2022 at 22:04
3
$\begingroup$

A slightly generalized supplement of @DanielHuber's nice answer, which describes the line in vector form {p1,p2}+t {e1,e2} (includes possible case of infinte slope)

e = {e1, e2};
p = {p1, p2};

This approach minimizes the projection perpendicular to the unknown direction e and is therefore independent of the underlying xy-coordinate system

\[CapitalDelta] = Map[Cross[e] . (# - p) &, sTest];
mini = NMinimize[{\[CapitalDelta] . \[CapitalDelta],Map[# < 0 &, \[CapitalDelta]], e . e == 1}, {e1, e2, p1, p2}]

Show[{ListPlot[sTest],Graphics[InfiniteLine[{p, p + e}] /. mini[[2]]]   }] 

enter image description here

$\endgroup$
2
$\begingroup$

A variation of Daniel Huber's answer is to use matrix-form LinearOptimization directly. Speed doesn't matter on the example case much, but this is 2.7 times faster on my system, and with ten million points, this grows to a tenfold improvement:

line = {1, x} . LinearOptimization[
    {Length[sTest], Total[sTest[[All, 1]]]},
    Transpose[{{1, #[[1]]}, -#[[2]]} & /@ sTest]];
Plot[line, {x, sTest[[1, 1]], sTest[[-1, 1]]}, Epilog -> Point[sTest]]
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.