39.3: Ejemplo LSF - Estimando las mejores elipses

%matplotlib inline
import matplotlib.pylab as plt
import numpy as np
import sympy as sym
sym.init_printing(use_unicode=True)

Consideremos ahora los siguientes conjuntos de puntos. Piense en estos como observaciones del planeta moviéndose en una órbita elíptica.

x=[0, 1.0, 1.1, -1.1, -1.2, 1.3]
y =[2*1.5, 2*1.0, 2*-0.99, 2*-1.02, 2*1.2, 2*0]

plt.scatter(x,y)
plt.axis('equal')

En este problema queremos intentar ajustar una elipse a los datos anteriores. Primero veamos una ecuación general para una elipse:

\[ \frac{(u+x)^2}{a^2} + \frac{(v+y)^2}{b^2} = 1 \]

Dónde\(u\) y\(v\) son las\(y\) coordenadas\(x\) y para el centro de la elipse y\(a\) y\(b\) son las longitudes de los tamaños de los ejes de la elipse. Una búsqueda rápida sobre cómo trazar una elipse en python surge con el siguiente ejemplo:

# Code from: https://stackoverflow.com/questions/10952060/plot-ellipse-with-matplotlib-pyplot-python

u=1.     #x-position of the center
v=0.5    #y-position of the center
a=2.     #radius on the x-axis
b=1.5    #radius on the y-axis

t = np.linspace(0, 2*np.pi, 100)
plt.plot( u+a*np.cos(t) , v+b*np.sin(t) )
plt.grid(color='lightgray',linestyle='--')
plt.show()

Observe que este ejemplo utiliza ecuaciones de la forma:

\[t = [0, \dots, 2\pi] \nonumber \]

\[x = u+a\cos(t) \nonumber \]

\[y = v+b\sin(t) \nonumber \]

Resulta que esta forma de la ecuación es más fácil de trazar y las variables\(u\),\(v\),\(a\),\(b\) son las mismas que nuestra ecuación original.

Ahora vamos a expandir la ecuación original (ecuación\(\PageIndex{1}\) desde arriba) y obtenemos lo siguiente:

\[x^2−2ux-u^2+y^2−2vy+v^2=r^2 \]

Pregunta

¿Por qué no podemos\(\PageIndex{2}\) convertir la ecuación en la forma\(Ax=b\) y resolverla usando Ajuste de Mínimos Cuadrados? Discuta con tu grupo y prepárate para compartir tu pensamiento con la clase.

Si miramos nuestros datos más de cerca podemos simplificar las ecuaciones\(\PageIndex{1}\) y\(\PageIndex{2}\) asumiendo que el centroide\((u,v)\) está en el origen. Esta suposición da como resultado la siguiente ecuación:

\[ \frac{x^2}{a^2} + \frac{y^2}{b^2} = 1 \]

Observe que podemos reescribir esto en un conjunto lineal estándar de ecuaciones definiendo\(c_o = \frac{1}{a^2}\)\(c_1 = \frac{1}{b^2}\) y reescribiendo la ecuación de la siguiente manera:

\[ c_0x^2 + c_1y^2 = 1 \]

Hacer esto

Dado que conocemos los\(y\) valores\(x\) y de nuestras observaciones puntuales, la ecuación ahora\(\PageIndex{4}\) es lineal y se puede resolver usando Ajuste de mínimos cuadrados. Usar los puntos de observación de los constructos anteriores A y b como matrices numpy para el sistema sobredefinido\(Ax=b\)

sym.Matrix(A)

sym.Matrix(b)

Hacer esto

Resolver lo anterior sobre definido sistema de ecuaciones lineales para\(c_0\) y\(c_1\) usando LSF.

# Put your answer to the above question here

Suponiendo que tenemos\(c\) en el formato correcto, ahora podemos calcular\(a\) y\(b\) a partir de la solución para\(c_0\) y\(c_1\) calculado en el paso anterior y graficar usando nuestro código de trazado:

c = 1/np.sqrt(np.abs(c))
b=c[1,0] 
a=c[0,0]
print(a,b)

u=0     #x-position of the center
v=0     #y-position of the center

t = np.linspace(0, 2*np.pi, 100)
plt.plot(u+a*np.cos(t) , v+b*np.sin(t) )
plt.scatter(x,y)
plt.grid(color='lightgray',linestyle='--')
plt.axis('equal');