Аппроксимация нелинейных отношений

Question

Это график рассеяния, который иллюстрирует взаимосвязь между двумя переменными:

Очевидно, что это нелинейная связь. Обе переменные являются временными рядами - точки - это разные наблюдения.

Как я могу подогнать кривую (в python), которая будет приближаться к ней?

РЕДАКТИРОВАТЬ:

Обратите внимание, что это не 2D-отношения (как указал JamesPhilips ниже).

Как я уже говорил, это два временных ряда. Я полагаю, что правильной вещью будет go для 3D-подгонки (включая время в третьем измерении). Таким образом, функция будет принимать два входа (х и время). Как это сделать?

EDIT2:

Я прилагаю образец этого набора данных здесь

EDIT3:

I мне повезло получить два высококачественных ответа от norok2 и JamesPhilips (большое спасибо им обоим!), и я буду изучать их. Однако у меня сложилось впечатление, что ни одна из идей, предложенных до сих пор, не использует в значительной степени тот факт, что это временные ряды. Моя интуиция заключается в том, что там есть какой-то сигнал (я знаю, что отсутствие временных отметок усложняет ситуацию). Так что я буду держать вопрос открытым некоторое время на случай, если кто-то захочет добавить другие идеи.

Кроме того, похоже, что набор данных, на который я поместил ссылку, был отсортирован не с использованием исходного индекса (мой плохой !) - Я помещаю ссылку на правильно отсортированный набор данных здесь .

Answer 1

В соответствии с вопросом, два столбца 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) предполагаются бесшумными, то простой 1D-интерполятор, как предусмотрено 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)), и простым, но эффективным подходом по снижению шума будет фильтрация Гаусса (как предусмотрено * 1054). *:

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')

Это поможет вам Интерполятор без необходимости явно подавлять шум на вашем целевом сигнале, который может быть более или менее желательным, в зависимости от приложения.

---

Повторная параметризация x(t') и y(t') с t'(t) (переупорядочение)

Если мы ослабим требование, что x(t) и y(t) относятся к временной серии, мы могли бы исследовать x(t') и y(t') для данного t'(t) преобразования.

Возможное преобразование, которое приводит к несколько интересному результату, получается путем сортировки данных CSV по y (временные ряды отсортированы по x):

data = data[data[:, 1].argsort()]
x = data[:, 0]
y = data[:, 1]

с этим преобразованием мы получаем следующий интерполятор для ANN подход:

и это для сглаженных x(t') и y(t'):

Возможно, существуют более эффективные переупорядочения, но их может быть непросто сформулировать. Относительно простая формулировка может включать кластеризацию, но я считаю, что этот ответ уже довольно длинный.

(полный код доступен здесь )

Answer 2

В комментариях приведен графический Python сборщик поверхности, считывающий данные из файла CSV. Вы должны быть в состоянии щелкнуть мышью, перетащить трехмерные уловки в 3-пространстве для проверки.

В этом примере я угадал простое уравнение плоской плоскости "csv_column_two = (a * index) + (b * csv_column_one) + c ", потому что трехмерная диаграмма рассеяния и трехмерная плоскость поверхности показывают, что могут быть выбросы с левой стороны в виде графика. С этим примером в руках вы можете легко попробовать вариации данных и уравнения. Установщик также печатает среднеквадратическое значение и значение 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)