1
$\begingroup$

I have noisy measurements of movement along a circle. I want to fit a circle to these measurements. I tried two methods, a straight forward moment fit, and then an ODR fit (described here. However even after a large number of measurements the fit isn't very good. I Googled for alternative methods, and I came across this Bayesian method (described here and here).

Now I am not the biggest expert in the field of Statistics. I have a few questions:

1)

How does this expression

$$ \int_0^{2\pi} exp\left(-\frac{([x_i-a-Rcos(\theta)]^2+[y_i-b-Rsin(\theta)]^2)}{2\sigma^2}\right)\left(\frac{1}{2\pi R}\right)\left(Rd\theta\right) $$

become

$$ \int_{x^2+y^2=r^2} exp\left(-[(x-a)^2+(y-b)^2] \right)ds $$ and then finally $$ \frac{\sum_{n\geq0}[(a^2+b^2)r^2]^n}{r(n!)^2exp(r^2+a^2+b^2)}? $$

The first one is straight forward, the numerator measures the distance from the observation to the model circle. I understand the two latter expressions might be mathematically correct, but where are the measurements accounted for? Is $r$ the distance from $a$ and $b$ the measurement? If so, what happened to $R$, the radius of the model circle?

2)

I implemented the first expression in the below script. I effectively just made a very crude implementation of the first expression. To get something to minimize I just negated the result from the f function.

When I run the script I get this result:

Optimization terminated successfully.
         Current function value: -0.000000
         Iterations: 1
         Function evaluations: 36
True     : (-0.75, 17.1357733407, 7.17507917513)
ODR      : (-1.39004510002, 17.6456897937, 7.84162189191)
Moment   : (-0.718137104062, 17.331971949, 7.79343311432)
Bayesian : (-1.36598860985, 17.6329433542, 7.8070304964)

The fit is very close to the ODR fit. Is that really correct? Have I done something wrong? If I use a different starting point I get this result:

Optimization terminated successfully.
         Current function value: -0.000000
         Iterations: 1
         Function evaluations: 37
True     : (-0.75, 17.1357733407, 7.17507917513)
ODR      : (-1.39004510002, 17.6456897937, 7.84162189191)
Moment   : (-0.718137104062, 17.331971949, 7.79343311432)
Bayesian : (3.58792896155, 3.58792896155, 3.58792896155)

Which is obviously completely off. I don't know if something is off with my implementation, but it doesn't seem very stable.

I'd appreciate if someone could explain and perhaps give some hints how to implement the method correctly.

Here is the script:

import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import matplotlib.patches as mp
import matplotlib.lines as ml

measured = [
        (-0.688178755078, 15.8677171478),
        (-0.151151943426, 13.4866066925),
        (3.78976457129, 6.56969802333),
        (1.96316792749, 17.1880280513),
        (7.05896388477, 8.23926931886),
        (5.49530775317, 11.7184387005),
        (0.0423668803563, 17.6102815224),
        (6.18894286542, 20.3315009667),
        (4.28791104009, 17.5909328728),
        (2.25742233902, 14.3001354281),
        (1.5785891502, 21.8207081174),
        (3.86293369662, 17.8221320887),
        (2.23104933998, 29.8240525695),
        (3.51753336186, 18.1992607751),
        (3.69531639237, 23.7874420849),
        (-3.97255840552, 18.9568227347),
        (-4.35533625142, 26.5937899877),
        (-2.33999672523, 26.109688906),
        (-9.90157938056, 25.671717866),
        (-9.66152382784, 25.1489873874),
        (-1.83066740498, 20.9947460319),
        (-7.51964394487, 21.3606860471),
        (-7.16634412064, 15.0653254758),
        (-5.50410583327, 13.4459567027),
        (-8.96098586367, 14.0002843465),
        (-7.23686894667, 11.4199134054),
        (-10.1421616823, 12.5158925921),
        (-2.75614546129, 7.23673165366),
        (-3.56473630629, 9.81999032758),
        (-3.66950296152, 10.9589958433),
        (3.48133825317, 9.50925034721),
        (0.902034680549, 10.0559011217),
        (13.2432293744, 7.32620811258),
        (-3.85694247624, 12.0363962945),
        (10.2398103913, 7.33147515435),
        (10.1970422769, 16.8407961967),
        (7.39563335973, 11.5658979851),
        (5.69723738256, 16.1450897403),
        (6.30985342424, 19.3627554915),
        (4.05577574122, 22.7832781348),
        (4.04862541862, 16.8379949569),
        (4.11302641793, 24.8862198091),
        (0.295819734192, 21.5039718691),
        (3.00328469289, 22.8067409838),
        (1.69354198287, 28.324801484),
        (-4.59239755288, 26.0854351314),
        (0.0672279344978, 23.8946303345),
        (-1.17798319268, 25.2672845126),
        (-5.24007207626, 28.631765062),
        (-1.66445200451, 25.6904112081),
        (-10.7664748545, 17.3418198766),
        (-8.10668456322, 16.4001835089),
        (-7.71037072225, 17.8466977569),
        (-7.02580514688, 18.3018455867),
        (-3.9871884538, 11.9744631576),
        (-4.72464581963, 8.44923721752),
        (-5.44475445811, 13.456987894),
        (-9.81924966883, 20.9620348789),
        (0.725059600915, 12.6813331179),
        (-0.987527307949, 11.9616763468)
]

