近似非线性关系
Approximating a non-linear relationship
这是一个散点图，说明了两个变量之间的关系：
很明显这是一个非线性关系。两个变量都是时间序列 - 点是不同的观察结果。
如何拟合一条曲线（在 python 中）使其近似？
编辑：
请注意，这不是二维关系（正如 JamesPhilips 在下面指出的那样）。
正如我提到的，这是两个时间序列。我想正确的做法是进行 3D 拟合（包括作为三维的时间）。因此该函数需要两个输入（x 和时间）。怎么做？
编辑 2：
我附上了该数据集的样本 here
编辑 3：
我很幸运收到了 norok2 和 JamesPhilips 的两个高质量答案（非常感谢他们！），我将探索这些。然而，我的印象是 none 到目前为止提出的想法都充分利用了这些是时间序列这一事实。我的直觉是那里有一些信号（我知道，没有时间戳会使事情变得复杂）。所以我会把这个问题搁置一段时间，以防有人想提出其他想法。
此外，我放入 link 的数据集似乎没有使用原始索引进行排序（我的错！）- 我将 link 放入正确排序的数据集 here.
根据评论，这是一个图形 Python 表面装配器从 csv 文件中读取数据。您应该能够在 3-space 中单击鼠标并拖动并旋转 3D 策略以进行检查。
在这个例子中，我猜测了一个简单的平面方程 "csv_column_two = (a * index) + (b * csv_column_one) + c"，因为 3D 散点图和 3D 表面平面显示了绘制的左侧可能是异常值。有了这个例子，您可以轻松地尝试数据和方程式的变化。拟合器还打印 RMSE 和 R 平方值以帮助模型评估和比较。
import numpy, scipy, scipy.optimize
import matplotlib
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm # to colormap 3D surfaces from blue to red
import matplotlib.pyplot as plt
graphWidth = 800 # units are pixels
graphHeight = 600 # units are pixels
# 3D contour plot lines
numberOfContourLines = 16
def SurfacePlot(func, data, fittedParameters):
f = plt.figure(figsize=(graphWidth/100.0, graphHeight/100.0), dpi=100)
matplotlib.pyplot.grid(True)
axes = Axes3D(f)
x_data = data[0]
y_data = data[1]
z_data = data[2]
xModel = numpy.linspace(min(x_data), max(x_data), 20)
yModel = numpy.linspace(min(y_data), max(y_data), 20)
X, Y = numpy.meshgrid(xModel, yModel)
Z = func(numpy.array([X, Y]), *fittedParameters)
axes.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.coolwarm, linewidth=1, antialiased=True)
axes.scatter(x_data, y_data, z_data) # show data along with plotted surface
axes.set_title('Surface Plot (click-drag with mouse)') # add a title for surface plot
axes.set_xlabel('Index') # X axis data label
axes.set_ylabel('CSV file column 1') # Y axis data label
axes.set_zlabel('CSV file column 2') # Z axis data label
plt.show()
plt.close('all') # clean up after using pyplot or else there can be memory and process problems
def ContourPlot(func, data, fittedParameters):
f = plt.figure(figsize=(graphWidth/100.0, graphHeight/100.0), dpi=100)
axes = f.add_subplot(111)
x_data = data[0]
y_data = data[1]
z_data = data[2]
xModel = numpy.linspace(min(x_data), max(x_data), 20)
yModel = numpy.linspace(min(y_data), max(y_data), 20)
X, Y = numpy.meshgrid(xModel, yModel)
Z = func(numpy.array([X, Y]), *fittedParameters)
axes.plot(x_data, y_data, 'o')
axes.set_title('Contour Plot') # add a title for contour plot
axes.set_xlabel('Index') # X axis data label
axes.set_ylabel('CSV file column 1') # Y axis data label
CS = matplotlib.pyplot.contour(X, Y, Z, numberOfContourLines, colors='k')
matplotlib.pyplot.clabel(CS, inline=1, fontsize=10) # labels for contours
plt.show()
plt.close('all') # clean up after using pyplot or else there can be memory and process problems
def ScatterPlot(data):
f = plt.figure(figsize=(graphWidth/100.0, graphHeight/100.0), dpi=100)
matplotlib.pyplot.grid(True)
axes = Axes3D(f)
x_data = data[0]
y_data = data[1]
z_data = data[2]
axes.scatter(x_data, y_data, z_data)
axes.set_title('Scatter Plot (click-drag with mouse)')
axes.set_xlabel('Index')
axes.set_ylabel('CSV file column 1')
axes.set_zlabel('CSV file column 2')
plt.show()
plt.close('all') # clean up after using pyplot or else there can be memory and process problems
def func(data, a, b, c):
x = data[0]
y = data[1]
return (a * x) + (b * y) + c
if __name__ == "__main__":
filename = 'test_bfa_corr.csv'
filetext = open(filename, 'rt').read()
lines = filetext.split('\n')
xData = []
yData = []
zData = []
for i in range(len(lines)):
line = lines[i]
spl = line.split(',')
xData.append(i+1)
yData.append(spl[0])
zData.append(spl[1])
xData = numpy.array(xData, dtype=float)
yData = numpy.array(yData, dtype=float)
zData = numpy.array(zData, dtype=float)
data = [xData, yData, zData]
initialParameters = [1.0, 1.0, 1.0] # these are the same as scipy default values in this example
# here a non-linear surface fit is made with scipy's curve_fit()
fittedParameters, pcov = scipy.optimize.curve_fit(func, [xData, yData], zData, p0 = initialParameters)
ScatterPlot(data)
SurfacePlot(func, data, fittedParameters)
ContourPlot(func, data, fittedParameters)
print('fitted parameters', fittedParameters)
modelPredictions = func(data, *fittedParameters)
absError = modelPredictions - zData
SE = numpy.square(absError) # squared errors
MSE = numpy.mean(SE) # mean squared errors
RMSE = numpy.sqrt(MSE) # Root Mean Squared Error, RMSE
Rsquared = 1.0 - (numpy.var(absError) / numpy.var(zData))
print('RMSE:', RMSE)
print('R-squared:', Rsquared)
根据问题，x和y这两列是两个时间序列，即x(t)和y(t)。 t 参数由索引
表示
首先，让我们加载数据：
import io
import requests
import numpy as np
import scipy as sp
import matplotlib as mpl
import scipy.interpolate
import scipy.ndimage
import matplotlib.pyplot as plt
file_id = '1q4zY7B-BwG8bmbQJT3QvRt6B2MD4k0a0'
url = requests.get('https://drive.google.com/uc?export=download&id=' + file_id)
csv_file = io.StringIO(url.text)
data = np.loadtxt(csv_file, delimiter=',')
x = data[:, 0]
y = data[:, 1]
t = np.arange(len(x))
现在，y(x) 通常可能没有明确定义。通过绘制 x(t) 和 y(t)（可能沿着 y(x)）获得更有用的数据表示：
fig, ax = plt.subplots(1, 3, figsize=(15, 4), squeeze=False)
ax[0, 0].scatter(t, x, color='k', s=8.0)
ax[0, 1].scatter(t, y, color='k', s=8.0)
ax[0, 2].scatter(x, y, color='k', s=8.0)
ax[0, 0].plot(t, x, color='b')
ax[0, 1].plot(t, y, color='b')
ax[0, 2].plot(x, y, color='b')
请注意，虽然 y(x) 可视化得到两个聚类，即拉伸螺旋和直线，但没有更多信息，不应过度解释此观察结果。
现在，没有适合的模型，我们可以做的是为 x(t) 和 y(t) 设置一个插值数值函数。
如果假定 x(t) 和 y(t) 是无噪声的，那么一个简单的一维插值器，由 scipy.interpolate.interp1d():
提供
func_x_t = sp.interpolate.interp1d(t, x, kind='cubic', assume_sorted=True)
func_y_t = sp.interpolate.interp1d(t, y, kind='cubic', assume_sorted=True)
x_interp = func_x_t(t)
y_interp = func_y_t(t)
fig, ax = plt.subplots(1, 3, figsize=(15, 4), squeeze=False)
ax[0, 0].scatter(t, x, color='k', s=8.0)
ax[0, 1].scatter(t, y, color='k', s=8.0)
ax[0, 2].scatter(x, y, color='k', s=8.0)
ax[0, 0].plot(t, x_interp, color='r')
ax[0, 1].plot(t, y_interp, color='r')
ax[0, 2].plot(x_interp, y_interp, color='r')
请注意，红线现在由插值器生成。 SciPy 提供了多种不同的插值器，可能值得探索。
如果 x(t) 和 y(t) 是有噪声的测量值，则可以按上述方法获得更有用的插值器，但使用去噪的 x(t) 和 y(t)。在这里，我假设观察到的高频振荡是由噪声驱动的（在 x(t) 和 y(t) 中），一种简单但有效的去噪方法是高斯滤波（由 scipy.ndimage.gaussian_filter1d():
smooth_x = sp.ndimage.gaussian_filter1d(x, 12.0, mode='nearest')
smooth_y = sp.ndimage.gaussian_filter1d(y, 12.0, mode='nearest')
func_x_t = sp.interpolate.interp1d(t, smooth_x, kind='cubic', assume_sorted=True)
func_y_t = sp.interpolate.interp1d(t, smooth_y, kind='cubic', assume_sorted=True)
x_smooth_interp = func_x_t(t)
y_smooth_interp = func_y_t(t)
fig, ax = plt.subplots(1, 3, figsize=(15, 4), squeeze=False)
ax[0, 0].scatter(t, x, color='k', s=8.0)
ax[0, 1].scatter(t, y, color='k', s=8.0)
ax[0, 2].scatter(x, y, color='k', s=8.0)
ax[0, 0].plot(t, smooth_x, color='g')
ax[0, 1].plot(t, smooth_y, color='g')
ax[0, 2].plot(smooth_x, smooth_y, color='g')
ax[0, 0].plot(t, x_smooth_interp, color='r')
ax[0, 1].plot(t, y_smooth_interp, color='r')
ax[0, 2].plot(x_smooth_interp, y_smooth_interp, color='r')
请注意，*_smooth 和 *_smooth_interp 会叠加在一起。
另一种方法是使用人工神经网络，例如来自 scikit-learn:
import sklearn as skl
import sklearn.neural_network as skl_nn
import sklearn.preprocessing
x_train = t.reshape(-1, 1)
y_train = data
reg = skl_nn.MLPRegressor(
solver='adam', hidden_layer_sizes=(24, 8), activation='tanh',
learning_rate='adaptive', max_iter=1024)
reg.fit(x_train, y_train)
y_predict = reg.predict(x_train)
x_ann = y_predict[:, 0]
y_ann = y_predict[:, 1]
fig, ax = plt.subplots(1, 3, figsize=(15, 4), squeeze=False)
ax[0, 0].scatter(t, x, color='k', s=8.0)
ax[0, 1].scatter(t, y, color='k', s=8.0)
ax[0, 2].scatter(x, y, color='k', s=8.0)
ax[0, 0].plot(t, x, color='b')
ax[0, 1].plot(t, y, color='b')
ax[0, 2].plot(x, y, color='b')
ax[0, 0].plot(t, x_ann, color='r')
ax[0, 1].plot(t, y_ann, color='r')
ax[0, 2].plot(x_ann, y_ann, color='r')
这可以让您使用插值器，而无需明确地对目标信号进行降噪，这可能或多或少需要，具体取决于应用程序。
用 t'(t) 重新参数化 x(t') 和 y(t')（重新排序）
如果我们放宽 x(t) 和 y(t) 来自时间序列的要求，我们可以研究 x(t') 和 y(t') 对于给定的 t'(t) 转换.
通过按 y 对 CSV 数据进行排序（时间序列按 x 排序）：
data = data[data[:, 1].argsort()]
x = data[:, 0]
y = data[:, 1]
通过这种转换，我们获得了以下用于 ANN 方法的插值器：
这是平滑的 x(t') 和 y(t'):
可能有更有效的重新排序，但制定它们可能并不简单。
一个相对简单的公式可能涉及聚类，但我相信这个答案已经很长了。
（完整代码可用 here）
这是一个散点图，说明了两个变量之间的关系：
很明显这是一个非线性关系。两个变量都是时间序列 - 点是不同的观察结果。
如何拟合一条曲线（在 python 中）使其近似？
编辑：
请注意，这不是二维关系（正如 JamesPhilips 在下面指出的那样）。
正如我提到的，这是两个时间序列。我想正确的做法是进行 3D 拟合（包括作为三维的时间）。因此该函数需要两个输入（x 和时间）。怎么做？
编辑 2：
我附上了该数据集的样本 here
编辑 3：
我很幸运收到了 norok2 和 JamesPhilips 的两个高质量答案（非常感谢他们！），我将探索这些。然而，我的印象是 none 到目前为止提出的想法都充分利用了这些是时间序列这一事实。我的直觉是那里有一些信号（我知道，没有时间戳会使事情变得复杂）。所以我会把这个问题搁置一段时间，以防有人想提出其他想法。
此外，我放入 link 的数据集似乎没有使用原始索引进行排序（我的错！）- 我将 link 放入正确排序的数据集 here.
根据评论，这是一个图形 Python 表面装配器从 csv 文件中读取数据。您应该能够在 3-space 中单击鼠标并拖动并旋转 3D 策略以进行检查。
在这个例子中，我猜测了一个简单的平面方程 "csv_column_two = (a * index) + (b * csv_column_one) + c"，因为 3D 散点图和 3D 表面平面显示了绘制的左侧可能是异常值。有了这个例子，您可以轻松地尝试数据和方程式的变化。拟合器还打印 RMSE 和 R 平方值以帮助模型评估和比较。
import numpy, scipy, scipy.optimize
import matplotlib
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm # to colormap 3D surfaces from blue to red
import matplotlib.pyplot as plt
graphWidth = 800 # units are pixels
graphHeight = 600 # units are pixels
# 3D contour plot lines
numberOfContourLines = 16
def SurfacePlot(func, data, fittedParameters):
f = plt.figure(figsize=(graphWidth/100.0, graphHeight/100.0), dpi=100)
matplotlib.pyplot.grid(True)
axes = Axes3D(f)
x_data = data[0]
y_data = data[1]
z_data = data[2]
xModel = numpy.linspace(min(x_data), max(x_data), 20)
yModel = numpy.linspace(min(y_data), max(y_data), 20)
X, Y = numpy.meshgrid(xModel, yModel)
Z = func(numpy.array([X, Y]), *fittedParameters)
axes.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.coolwarm, linewidth=1, antialiased=True)
axes.scatter(x_data, y_data, z_data) # show data along with plotted surface
axes.set_title('Surface Plot (click-drag with mouse)') # add a title for surface plot
axes.set_xlabel('Index') # X axis data label
axes.set_ylabel('CSV file column 1') # Y axis data label
axes.set_zlabel('CSV file column 2') # Z axis data label
plt.show()
plt.close('all') # clean up after using pyplot or else there can be memory and process problems
def ContourPlot(func, data, fittedParameters):
f = plt.figure(figsize=(graphWidth/100.0, graphHeight/100.0), dpi=100)
axes = f.add_subplot(111)
x_data = data[0]
y_data = data[1]
z_data = data[2]
xModel = numpy.linspace(min(x_data), max(x_data), 20)
yModel = numpy.linspace(min(y_data), max(y_data), 20)
X, Y = numpy.meshgrid(xModel, yModel)
Z = func(numpy.array([X, Y]), *fittedParameters)
axes.plot(x_data, y_data, 'o')
axes.set_title('Contour Plot') # add a title for contour plot
axes.set_xlabel('Index') # X axis data label
axes.set_ylabel('CSV file column 1') # Y axis data label
CS = matplotlib.pyplot.contour(X, Y, Z, numberOfContourLines, colors='k')
matplotlib.pyplot.clabel(CS, inline=1, fontsize=10) # labels for contours
plt.show()
plt.close('all') # clean up after using pyplot or else there can be memory and process problems
def ScatterPlot(data):
f = plt.figure(figsize=(graphWidth/100.0, graphHeight/100.0), dpi=100)
matplotlib.pyplot.grid(True)
axes = Axes3D(f)
x_data = data[0]
y_data = data[1]
z_data = data[2]
axes.scatter(x_data, y_data, z_data)
axes.set_title('Scatter Plot (click-drag with mouse)')
axes.set_xlabel('Index')
axes.set_ylabel('CSV file column 1')
axes.set_zlabel('CSV file column 2')
plt.show()
plt.close('all') # clean up after using pyplot or else there can be memory and process problems
def func(data, a, b, c):
x = data[0]
y = data[1]
return (a * x) + (b * y) + c
if __name__ == "__main__":
filename = 'test_bfa_corr.csv'
filetext = open(filename, 'rt').read()
lines = filetext.split('\n')
xData = []
yData = []
zData = []
for i in range(len(lines)):
line = lines[i]
spl = line.split(',')
xData.append(i+1)
yData.append(spl[0])
zData.append(spl[1])
xData = numpy.array(xData, dtype=float)
yData = numpy.array(yData, dtype=float)
zData = numpy.array(zData, dtype=float)
data = [xData, yData, zData]
initialParameters = [1.0, 1.0, 1.0] # these are the same as scipy default values in this example
# here a non-linear surface fit is made with scipy's curve_fit()
fittedParameters, pcov = scipy.optimize.curve_fit(func, [xData, yData], zData, p0 = initialParameters)
ScatterPlot(data)
SurfacePlot(func, data, fittedParameters)
ContourPlot(func, data, fittedParameters)
print('fitted parameters', fittedParameters)
modelPredictions = func(data, *fittedParameters)
absError = modelPredictions - zData
SE = numpy.square(absError) # squared errors
MSE = numpy.mean(SE) # mean squared errors
RMSE = numpy.sqrt(MSE) # Root Mean Squared Error, RMSE
Rsquared = 1.0 - (numpy.var(absError) / numpy.var(zData))
print('RMSE:', RMSE)
print('R-squared:', Rsquared)
根据问题，x和y这两列是两个时间序列，即x(t)和y(t)。 t 参数由索引
首先，让我们加载数据：
import io
import requests
import numpy as np
import scipy as sp
import matplotlib as mpl
import scipy.interpolate
import scipy.ndimage
import matplotlib.pyplot as plt
file_id = '1q4zY7B-BwG8bmbQJT3QvRt6B2MD4k0a0'
url = requests.get('https://drive.google.com/uc?export=download&id=' + file_id)
csv_file = io.StringIO(url.text)
data = np.loadtxt(csv_file, delimiter=',')
x = data[:, 0]
y = data[:, 1]
t = np.arange(len(x))
现在，y(x) 通常可能没有明确定义。通过绘制 x(t) 和 y(t)（可能沿着 y(x)）获得更有用的数据表示：
fig, ax = plt.subplots(1, 3, figsize=(15, 4), squeeze=False)
ax[0, 0].scatter(t, x, color='k', s=8.0)
ax[0, 1].scatter(t, y, color='k', s=8.0)
ax[0, 2].scatter(x, y, color='k', s=8.0)
ax[0, 0].plot(t, x, color='b')
ax[0, 1].plot(t, y, color='b')
ax[0, 2].plot(x, y, color='b')
请注意，虽然 y(x) 可视化得到两个聚类，即拉伸螺旋和直线，但没有更多信息，不应过度解释此观察结果。
现在，没有适合的模型，我们可以做的是为 x(t) 和 y(t) 设置一个插值数值函数。
如果假定 x(t) 和 y(t) 是无噪声的，那么一个简单的一维插值器，由 scipy.interpolate.interp1d():
func_x_t = sp.interpolate.interp1d(t, x, kind='cubic', assume_sorted=True)
func_y_t = sp.interpolate.interp1d(t, y, kind='cubic', assume_sorted=True)
x_interp = func_x_t(t)
y_interp = func_y_t(t)
fig, ax = plt.subplots(1, 3, figsize=(15, 4), squeeze=False)
ax[0, 0].scatter(t, x, color='k', s=8.0)
ax[0, 1].scatter(t, y, color='k', s=8.0)
ax[0, 2].scatter(x, y, color='k', s=8.0)
ax[0, 0].plot(t, x_interp, color='r')
ax[0, 1].plot(t, y_interp, color='r')
ax[0, 2].plot(x_interp, y_interp, color='r')
请注意，红线现在由插值器生成。 SciPy 提供了多种不同的插值器，可能值得探索。
如果 x(t) 和 y(t) 是有噪声的测量值，则可以按上述方法获得更有用的插值器，但使用去噪的 x(t) 和 y(t)。在这里，我假设观察到的高频振荡是由噪声驱动的（在 x(t) 和 y(t) 中），一种简单但有效的去噪方法是高斯滤波（由 scipy.ndimage.gaussian_filter1d():
smooth_x = sp.ndimage.gaussian_filter1d(x, 12.0, mode='nearest')
smooth_y = sp.ndimage.gaussian_filter1d(y, 12.0, mode='nearest')
func_x_t = sp.interpolate.interp1d(t, smooth_x, kind='cubic', assume_sorted=True)
func_y_t = sp.interpolate.interp1d(t, smooth_y, kind='cubic', assume_sorted=True)
x_smooth_interp = func_x_t(t)
y_smooth_interp = func_y_t(t)
fig, ax = plt.subplots(1, 3, figsize=(15, 4), squeeze=False)
ax[0, 0].scatter(t, x, color='k', s=8.0)
ax[0, 1].scatter(t, y, color='k', s=8.0)
ax[0, 2].scatter(x, y, color='k', s=8.0)
ax[0, 0].plot(t, smooth_x, color='g')
ax[0, 1].plot(t, smooth_y, color='g')
ax[0, 2].plot(smooth_x, smooth_y, color='g')
ax[0, 0].plot(t, x_smooth_interp, color='r')
ax[0, 1].plot(t, y_smooth_interp, color='r')
ax[0, 2].plot(x_smooth_interp, y_smooth_interp, color='r')
请注意，*_smooth 和 *_smooth_interp 会叠加在一起。
另一种方法是使用人工神经网络，例如来自 scikit-learn:
import sklearn as skl
import sklearn.neural_network as skl_nn
import sklearn.preprocessing
x_train = t.reshape(-1, 1)
y_train = data
reg = skl_nn.MLPRegressor(
solver='adam', hidden_layer_sizes=(24, 8), activation='tanh',
learning_rate='adaptive', max_iter=1024)
reg.fit(x_train, y_train)
y_predict = reg.predict(x_train)
x_ann = y_predict[:, 0]
y_ann = y_predict[:, 1]
fig, ax = plt.subplots(1, 3, figsize=(15, 4), squeeze=False)
ax[0, 0].scatter(t, x, color='k', s=8.0)
ax[0, 1].scatter(t, y, color='k', s=8.0)
ax[0, 2].scatter(x, y, color='k', s=8.0)
ax[0, 0].plot(t, x, color='b')
ax[0, 1].plot(t, y, color='b')
ax[0, 2].plot(x, y, color='b')
ax[0, 0].plot(t, x_ann, color='r')
ax[0, 1].plot(t, y_ann, color='r')
ax[0, 2].plot(x_ann, y_ann, color='r')
这可以让您使用插值器，而无需明确地对目标信号进行降噪，这可能或多或少需要，具体取决于应用程序。
用 t'(t) 重新参数化 x(t') 和 y(t')（重新排序）
如果我们放宽 x(t) 和 y(t) 来自时间序列的要求，我们可以研究 x(t') 和 y(t') 对于给定的 t'(t) 转换.
通过按 y 对 CSV 数据进行排序（时间序列按 x 排序）：
data = data[data[:, 1].argsort()]
x = data[:, 0]
y = data[:, 1]
通过这种转换，我们获得了以下用于 ANN 方法的插值器：
这是平滑的 x(t') 和 y(t'):
可能有更有效的重新排序，但制定它们可能并不简单。 一个相对简单的公式可能涉及聚类，但我相信这个答案已经很长了。
（完整代码可用 here）