Как мне совместить две линейные функции с набором точек данных?

Question

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

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

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

Одна из идей, которые у меня были, - это выбрать «точку отсечения», скажем, вt=100 (ось X) и подгоните точки слева к линии, а точки справа к другой линии.Затем я мог бы вычислить chi ^ 2 для обеих строк, чтобы получить «доброту» подгонок.Тогда я мог бы просто продолжать делать одно и то же много раз с немного разными точками среза, пока не найду пару линий, которая дает наилучшие общие соответствия.

Другая идея, которая кажется более сложной, состоит в описанииэти точки данных представляют собой комбинацию двух линий, y= m1*t + m2*t + b1 + b2, где m s - наклоны, а b s - y-перехватчики.Тогда, взяв производную полной кривой, я бы получил dy/dt = m1+m2.Тогда, возможно, я мог бы циклически проходить через различные «точки отсечения» и подгоняться к линиям, пока не получу комбинацию, где производная наиболее близка к m1+m2.Но я не уверен, как это сделать, так как изначально я не работаю с одной кривой, а просто с группой дискретных точек данных.

Как лучше всего провести оптимизацию двух подгонок?в Python?

Answer 1

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

Сначала я пытаюсь воссоздать данные:

import numpy as np; import matplotlib.pyplot as plt
y1 = np.linspace(5.5, 3.7, num=100)
y1 = y1 + np.random.rand(y1.shape[0]) * np.linspace(0, .3, num=y1.shape[0])
y2 = np.linspace(3.7, 1.1, num=500)
y2 = y2 + np.random.rand(y2.shape[0]) * np.linspace(0.3, 1.9, num=y2.shape[0])
y = np.append(y1, y2)
x = np.array(range(len(y)))

ЗатемЯ делаю два линейных подгонки, используя numpy.linalg.lstsq, который сам основан на методе наименьших квадратов :

x1 = x[:100]
y1 = y[:100]
A1 = np.vstack([x1, np.ones(len(x1))]).T
m1, c1 = np.linalg.lstsq(A1, y1, rcond=None)[0]

x2 = x[100:]
y2 = y[100:]
A2 = np.vstack([x2, np.ones(len(x2))]).T
m2, c2 = np.linalg.lstsq(A2, y2, rcond=None)[0]

Построение этой диаграммы приводит к следующему изображению:

plt.scatter(x, y)
plt.plot(x1, m1 * x1 + c1, 'r')
plt.plot(x2, m2 * x2 + c2, 'r')
plt.show()

Теперь вы можете использовать алгоритм автоматической сегментации, например SWAB , для замены срезов [100:] и [:100], но я бы посоветовал против этого, если вы сможете вручную решить, в какой момент разбивать данные.Посмотрите на предоставленную ссылку, если вы ищете реализацию.

Answer 2

Ниже приведен пример подгонки двух прямых линий к одному набору данных, причем точка пересечения между двумя линиями также устанавливается в качестве параметра.В этом примере используется генетический алгоритм дифференциальной эволюции (DE) scipy для определения начальных оценок параметров.Реализация scipy в DE использует алгоритм Latin Hypercube для обеспечения тщательного поиска пространства параметров, и для этого алгоритма требуются границы, в которых производится поиск - в этом примере эти границы взяты из максимальных и минимальных значений данных.

import numpy, scipy, matplotlib
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from scipy.optimize import differential_evolution
import warnings

xData = numpy.array([19.1647, 18.0189, 16.9550, 15.7683, 14.7044, 13.6269, 12.6040, 11.4309, 10.2987, 9.23465, 8.18440, 7.89789, 7.62498, 7.36571, 7.01106, 6.71094, 6.46548, 6.27436, 6.16543, 6.05569, 5.91904, 5.78247, 5.53661, 4.85425, 4.29468, 3.74888, 3.16206, 2.58882, 1.93371, 1.52426, 1.14211, 0.719035, 0.377708, 0.0226971, -0.223181, -0.537231, -0.878491, -1.27484, -1.45266, -1.57583, -1.61717])
yData = numpy.array([0.644557, 0.641059, 0.637555, 0.634059, 0.634135, 0.631825, 0.631899, 0.627209, 0.622516, 0.617818, 0.616103, 0.613736, 0.610175, 0.606613, 0.605445, 0.603676, 0.604887, 0.600127, 0.604909, 0.588207, 0.581056, 0.576292, 0.566761, 0.555472, 0.545367, 0.538842, 0.529336, 0.518635, 0.506747, 0.499018, 0.491885, 0.484754, 0.475230, 0.464514, 0.454387, 0.444861, 0.437128, 0.415076, 0.401363, 0.390034, 0.378698])


