CODEWITHSUNDEEP
pythonnumpymatplotlibscipylinear-regression
Tom Kurushingal
Tom Kurushingal

Reputation: 6476

How to find the appropriate linear fit in Python?

I am trying to find the most appropriate linear fit for a large amount of data that has linear behaviour for most of samples. The data (link) when plotted in the raw form is as shown below:

enter image description here

I need the linear fit that encompasses most of the points as shown by the thick orange line in the figure below:

enter image description here

I tried calculating the mean of the points but how do I extract the linear region using Python?

Reproducible code

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


data = np.loadtxt('linear.dat', skiprows = 1, delimiter = '\t')
print data
x = data[:, 0]
y = data[:, 1:]
m = y.shape[0]
n = y.shape[1]


def linear_fit(x, a, b):
    return a * x + b

y_fit = np.empty(shape=(m, n))


for i in range(n):
    fit_y_fit_a, fit_y_fit_b = optimize.curve_fit(linear_fit, x, y[:, i])[0]
    y_fit[:, i] = fit_y_fit_a * x + fit_y_fit_b

y[~np.isfinite(y)] = 0
y_mean = np.mean(y, axis = 1)

fig = plt.figure(figsize=(5, 5))
fig.clf()
plot_y_vs_x = plt.subplot(111)
markers = itertools.cycle(('o', '^', 's', 'v', 'h', '>', 'p', '<'))
for i in range(n):
    plot_y_vs_x.plot(x, y, linestyle = '', marker = markers.next(), alpha = 1, zorder = 2)
    # plot_y_vs_x.plot(x, y_fit, linestyle = ':', color = 'gray', linewidth = 0.5, zorder = 1)
    plot_y_vs_x.plot(x, y_mean, linestyle = '-', linewidth = 3.0, color = 'red', zorder = 3)
plot_y_vs_x.set_ylim([-10, 10])
plot_y_vs_x.set_ylabel('Y', labelpad = 6)
plot_y_vs_x.set_xlabel('X', labelpad = 6)
fig.savefig('plot.pdf')
plt.close()

Upvotes: 0

Views: 6809

Answers (2)

Lee
Lee

Reputation: 31040

I'd suggest using polyfit and the poly1d class:

polyfit gives a least squares polynomial best fit, with an order of your choosing (1 for linear).

a=np.genfromtxt('linear.dat',skiprows=1)
x=a[:,0]
y=a[:,1:]
k = np.linspace(700,1000,50)
plt.clf()
for z in y.T:
    plt.scatter(x,z)
for z in y.T:
    fit = np.polyfit(x,z,2) # increase order to get better fit
    fit_fn = np.poly1d(fit)
    plt.plot(k,fit_fn(k))
plt.xlim(700,1000)

enter image description here

The data doesn't seem very linear to me. But you could just use the first 10 points: 

k = np.linspace(700,900,50)
plt.clf()
plt.scatter(x,y[:,5]) # e.g. line 5
fit = np.polyfit(x[-10:],y[-10:,5],1) # increase order to get better fit
fit_fn = np.poly1d(fit)
plt.plot(k,fit_fn(k))
plt.xlim(700,1000)

enter image description here

Upvotes: 1

kponz
kponz

Reputation: 528

You're looking to calculate the linear regression of your points. To do that, 

import numpy as np
x = np.array([0, 1, 2, 3])
y = np.array([-1, 0.2, 0.9, 2.1])
A = np.vstack([x, np.ones(len(x))]).T
m, c = np.linalg.lstsq(A, y)[0]

This will give you values m and c that fit to y = mx + c. Simply replace x and y here with your own values as numpy arrays.

http://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.lstsq.html

Upvotes: 3

Related Questions