Python: Pinpointing the Linear Part of a Slope [closed]

Question

I have several plots that look like the following:

I am wondering what kind of methods there might be for finding the slope between approximately 5.5 and 8 for the x-axis. Where there are several plots like this, I am more so wondering if there is a way to automatically find the slope value.

Any suggestions?

I am thinking ployfit(), or a linear regression. The problem is that I am unsure of how to find the values automatically.

Answer 1

A generic way to find linear parts in data sets is to calculate the second derivative of the function, and see where it is (close to) zero. There are several things to consider on the way to the solution:

• How to calculate the second derivative of noisy data? One fast and simple method, that can easily be adapted to different noise levels, data set sizes and expected lengths of the linear patch, is to convolve the data with a convolution kernel equal to the second derivative of a Gaussian. The adjustable part is the width of the kernel.

• What does "close to zero" mean in your context? To answer this question, you have to experiment with your data.

• The results of this method could be used as an input to the chi^2-method described above, to identify candidate regions in the data set.

Here some source code that will get you started:

from matplotlib import pyplot as plt

import numpy as np

# create theoretical data
x_a = np.linspace(-8,0, 60)
y_a = np.sin(x_a)
x_b = np.linspace(0,4,30)[1:]
y_b = x_b[:]
x_c = np.linspace(4,6,15)[1:]
y_c = np.sin((x_c - 4)/4*np.pi)/np.pi*4. + 4
x_d = np.linspace(6,14,120)[1:]
y_d = np.zeros(len(x_d)) + 4 + (4/np.pi)

x = np.concatenate((x_a, x_b, x_c, x_d))
y = np.concatenate((y_a, y_b, y_c, y_d))


# make noisy data from theoretical data
y_n = y + np.random.normal(0, 0.27, len(x))

# create convolution kernel for calculating
# the smoothed second order derivative
smooth_width = 59
x1 = np.linspace(-3,3,smooth_width)
norm = np.sum(np.exp(-x1**2)) * (x1[1]-x1[0]) # ad hoc normalization
y1 = (4*x1**2 - 2) * np.exp(-x1**2) / smooth_width *8#norm*(x1[1]-x1[0])



# calculate second order deriv.
y_conv = np.convolve(y_n, y1, mode="same")

# plot data
plt.plot(x,y_conv, label = "second deriv")
plt.plot(x, y_n,"o", label = "noisy data")
plt.plot(x, y, label="theory")
plt.plot(x, x, "0.3", label = "linear data")
plt.hlines([0],-10, 20)
plt.axvspan(0,4, color="y", alpha=0.2)
plt.axvspan(6,14, color="y", alpha=0.2)
plt.axhspan(-1,1, color="b", alpha=0.2)
plt.vlines([0, 4, 6],-10, 10)
plt.xlim(-2.5,12)
plt.ylim(-2.5,6)
plt.legend(loc=0)
plt.show()

This is the result:

smooth_width is the width of the convolution kernel. In order to adjust the amount of noise, change the value 0.27 in random.normal to different values. And please note, that this method does not work well close to the border of the data space.

As you can see, the "close-to-zero" requirement for the second derivative (blue line) holds quite well for the yellow parts, where the data is linear.