5

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.

m_goldberg
  • 107,779
  • 16
  • 103
  • 257
T-Ray
  • 321
  • 2
  • 9
  • Welcome to the site, John! I have formatted your post for you. Your problem is a common syntax error. See http://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. – Verbeia Nov 14 '14 at 06:25
  • 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. – T-Ray Nov 14 '14 at 15:12
  • 1
    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. – Sjoerd C. de Vries Nov 14 '14 at 15:52
  • 2
    Here are some relevant posts from MSE and other forums.1 2 3 4 – Daniel Lichtblau Nov 14 '14 at 17:24
  • 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.}."". – T-Ray Nov 14 '14 at 17:26
  • It appears that x is too big (the first coordinates of your data) for the starting point {x0,y0,R} = {1.,1.,1.}. Try specifying a starting point: FindFit[.., {{x0,x1}, {y0,y1}, {R,R1}}, x], where x1, y1, and R1 are well-chosen numbers. For x1, I'd try x1 = Mean[Data[[All, 1]]]. You want R1 > Max[Abs[Data[[All, 1]] - x1]]. Or perhaps first try a constraint R - (x - x0)^2 > 0. I can't really check it out, without your data. – Michael E2 Nov 14 '14 at 19:40

3 Answers3

14

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

SPPearce
  • 5,653
  • 2
  • 18
  • 40
ubpdqn
  • 60,617
  • 3
  • 59
  • 148
  • 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. – Daniel Lichtblau Nov 16 '14 at 01:24
  • @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... – ubpdqn Nov 16 '14 at 03:02
  • 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. – T-Ray Nov 17 '14 at 17:37
  • 1
    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. – JimB Jun 12 '15 at 16:43
  • 1
    @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 – Histograms Jun 12 '15 at 18:10
6

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

Mr.Wizard
  • 271,378
  • 34
  • 587
  • 1,371
MikeY
  • 7,153
  • 18
  • 27
1

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}]
MarcoB
  • 67,153
  • 18
  • 91
  • 189
Tim Daughters
  • 11
  • 1
  • for this to work reliably you need to supply a good initial guess for a,b (eg use the mean value of the data ) – george2079 Jun 05 '17 at 16:54