def func(xArray, breakpoint, slopeA, offsetA, slopeB, offsetB):
    returnArray = []
    for x in xArray:
        if x < breakpoint:
            returnArray.append(slopeA * x + offsetA)
        else:
            returnArray.append(slopeB * x + offsetB)
    return returnArray


# function for genetic algorithm to minimize (sum of squared error)
def sumOfSquaredError(parameterTuple):
    warnings.filterwarnings("ignore") # do not print warnings by genetic algorithm
    val = func(xData, *parameterTuple)
    return numpy.sum((yData - val) ** 2.0)


def generate_Initial_Parameters():
    # min and max used for bounds
    maxX = max(xData)
    minX = min(xData)
    maxY = max(yData)
    minY = min(yData)
    slope = 10.0 * (maxY - minY) / (maxX - minX) # times 10 for safety margin

    parameterBounds = []
    parameterBounds.append([minX, maxX]) # search bounds for breakpoint
    parameterBounds.append([-slope, slope]) # search bounds for slopeA
    parameterBounds.append([minY, maxY]) # search bounds for offsetA
    parameterBounds.append([-slope, slope]) # search bounds for slopeB
    parameterBounds.append([minY, maxY]) # search bounds for offsetB


    result = differential_evolution(sumOfSquaredError, parameterBounds, seed=3)
    return result.x

# by default, differential_evolution completes by calling curve_fit() using parameter bounds
geneticParameters = generate_Initial_Parameters()

# call curve_fit without passing bounds from genetic algorithm
fittedParameters, pcov = curve_fit(func, xData, yData, geneticParameters)
print('Parameters:', fittedParameters)
print()

modelPredictions = func(xData, *fittedParameters) 

absError = modelPredictions - yData

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(yData))

print()
print('RMSE:', RMSE)
print('R-squared:', Rsquared)

print()


##########################################################
# graphics output section
def ModelAndScatterPlot(graphWidth, graphHeight):
    f = plt.figure(figsize=(graphWidth/100.0, graphHeight/100.0), dpi=100)
    axes = f.add_subplot(111)

    # first the raw data as a scatter plot
    axes.plot(xData, yData,  'D')

    # create data for the fitted equation plot
    xModel = numpy.linspace(min(xData), max(xData))
    yModel = func(xModel, *fittedParameters)

    # now the model as a line plot
    axes.plot(xModel, yModel)

    axes.set_xlabel('X Data') # X axis data label
    axes.set_ylabel('Y Data') # Y axis data label

    plt.show()
    plt.close('all') # clean up after using pyplot

graphWidth = 800
graphHeight = 600
ModelAndScatterPlot(graphWidth, graphHeight)

Answer 3

Формально, второй подход, который вы упомянули, - это попытка разместить ваши данные в конической части.По сути, пара прямых линий (POSL) является типом конического сечения , и, следовательно, может быть представлена ​​общим уравнением конической формы ax^2 + by^2 + 2hxy + 2gx + 2fy + c = 0 (вместо того, которое вы упомянули в своем вопросе, которое простов итоге получается одна прямая линия склона m1 + m2).Чтобы уточнить, , безусловно, будет уравнением приведенной выше формы, которое при построении дает POSL, который соответствует вашим данным.Может ли ваш алгоритм сходиться к нему - это отдельная история.

Вы можете попытаться использовать методы, чтобы найти коэффициенты a, b, h, g, f & c.В идеале, коэффициенты, которые вы получили бы для конического сечения, должны образовывать POSL, который бы точно подходил вашему набору данных.

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

Этот метод математически очень аккуратный, но я бы поспорил, что попытка заставить вашу обученную модель сходиться к POSL будет эквивалентна балансировке чего-либо на острие ножа (условиедля POSL очень узко, по сути).Скорее всего, в итоге получится уравнение параболы, эллипса или гиперболы (что может привести к тому, что вы в достаточной степени подберете свой набор данных, что сделает даже неоптимальное решение конической регрессии целесообразным).

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