true = [
        (0.0, 10.0),
        (1.4672214011, 10.3118675362),
        (2.83753958756, 10.9219725008),
        (4.05106507913, 11.8036503793),
        (5.05476098867, 12.9183676175),
        (5.80476098867, 14.2174057232),
        (6.26828648023, 15.6439904976),
        (6.42507917513, 17.1357733407),
        (6.26828648023, 18.6275561837),
        (5.80476098867, 20.0541409582),
        (5.05476098867, 21.3531790638),
        (4.05106507913, 22.4678963021),
        (2.83753958756, 23.3495741805),
        (1.4672214011, 23.9596791451),
        (-8.881784197e-16, 24.2715466813),
        (-1.5, 24.2715466813),
        (-2.9672214011, 23.9596791451),
        (-4.33753958756, 23.3495741805),
        (-5.55106507913, 22.4678963021),
        (-6.55476098867, 21.3531790638),
        (-7.30476098867, 20.0541409582),
        (-7.76828648023, 18.6275561837),
        (-7.92507917513, 17.1357733407),
        (-7.76828648023, 15.6439904976),
        (-7.30476098867, 14.2174057232),
        (-6.55476098867, 12.9183676175),
        (-5.55106507913, 11.8036503793),
        (-4.33753958756, 10.9219725008),
        (-2.9672214011, 10.3118675362),
        (-1.5, 10.0),
        (-4.10782519111e-14, 10.0),
        (1.4672214011, 10.3118675362),
        (2.83753958756, 10.9219725008),
        (4.05106507913, 11.8036503793),
        (5.05476098867, 12.9183676175),
        (5.80476098867, 14.2174057232),
        (6.26828648023, 15.6439904976),
        (6.42507917513, 17.1357733407),
        (6.26828648023, 18.6275561837),
        (5.80476098867, 20.0541409582),
        (5.05476098867, 21.3531790638),
        (4.05106507913, 22.4678963021),
        (2.83753958756, 23.3495741805),
        (1.4672214011, 23.9596791451),
        (-6.8167693712e-14, 24.2715466813),
        (-1.5, 24.2715466813),
        (-2.9672214011, 23.9596791451),
        (-4.33753958756, 23.3495741805),
        (-5.55106507913, 22.4678963021),
        (-6.55476098867, 21.3531790638),
        (-7.30476098867, 20.0541409582),
        (-7.76828648023, 18.6275561837),
        (-7.92507917513, 17.1357733407),
        (-7.76828648023, 15.6439904976),
        (-7.30476098867, 14.2174057232),
        (-6.55476098867, 12.9183676175),
        (-5.55106507913, 11.8036503793),
        (-4.33753958756, 10.9219725008),
        (-2.9672214011, 10.3118675362),
        (-1.5, 10.0),
]

def odr(xi, yi):

        from scipy import odr
        import scipy

        x = scipy.r_[xi]
        y = scipy.r_[yi]

        def calc_R(xc, yc):
                return scipy.sqrt((x-xc)**2 + (y-yc)**2)

        def f_3(beta, x):
            """ implicit definition of the circle """
            return (x[0]-beta[0])**2 + (x[1]-beta[1])**2 -beta[2]**2

        x_m = scipy.mean(x)
        y_m = scipy.mean(y)

        # initial guess for parameters
        R_m = calc_R(x_m, y_m).mean()
        beta0 = [ x_m, y_m, R_m]

        # for implicit function :
        #       data.x contains both coordinates of the points (data.x = [x, y])
        #       data.y is the dimensionality of the response
        lsc_data  = odr.Data(scipy.row_stack([x, y]), y=1)
        lsc_model = odr.Model(f_3, implicit=True)
        lsc_odr   = odr.ODR(lsc_data, lsc_model, beta0)
        lsc_out   = lsc_odr.run()

        xc_3, yc_3, R_3 = lsc_out.beta

        return (xc_3, yc_3), R_3

