0
$\begingroup$

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.

$\endgroup$
6
  • 1
    $\begingroup$ 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. $\endgroup$
    – Tim Laska
    Aug 26, 2020 at 23:48
  • $\begingroup$ @TimLaska Hi, data can be obtained from the first code data $\endgroup$
    – Rupesh
    Aug 27, 2020 at 0:34
  • $\begingroup$ Sorry, I missed that. $\endgroup$
    – Tim Laska
    Aug 27, 2020 at 0:35
  • $\begingroup$ 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. $\endgroup$ Aug 27, 2020 at 2:29
  • 1
    $\begingroup$ 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}}] $\endgroup$
    – flinty
    Aug 28, 2020 at 0:54

2 Answers 2

1
$\begingroup$

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`} *)
$\endgroup$
4
  • $\begingroup$ 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 $\endgroup$
    – Rupesh
    Aug 27, 2020 at 2:15
  • $\begingroup$ @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. $\endgroup$
    – JimB
    Aug 27, 2020 at 2:51
  • $\begingroup$ @David Keith, thanks man you are awesome, worked like a charm $\endgroup$
    – Rupesh
    Aug 27, 2020 at 3:08
  • $\begingroup$ @JimB Thanks Jim, you are a life saver. $\endgroup$
    – Rupesh
    Aug 27, 2020 at 3:09
1
$\begingroup$

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

$\endgroup$
3
  • $\begingroup$ 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. $\endgroup$
    – Rupesh
    Aug 27, 2020 at 3:19
  • 1
    $\begingroup$ @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. $\endgroup$
    – Tim Laska
    Aug 27, 2020 at 4:29
  • $\begingroup$ 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 $\endgroup$
    – Rupesh
    Aug 27, 2020 at 13:11

Your Answer

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

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