Non-linear fitting to an ellipse

In this example, we are given a noisy series of data points which we want to fit to an ellipse. The equation for an ellipse may be written as a nonlinear function of angle, $\theta$ ($0 \le \theta < 2\pi$), which depends on the parameters $a$ (the semi-major axis) and $e$ (the eccentricity): $$ r(\theta; a,e) = \frac{a(1-e^2)}{1-e\cos\theta}. $$ To fit a sequence of data points $(\theta, r)$ to this function, we first code it as a Python function taking two arguments: the independent variable, theta, and a tuple of the parameters, p = (a, e). The function we wish to minimise is the difference between this model function and the data, r, defined as the method residuals:

def f(theta, p):
    a, e = p
    return a * (1 - e**2)/(1 - e*np.cos(theta))

def residuals(p, r, theta):
    return r - f(theta, p)

We also need to give leastsq an initial guess for the fit parameters, say p0 = (1,0.5). The simplest call to fit the function would then pass to leastsq the objects residuals, p0 and args=(r, theta) (the additional arguments needed by the residuals function):

plsq = leastsq(residuals, p0, args=(r, theta))

If at all possible, however, it is better to also provide the Jacobian (the first derivative of the fit function with respect to the parameters to be fitted). Expressions for these are straightforward to calculate and implement:

\begin{align*} \frac{\partial f}{\partial a} &= \frac{(1-e^2)}{1-e\cos\theta},\\ \frac{\partial f}{\partial e} &= \frac{a(1-e^2)\cos\theta - 2ae(1-e\cos\theta)}{(1-e\cos\theta)^2}. \end{align*}

However, the function we wish to minimise is the residuals function, $r - f$ so we need the negatives of these derivatives. Here is the working code and the fit result.

import numpy as np
from scipy import optimize
import matplotlib.pyplot as plt

def f(theta, p):
    a, e = p
    return a * (1 - e**2)/(1 - e*np.cos(theta))

# The data to fit
theta = np.array([0.0000,0.4488,0.8976,1.3464,1.7952,2.2440,2.6928,
                  3.1416,3.5904,4.0392,4.4880,4.9368,5.3856,5.8344,6.2832])
r = np.array([4.6073, 2.8383, 1.0795, 0.8545, 0.5177, 0.3130, 0.0945, 0.4303,
              0.3165, 0.4654, 0.5159, 0.7807, 1.2683, 2.5384, 4.7271])

def residuals(p, r, theta):
    """ Return the observed - calculated residuals using f(theta, p). """
    return r - f(theta, p)

def jac(p, r, theta):
    """ Calculate and return the Jacobian of residuals. """
    a, e = p
    da = (1 - e**2)/(1 - e*np.cos(theta))
    de = (-2*a*e*(1-e*np.cos(theta)) + a*(1-e**2)*np.cos(theta))/(1 -
                                                        e*np.cos(theta))**2
    return -da,  -de

# Initial guesses for a, e
p0 = (1, 0.5)
plsq = optimize.leastsq(residuals, p0, Dfun=jac, args=(r, theta), col_deriv=True)
print(plsq)

plt.polar(theta, r, 'x')
theta_grid = np.linspace(0, 2*np.pi, 200)
plt.polar(theta_grid, f(theta_grid, plsq[0]), lw=2)
plt.show()