0

I want to fit an ellipse model to a data. This sample data I extracted from a simple parametric plot of 2 Sin[t], Cos[t](Let's say, I don't know that) and I want to fit my model of an ellipse to this data. This is what I have tried:

Sample data:

data = Flatten[ Cases[ParametricPlot[{Cos[t], 2 Sin[t]}, {t, 0, 2 Pi}], Line[data_] :> data, Infinity], 1];

Model:

x^2/a^2 + y^2/b^2 == 1
f[x_]:=Sqrt[(1 - x^2/a^2) b^2];

and,

fit = NonlinearModelFit[data,f[x], {a,b}, x]

It throws me complex infinity error and I am not able to solve it. Any help is appreciated.

Rupesh
  • 887
  • 5
  • 10
  • 1
    You did not post your data so it is difficult to answer your question, but you should look at SingularValueDecomposition and at this answer 51549. – Tim Laska Aug 26 '20 at 23:48
  • @TimLaska Hi, data can be obtained from the first code data – Rupesh Aug 27 '20 at 00:34
  • Sorry, I missed that. – Tim Laska Aug 27 '20 at 00:35
  • The function form indicated will only be correct for the upper half of the ellipse. If you use the implicit form then there is no issue with negative y values. – Daniel Lichtblau Aug 27 '20 at 02:29
  • @DanielLichtblau yes, daniel it's for the positive half as I haven't put the other solution of y. However, Still throws me complexinfinity error without giving a fit. – Rupesh Aug 27 '20 at 02:31
  • 1
    You could use breg = BoundingRegion[data, "FastEllipse"] and get a quick and very close approximation to the ellipse which gives Ellipsoid[{-0.00227863, 0.00734437}, {{1.00538, 0.0000434228}, {0.0000434228, 4.07901}}] – flinty Aug 28 '20 at 00:54

2 Answers2

2

Since the model is of the form f(x,y) = 1, arrange the data as {{x,y,1} . . } for the fit:

(* make some data *)
(* the data is of the form { {x,y,1}, . . . } *)

eq = x^2/4 + y^2/9 == 1;


y[xx_] := y /. Solve[eq /. x -> xx, y]


points = Union[
   Flatten[Table[{x, y[x]}, {x, 
       Range[-2, 2, .1]}] /. {x_, {y1_, y2_}} -> {{x, y1, 1}, {x, y2, 
        1}}, 1]];


(* plot the data *)
ListPlot[points[[All, {1, 2}]]]


(* the general model *)


model = x^2/a^2 + y^2/b^2;


(* fit the data *)


fit = NonlinearModelFit[points, model, {a, b}, {x, y}];


fit["BestFitParameters"]


(* {a[Rule]2.,b\[Rule]2.999999999953799} *)
David Keith
  • 4,340
  • 1
  • 12
  • 28
  • it works for the sample file but If I am to do it for my data(which isn't too big) it gives me a fitting parameter of a circle that is not consistent with the data plot. If you have time could you please have a look at it, here is the file – Rupesh Aug 27 '20 at 02:15
  • @Rupesh Your data has exponents of -29 and -30. If you multiply the $x$ and $y$ values by 10^30, the code above works just fine. – JimB Aug 27 '20 at 02:51
  • @David Keith, thanks man you are awesome, worked like a charm – Rupesh Aug 27 '20 at 03:08
  • @JimB Thanks Jim, you are a life saver. – Rupesh Aug 27 '20 at 03:09
2

Here is a SingularValueDecomposition approach following @Danial Lichtblau's answer to 51549.

mean = Mean[data];
newpts = Map[# - mean &, data];
{uu, ww, vv} = SingularValueDecomposition[newpts, 2];
ListPlot[uu, AspectRatio -> Automatic]
rsqr = Mean[Map[#.# &, uu]];
{nx, ny} = Inverse[vv.ww].({x, y} - mean);
expr = Expand[nx^2 + ny^2] == rsqr;
expr = MultiplySides[expr, 1/expr[[2]]] // Expand
reg = ImplicitRegion[expr[[1]] <= expr[[2]] 1.1, {x, y}];
ContourPlot[Evaluate@expr, {x, y} \[Element] reg, 
 Epilog -> {Red, PointSize[Medium], Point[data]}, 
 AspectRatio -> Automatic]

SVD Images

Tim Laska
  • 16,346
  • 1
  • 34
  • 58
  • it's a nice approach and I appreciate you spending time on this, but it still doesn't fit to my main data. The solution given by David after using Jim's suggestion works. – Rupesh Aug 27 '20 at 03:19
  • 1
    @Rupesh I could not get your data file to import cleanly. If you edit your post to show how to import your file, then I could take a look. For example, you could Compress[data] and upload to pastebin. Then, you can regenerate the data by data2 = Uncompress[ Import["https://pastebin.com/raw/GmDR0nRH", "String"]]; data == data2 (* True *) as I did for your original dataset. – Tim Laska Aug 27 '20 at 04:29
  • I usually import these dat files I generate as TSV. Please try this. It should work.Import["data.dat", "TSV"] /. s_String :> ToExpression[s]. I will look option for Pastebin too, just woke up – Rupesh Aug 27 '20 at 13:11