CODEWITHSUNDEEP
pythonoptimizationleast-squares
lesterdisaster
lesterdisaster

Reputation: 11

How do I solve this problem using linear least squares with scipy?

I had the following problem on a test and my code didn't give me the exact results I needed, but I can't really find what went wrong.

We discovered a new comet whose elliptical orbit can be represented in a Cartesian $(x, y)$ coordinate system by the equation

$$ay^2+ bxy + cx + dy + e = x^2$$ Use a SciPy routine to solve the linear least squares problem to determine the orbital parameters a, b, c, >d, e, given the following observations of the comet’s position

Observation 1 2 3 4 5 6 7 8 9 10
x 1.02 0.95 0.87 0.77 0.67 0.56 0.44 0.3 0.16 0.01
y 0.39 0.32 0.27 0.22 0.18 0.15 0.13 0.12 0.13 0.15

I wrote this code, but when i plot the equation using matplotlib.contour the curve doesn't match the data points.

def fit(x, y):
    n = (np.shape(x)[0])
    A = np.array([y**2, x * y, x, y, np.ones(n)]).T
    b = x**2 

    return linalg.lstsq(A, b)[0]

obx = np.array([1.02, 0.95, 0.87, 0.77, 0.67, 0.56, 0.44, 0.3, 0.16, 0.01])
oby = np.array([0.39, 0.32, 0.27, 0.22, 0.18, 0.15, 0.13, 0.12, 0.13, 0.15])

fit(obx, oby)

Does somebody know what I am doing wrong here, should i maybe use curvefit instead of lstsq , or is my mistake in the plotting code?

Some follow-up clarification, the code I wrote gave this output for the constants a to e.

array([-2.63562548, 0.14364618, 0.55144696, 3.22294034, -0.43289427])

I plotted the result with this code

obx = np.array([1.02, 0.95, 0.87, 0.77, 0.67, 0.56, 0.44, 0.3, 0.16, 0.01])
oby = np.array([0.39, 0.32, 0.27, 0.22, 0.18, 0.15, 0.13, 0.12, 0.13, 0.15])

def data_plot(x, y, a, b, c, d, e):

    def f(x, y):
    
        return a * y**2 + b * x * y + c * x + d * y + e
    plt.close() 

    size = 100
    xrang = np.linspace(0,0.5, size)
    yrang = np.linspace(0,90, size)
    X, Y = np.meshgrid(xrang, yrang)

    F = f(X, Y)
    G = X**2

    plt.contour( (F-G), [0])

    plt.scatter(x, y)
    plt.xlim([-0.5, 1.5])
    plt.ylim([0, 0.5])

    plt.xlabel('x-coordinate')
    plt.ylabel('y-coordinate')
    plt.show()
    return None
data_plot(obx, oby, -2.63562548, 0.14364618, 0.55144696, 3.22294034, -0.43289427 
)

which give this, obviously wrong, result.plot

Upvotes: 0

Views: 92

Answers (1)

Nick ODell
Nick ODell

Reputation: 25364

Does somebody know what I am doing wrong here, should i maybe use curvefit instead of lstsq , or is my mistake in the plotting code?

I think it's a mistake in your plotting code. I plotted this in a different manner, and it agrees with the initial points.

from sympy import plot_implicit, symbols, Eq, And
from sympy.plotting.plot import List2DSeries
import numpy as np

obx = np.array([1.02, 0.95, 0.87, 0.77, 0.67, 0.56, 0.44, 0.3, 0.16, 0.01])
oby = np.array([0.39, 0.32, 0.27, 0.22, 0.18, 0.15, 0.13, 0.12, 0.13, 0.15])
x, y = symbols('x y')
a, b, c, d, e = -2.63562548,  0.14364618,  0.55144696,  3.22294034, -0.43289427
p1 = plot_implicit(Eq(a*y**2+ b*x*y + c*x + d*y + e, x**2), (x, -2, 2), (y, -2, 2), line_color='red')
p1.append(List2DSeries(obx, oby))
p1.show()

(Blue is initial points, red is the least-squares fit.)

ellipse plot

Upvotes: 0

Related Questions