CODEWITHSUNDEEP
pythonpython-2.7matplotlibscipylinear-regression
curious_cosmo
curious_cosmo

Reputation: 1214

Calculating slope through discrete points in Python

I've graphed four separate points using matplotlib, and would like to find the slope of the best fit line through them. Alternatively, I'd like to graph these points as one line and find the slope if possible. I'm also trying to do this on a logarithmic scale. Here's what I have (the period calculations came from elsewhere in my code):

import matplotlib.pyplot as plt
# First orbit
x_0 = 10.0
a1 = x_0
T1 = len(range(1, 98))

# Second orbit
x_0 = 5.0
a2 = x_0
T2 = len(range(1, 63))

# Third orbit
x_0 = 7.0
a3 = x_0
T3 = len(range(1, 81))

# Fourth orbit
x_0 = 13.0
a4 = x_0
T4 = len(range(1, 138))

smaxis = [a1, a2, a3, a4]
T = [T1, T2, T3, T4]

# Plot period versus semi-major axis
for i in range(len(T)):
    plt.plot(T[i], smaxis[i], markersize=3, marker='o')
plt.xlabel('Period (T)')
plt.ylabel('Semimajor Axis (a)')
plt.xscale('log')
plt.yscale('log')
plt.title('Logarithmic scale of T vs a')

And my graph is shown here.

I've tried using using linregress with this code:

from scipy.stats import linregress
linregress(T, smaxis)

But I'm not convinced that's correct, since T and smaxis are lists, and I need the slope between the line of best fit through the discrete points shown. How can I do this?

Upvotes: 0

Views: 7822

Answers (2)

harpan
harpan

Reputation: 8631

Consider below code which uses numpy's polyfit. 

x=T
y=smaxis
fit = np.polyfit(x, y, 1)
fit_fn = np.poly1d(fit)
s,i = fit
print("slope: ",s," intercept: ",i)
for i in range(len(T)):
    plt.plot(T[i], smaxis[i], markersize=3, marker='o')
plt.xlabel('Period (T)')
plt.ylabel('Semimajor Axis (a)')
plt.xscale('log')
plt.yscale('log')
plt.title('Logarithmic scale of T vs a')
plt.plot(x, fit_fn(x))
plt.show()

Output:

enter image description here

Upvotes: 3

Bill
Bill

Reputation: 11613

Here's how you catch and use the output of linregress.

from scipy.stats import linregress

slope, intercept, r_value, p_value, std_err = linregress(T, smaxis)

def a_predict(T):
    return intercept + slope*T

T_min, T_max = min(T), max(T)
a_min, a_max = a_predict(T_min), a_predict(T_max)

plt.plot([T_min, T_max], [a_min, a_max], 'r--')

print(slope, intercept, r_value, p_value, std_err)

output:

0.10753736192332683 -1.3585120207927215 0.9841584242334624 0.015841575766537552 0.013698301731763748

plot image

(I got this from the documentation).

But it may be more convenient for you to convert your lists into numpy arrays first.

import numpy as np
x = np.array(T)

Then you can do vectorized calculations as in the example in the documentation:

plt.plot(x, intercept + slope*x, 'r--')

Upvotes: 3

Related Questions