Services & ResourcesWolfram ForumsMathGroup Archive

[Date Index] [Thread Index] [Author Index]

Re: show residuals of circle fit

• To: mathgroup at smc.vnet.net
• Subject: [mg79944] Re: [mg79900] show residuals of circle fit
• From: Daniel Lichtblau <danl at wolfram.com>
• Date: Thu, 9 Aug 2007 05:22:36 -0400 (EDT)
• References: <200708080853.EAA05894@smc.vnet.net>

will wrote:
> Dear Mathforum,
> 
> I am having a problem with displaying residuals after fitting a circle to
> a set of points. I am using some code that I found on the mathgroup 
> site (see:http://forums.wolfram.com/mathgroup/archive/2003/Jan/msg00019.html)
>  to generate the points and fit the circle:
> 
> points generation and display:
> 
> points = Table[theta = Random[];
>   {2*Cos[theta] + 1 + .05*Random[], 
>    2*Sin[theta] - 3 + .05*Random[]}, {30}]
> 
> plot = ListPlot[points, AspectRatio -> Automatic]
> 
> circle fitting and display with points:
> 
> fit circle:
> 
> hypersphereFit[data_List] := 
>  Module[{mat, rhs, qq, rr, params, cen}, 
>   mat = Map[Append[#, 1] &, data];
>   rhs = Map[#.# &, data];
>   {qq, rr} = QRDecomposition[mat];
>   params = LinearSolve[rr, qq.rhs];
>   cen = Drop[params, -1]/2;
>   {Sqrt[Last[params] + cen.cen], cen}]
> 
> display circle (using Graphics[Circle[]] inputing the results from the above
> function) and points:
> 
> Show[Graphics[
>   Circle[Part[hypersphereFit[points], 2], 
>    Part[hypersphereFit[points], 1], {0, 1}]], plot]
> 
> now i really want to be able to display the residuals for the points, but can't
> figure out how to. 
> 
> If it is any help i use the following code when displaying residuals from 
> a quadratic curve fit:
> 
> quadratic fit:
> 
> func = Fit[points, {1, x, x^2}, x]
> 
> define range of plot:
> 
> Insert[Part[Sort[points], {1, -1}, 1], x, 1]
> 
> out: {x, 2.1305, 3.04856}
> 
> copy and paste this output into the range for the Plot function:
> 
> Show[Graphics[
>   Plot[func, {x, 2.1305, 3.04856}]], 
>  ListPlot[points, PlotStyle -> {Red, PointSize[Small]}], Axes -> True,
>   AspectRatio -> Automatic]
> 
> this shows the quadratic fitted line with the points
> 
> then calculate and display the residuals for the quadratic fit:
> 
> res = Transpose[{First /@ points, 
>     Last /@ points - (func /. x -> First /@ points)}];
> ListPlot[res]
> 
> I would really, really like to be able to do this residuals calculation
> and display for the circle fit.
> 
> Any help would be much appreciated,
> 
> William Parr


I do not see a way to get quite the same notion of residual, or rather, 
I don't think it would be what you want. In your second example you 
"fit" to a function y=f(x) rather than to an implicit relation 
f(x,y)==0, and your residuals are y-discrepancies. For the circle fit 
you would presumably want total-least-squares error, that is, residuals 
that measure distance from computed center of circle, relative to the 
computed radius (say, negative if inside, positive if outside). I'm not 
sure what would be best to plot against as x-axis: actual x coordinate, 
or maybe just equally spaced abscissa values 1,2,...

Anyway, here is simple code to do something along the lines I propose. I 
rewrote slightly to use random generation capabilities of version 6 
Mathematica.

points = Table[theta = RandomReal[];
   {2*Cos[theta] + 1, 2*Sin[theta] - 3} + RandomReal[{.05}, {2}], {30}]

{r, c} = hypersphereFit[points]

res = Map[Norm[# - c] &, points] - r;

ListPlot[res]

If you really want to plot vs x coordinates, instead do:

res = Transpose[{First /@ points, Map[Norm[# - c] &, points] - r}]

ListPlot[res]


Daniel Lichtblau
Wolfram Research

• References:
  • show residuals of circle fit
    • From: will <willpowers69@hotmail.com>