def moment(xi, yi):
        a, b = sum(xi) / len(xi), sum(yi) / len(yi)
        r = sum([np.sqrt((tx - a) ** 2 + (ty - b) ** 2) for tx, ty in zip(xi, yi)]) / len(xi)
        return (a, b), r

def bayesian(xi, yi, (ia, ib, ir)):
        from scipy.optimize import fmin_powell
        ts = np.linspace(0,2*np.pi,360)
        def f(args, xi, yi, ts):
                a = args[0]
                b = args[1]
                r = args[2]
                s = 1.
                for x, y in zip(xi, yi):
                        temp = 0.
                        for t in ts:
                                temp += np.exp(-((x - a - r*np.cos(t))**2 + (y - b - r*np.sin(t))**2))
                        s *= temp
                return -s
        res = fmin_powell(f, np.array([ia, ib, ir]), (xi, yi, ts))
        return (res[0], res[1]), res[2]

mxs, mys = zip(*measured)
txs, tys = zip(*true)

# True circle fit
(ta, tb), tr = odr(txs, tys)

# Measured circle fit odr

(mao, mbo), mro = odr(mxs, mys)

# Measured circle fit moment

(mam, mbm), mrm = moment(mxs, mys)

# Measured circle fit bayesian

(mab, mbb), mrb = bayesian(mxs, mys, (mam, mbm, mrm))
#(mab, mbb), mrb = bayesian(mxs, mys, (1., 1., 1.))

print "True     : ({}, {}, {})".format(ta, tb, tr)
print "ODR      : ({}, {}, {})".format(mao, mbo, mro)
print "Moment   : ({}, {}, {})".format(mam, mbm, mrm)
print "Bayesian : ({}, {}, {})".format(mab, mbb, mrb)

f, a = plt.subplots()
a.add_patch(mp.Circle((ta, tb), tr, fill=False, color='g'))
a.add_patch(mp.Circle((mao, mbo), mro, fill=False, color='r'))
a.add_patch(mp.Circle((mam, mbm), mrm, fill=False, color='b'))
a.add_patch(mp.Circle((mab, mbb), mrb, fill=False, color='y'))
a.plot(mxs, mys, '+')
a.plot(txs, tys, 'o', color='g')
plt.axis('scaled')
plt.show()
$\endgroup$
  • 6
    $\begingroup$ You basically seem to be asking people to study several pages of (dense) R code, compare it to many pages of Python code, and identify bugs or errors. That's really unfair to expect. But if you would explain your code in detail (in English) and describe the efforts you have made to test it and identify any differences, it might become answerable. The first part of your question is readily answered by identifying "$R$" with "$r$". $\endgroup$ – whuber Jan 20 '15 at 15:24
  • $\begingroup$ I think the move from the first integral to the second integral remains incorrect when $r=R$. It should be$$\int_{x^2+y^2=R^2}\exp\{-[(x_i-a-x)^2-(y_i-b-y)^2]/2\sigma^2]\}\text{d}x\text{d}y$$ $\endgroup$ – Xi'an Jan 20 '15 at 17:02
  • $\begingroup$ @whuber Well call me boorish, but this is stackexchange isn't it? People answers questions and get rewarded by points for good efforts. The paper I quoted is well referred to, I'd expect any person in the Bayesian know to help out and explain... (Besides, I still don't understand replace R with r in the latter two expressions in the first questions. Don't you lose the observations just using a, b, and R?) $\endgroup$ – c00kiemonster Jan 22 '15 at 13:20
  • $\begingroup$ Part (1) is certainly on topic on this site, but part (2) is not. Neither we--nor almost all other SE sites--are code review sites. I would like only to suggest that if you remove the second part you are likely to increase your chances of getting answers. You could replace the entire second part by the key phrase "... give some hints how to implement the method." $\endgroup$ – whuber Jan 22 '15 at 16:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.