1
$\begingroup$

I'm trying to come up with a fit that will give me x,y coordinates of a point, that best fit in 2 slopes, the 2 slopes intercepting at 90°.

At the moment I fit individually 2 lines to a 2 sets of measured points.

topline = LinearModelFit[{{666.638, 769.106}, {530.36, 566.72}, {499.235, 519.357}, {530.6, 
  567.42}}, x, x];
botline = LinearModelFit[{{791.76, 123.04}, {591.025, 259.925}, {545.76, 292.54}, {586.919, 
  263.629}}, x, x];
{x0 = x /. Solve[Normal[topline] == Normal[botline], x][[ 1]],
  z0 = Normal[botline[x0]]}

This trivially gives x0 and y0 = 409.01, 385.807 but does not forces the 2 lines to be perpendicular. Because of the measurement errors in the set of coordinates, I routinely end with values in the 89-91° range. However it has to be 90°, because it physically set up that way. I would like to find a fit function that will come up with the best x0,y0 obeying a right angle of the 2 slopes.

I tried using Solve, itself containing the 2 above LinearModelFit with an additional variable point, but this does not work, as LinearModelFit won't work with undefined points.

$\endgroup$
1
  • $\begingroup$ I suggest that you think about what your errors are. Least squares fitting (e.g. y=a x + c) assumes that you have errors in your y values, but not your x values. If you might have solutions close to y=c, then this model seems questionable. $\endgroup$
    – mikado
    Jan 15 at 19:11

1 Answer 1

1
$\begingroup$

If I correctly understand you try to find two perpendicular lines fitting your data:

That means first set of data is fitted by a1+m1 x and second set is fitted by a2-1/m1 x!

ptop = {{666.638, 769.106}, {530.36, 566.72}, {499.235,519.357}, {530.6,567.42}};
pbot = {{791.76, 123.04}, {591.025, 259.925}, {545.76,292.54}, {586.919,263.629}};

mini = NMinimize[
{Total[Map[(#[[2]] - (a1 + m1 #[[1]]))^2 &, ptop]] + 
Total[Map[(#[[2]] - (a2 - 1/m1 #[[1]]))^2 &, pbot]], m1 > 0}, {a1,a2, m1}]
(* {8.01103, {a1 -> -217.216, a2 -> 660.241, m1 -> 1.47809}}*)

intersection point

solx = Solve[a1 + m1 x == a2 - 1/m1 x  , x][[1]]
pt = {x, a1 + m1 x} /. solx /. mini[[2]] (*{407.241, 384.724}*)

result

Show[{ListPlot[{ptop, pbot}], 
Plot[Evaluate[{a1 + m1 x, a2 - 1/m1 x} /. mini[[2]]], {x, 100,800}],   
Graphics[{Red, PointSize[Large],Point[pt]}]}, AspectRatio -> 1 ,PlotRange -> {{ 0, 1000}, {0, 1000}}]

enter image description here

remark

I know the resourcefunction ResourceFunction["MultiNonlinearModelFit"] which is able to handle two datasets but I failed trying ResourceFunction["MultiNonlinearModelFit"][{ptop, pbot}, {a1 + m1 x, a2 - 1/m1 x}, {a1, a2, m1}, x]

addendum

An easy solution using Region-functionality

define a cross-region

reg = RegionUnion[InfiniteLine[{x0, y0}, {Cos[alfa], Sin[alfa]}], 
InfiniteLine[{x0, y0}, {Cos[alfa + 90 \[Degree]],Sin[alfa + 90 \[Degree]]}]]

mini = FindMinimum[Total@Map[RegionDistance[reg , #]^2 &, Join[ptop, pbot]], 
{{x0, 400}, {y0, 500}, {alfa, 5 \[Degree]}}] // Quiet
(*{4.0295, {x0 -> 407.106, y0 -> 385.437, alfa -> 0.974064}}*)
$\endgroup$
7
  • $\begingroup$ @SjoerdSmit Any idea what's wrong with my attempt ResourceFunction["MultiNonlinearModelFit"][{ptop, pbot}, {a1 + m1 x, a2 - 1/m1 x}, {a1, a2, m1}, x]? Mathematica v12.2 displays "Finding Resource..." but doesn't evaluate. Thanks! $\endgroup$ Jan 15 at 18:37
  • $\begingroup$ I have a similar problem with v12.2.0 on Win7-x64. $\endgroup$
    – Syed
    Jan 15 at 18:55
  • $\begingroup$ Works in version 14 $\endgroup$ Jan 16 at 8:51
  • $\begingroup$ @DanielHuber Thanks for your reply! Hope it gives the same result $\endgroup$ Jan 16 at 8:57
  • $\begingroup$ @Ulrich Neumann Yes it does. $\endgroup$ Jan 16 at 10:07

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.