What have I done wrong implementing this Bayesian method for fitting a circle to noisy data?

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()

• 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$". – whuber Jan 20 '15 at 15:24
• 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$$ – Xi'an Jan 20 '15 at 17:02