5
$\begingroup$

I'm currently trying to use the Mathematica's FindFit command to fit a curve of some data to the bottom half of a circle. The data is of the form

Data = {{a1, a2}, {b1, b2}, ...} 

My code looks like as follows

s = FindFit[Data, y0 - Sqrt[R^2 - (x - x0)^2], {x0, y0, R}, x]

The window is printing the below error,

FindFit::nrlnum: The function value {74724.4 + (-1. + y)^2, 75066.6 + (-1. + y)^2, 75409.6 + (-1. + y)^2, 75753.3 + (-1. + y)^2, 76097.9 + (-1. + y)^2, 76443.2 + (-1. + y)^2, <<40>>, 91274.4 + (-1. + y)^2, 91652.5 + (-1. + y)^2, 92031.4 + (1. + y)^2, 92411.1 + (-1. + y)^2, <<71>>} is not a list of real numbers with dimensions {121} at {x0, y0, R} = {1., 1., 1.}. >>

Any help would be appreciated.

$\endgroup$
6
  • $\begingroup$ Welcome to the site, John! I have formatted your post for you. Your problem is a common syntax error. See mathematica.stackexchange.com/a/18395/8 - underscore is a reserved character. For this reason, I have marked your post as a duplicate of that question. $\endgroup$
    – Verbeia
    Nov 14, 2014 at 6:25
  • $\begingroup$ Well my actual script contains no underscores in the nomenclature. I am new to both mathematica as well as stackexchange so I guess I should have stayed as close to my script as possible. My data name is simple "Data" in my script so there is another underlying error present. I have removed this underscore in my question above to avoid this confusion. $\endgroup$
    – T-Ray
    Nov 14, 2014 at 15:12
  • 1
    $\begingroup$ There is an undefined variable y in either your data or in the actual code that you used, but did not show here. Please inspect your data, and/or quit the kernel and try again. $\endgroup$ Nov 14, 2014 at 15:52
  • 2
    $\begingroup$ Here are some relevant posts from MSE and other forums.1 2 3 4 $\endgroup$ Nov 14, 2014 at 17:24
  • $\begingroup$ There does appear to be an undefined variable within the code. Quitting the kernel appears to have fixed this problem however the window is now printing "FindFit::nrlnum: "The function value {-7.44087-273.373\ I,-7.37627-273.998\ I,-7.29269-274.623\ I,-7.21934-275.248\ I,<<44>>,-2.97508-303.373\ I,-2.93673-303.998\ I,<<71>>} is not a list of real numbers with dimensions {121} at {x0,y0,R} = {1.,1.,1.}."". $\endgroup$
    – T-Ray
    Nov 14, 2014 at 17:26

3 Answers 3

14
$\begingroup$

Here is an approach using linear regression:

circfit[pts_] := Module[{reg, lm, bf, exp, center, rad},
  reg = {2 #1, 2 #2, #2^2 + #1^2} & @@@ pts;
  lm = LinearModelFit[reg, {1, x, y}, {x, y}];
  bf = lm["BestFitParameters"];
  exp = (x - #2)^2 + (y - #3)^2 - #1 - #2^2 - #3^2 & @@ bf;
  {center, rad} = {{#2, #3}, Sqrt[#2^2 + #3^2 + #1]} & @@ bf;
  circlefit[{"expression" -> exp, "center" -> center, 
    "radius" -> rad}]]; circlefit[list_][field_] := field /. list;
circlefit[list_]["Properties"] := list /. Rule[field_, _] :> field;
circlefit /: ReplaceAll[fields_, circlefit[list_]] := fields /. list;
Format[circlefit[list_], StandardForm] := 
 HoldForm[circlefit]["<" <> ToString@Length@list <> ">"]

This assumes that underlying data is a circle and aim is to find center and radius. There are doubtless much better ways.

For fun (and done quickly, so apologies for some inefficiencies and ugliness):

Manipulate[
 pt1 = Table[{x, cy + Sqrt[rad^2 - (x - cx)^2]}, {x, -4, 4, 0.1}];
 pt2 = Table[{x, cy - Sqrt[rad^2 - (x - cx)^2]}, {x, -4, 4, 0.1}];
 pt = Cases[Join[pt1, pt2], {_Real, _Real}];
 pts = pt + RandomReal[{-p, p}, {Length[pt], 2}];
 With[{fit = circfit[pts]},
  Column[{Show[
     ContourPlot[Evaluate@fit["expression"], {x, -1, 8}, {y, -1, 8}, 
      Contours -> {0}, ContourShading -> None, 
      ContourStyle -> Directive[Red, Thick]], ListPlot[pts], 
     AspectRatio -> Automatic, PlotLabel -> fit["expression"] == 0, 
     ImageSize -> 400, PlotRange -> {{-1, 8}, {-1, 8}}],
    Grid[{{"variable", "value", "model"}, {"center", {cx, cy}, 
       fit["center"]}, {"radius", rad, fit["radius"]}}]
    }]], {p, 0.05, 1},
 {cx, 0, 2},
 {cy, 0, 2},
 {rad, 1, 5}
 ]

Animated Gif showing the result of the manipulate

$\endgroup$
5
  • $\begingroup$ I have to wonder what those much better ways might be? I'd be surprised if there is something better than a least squares solution, which is what you have provided. $\endgroup$ Nov 16, 2014 at 1:24
  • $\begingroup$ @DanielLichtblau...I guess I mean better ways of implementing the idea. I'm still not sure what the underlying problem is,e.g. whether estimation of known circle or circular approximation of closed figure...where perhaps interpolation better for latter... $\endgroup$
    – ubpdqn
    Nov 16, 2014 at 3:02
  • $\begingroup$ Thank you. I was able to get it working by simply Taylor expanding the lower half of the circle and fitting the curve to the produced polynomial. Your method appears to provide more correct values however. $\endgroup$
    – T-Ray
    Nov 17, 2014 at 17:37
  • 1
    $\begingroup$ Knowing the underlying problem as ubpdqn points out is critical. If all one needs is a circle that looks good visually, then the above method works fine (and that's not to take away from the above method). But if statistical inferences need to be made about the parameters, then a model that mimics the data generation process is needed. Then maximum likelihood or a Bayesian approach could be used to estimate the model parameters and associated measures of precision. $\endgroup$
    – JimB
    Jun 12, 2015 at 16:43
  • 1
    $\begingroup$ @DanielLichtblau you could use MLESAC or RANSAC if there is extreme noise present. That's still least squares/MLE but more robust to outliers, also once appeared in this image processing answer $\endgroup$
    – Histograms
    Jun 12, 2015 at 18:10
6
$\begingroup$

Not exactly the problem, but I have this code languishing on my hard drive to compute a circle center and radius given a bunch of pts. From there it is straight-forward to get the polynomial.

First make some data...

pts =
  Table[RandomReal[{.9, 1.1}] {Sin[#], Cos[#]} & @ RandomReal[{0, 2 π}], {10}]

(* {{0.347046, -0.930661}, {0.218395, -1.07027}, {0.890144, -0.179376}, \
{-0.828852, -0.367061}, {0.770728, 0.607237}, {-0.118834, 
  0.934153}, {0.21855, -0.949232}, {0.136237, 
  1.02154}, {-0.881469, -0.285087}, {-0.878687, -0.377516}}*)

Define a penalty function

penalty[c_, radius_, pts_] := (Map[radius - Norm[# - c] &, pts])^2 // Total

Use Minimize[ ] to find the center of the circle and radius

res = Minimize[penalty[{x, y}, rr, pts], {x, y, rr}]

(* {0.0268756, {x -> 0.0247297, y -> -0.0119147, rr -> 0.969958}} *)

circ = Circle[{x, y}, rr] /. (res // Last)  

Show[circ // Graphics, ListPlot[pts]]

enter image description here

$\endgroup$
2
  • $\begingroup$ I made a small adjustment to your code to eliminate the non-localized Symbol theta and replace it with an anonymous function. I hope you do not mind. $\endgroup$
    – Mr.Wizard
    Jun 6, 2017 at 0:38
  • $\begingroup$ Nope, thanks for cleaning it up. $\endgroup$
    – MikeY
    Jun 6, 2017 at 14:28
1
$\begingroup$

I added a third element to the data list to give a "z" value that is a radius of 1.

Data = {{a1, a2, 1}, {b1, b2, 1}, ...} 

Then, the FindFit function works well.

FindFit[data, ((x - a)^2 + (y - b)^2)/r^2, {a, b, r}, {x, y}]
$\endgroup$
1
  • $\begingroup$ for this to work reliably you need to supply a good initial guess for a,b (eg use the mean value of the data ) $\endgroup$
    – george2079
    Jun 5, 2017 at 16:54

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.