Подгонка кривой экспоненциального распада в numpy и scipy

У меня возникли некоторые проблемы с подгонкой кривой к некоторым данным, но я не могу понять, в чем я ошибаюсь.

Раньше я делал это с помощью numpy.linalg.lstsq для экспоненциальных функций и scipy.optimize.curve_fit для сигмоидных функций. На этот раз я хотел создать сценарий, который позволил бы мне определять различные функции, определять параметры и проверять их соответствие данным. При этом я заметил, что Scipy leastsq и Numpy lstsq, похоже, дают разные ответы для одного и того же набора данных и одной и той же функции. Функция просто y = e^(l*x) и ограничена таким образом, что y=1 на x=0.

Линия тренда Excel согласуется с результатом Numpy lstsq, но, поскольку Scipy leastsq может выполнять любую функцию, было бы хорошо выяснить, в чем проблема.

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

## Sampled data
x = np.array([0, 14, 37, 975, 2013, 2095, 2147])
y = np.array([1.0, 0.764317544, 0.647136491, 0.070803763, 0.003630962,     0.001485394,     0.000495131])

# function
fp = lambda p, x: np.exp(p*x)

# error function
e = lambda p, x, y: (fp(p, x) - y)

# using scipy least squares
l1, s =  optimize.leastsq(e, -0.004, args=(x,y))
print l1
# [-0.0132281]


# using numpy least squares
l2 = np.linalg.lstsq(np.vstack([x, np.zeros(len(x))]).T,np.log(y))[0][0]
print l2
# -0.00313461628963 (same answer as Excel trend line)

# smooth x for plotting
x_ = np.arange(0, x[-1], 0.2)

plt.figure()
plt.plot(x, y, 'rx', x_, fp(l1, x_), 'b-', x_, fp(l2, x_), 'g-')
plt.show()

Редактировать - дополнительная информация

MWE выше включает небольшую выборку набора данных. При подборе фактических данных кривая scipy.optimize.curve_fit представляет R ^ 2, равное 0,82, а кривая numpy.linalg.lstsq, которая совпадает с рассчитанной по Excel, R ^ 2 составляет 0,41.


numpy scipy least-squares curve-fitting exponential
person StacyR    schedule 16.01.2013    source источник

Ответы (2)


arrow_upward
4
arrow_downward

Вы сводите к минимуму различные функции ошибок.

Когда вы используете numpy.linalg.lstsq, минимизируемая функция ошибок будет

np.sum((np.log(y) - p * x)**2)

а scipy.optimize.leastsq минимизирует функцию

np.sum((y - np.exp(p * x))**2)

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

Отдельно отметим, что я не могу проверить это прямо сейчас, но при использовании numpy.linalg.lstsq, мне не нужно vstack строку нулей, также работает следующее:

l2 = np.linalg.lstsq(x[:, None], np.log(y))[0][0]
person Jaime    schedule 16.01.2013
comment
Спасибо @Jaime - отличный ответ! К сожалению, мои познания в математике не так хороши; один пишет или ошибается [также см. правку выше], или они просто принципиально разные ...? Каковы последствия для других функций, например, если я хочу проверить соответствие сигмовидной или кривой Гомперца одним и тем же данным? - person StacyR; 17.01.2013
comment
@StacyR У меня нет знаний, чтобы правильно ответить на ваш вопрос, но я почти уверен, что подгонка экспоненты, как вы это сделали с np.linalg.lstsq, - это просто быстрый и грязный трюк, который не вычисляет ошибки должным образом. Здесь есть некоторые обсуждения (мне трудно следить): mathworld.wolfram.com/LeastSquaresFittingExponential.html Если вы не хотите глубоко погружаться в этот материал, я бы применил метод scipy для всего: он должен лучше подходить, и ваши результаты будут согласованными для всех функций. - person Jaime; 18.01.2013
comment
еще раз спасибо! Я провел еще несколько исследований по этому поводу и, как вы упомянули, обнаружил, что метод np.linalg.lstsq чрезмерно взвешивает y-ошибки при низких значениях x. Ссылка, которой вы поделились, и некоторые другие ресурсы, которые я нашел, позволили мне получить еще один аналитический метод (что делает его сложным, так это ограничение - все книги описывают метод для y = a e ^ b < / i> x, а не y = e ^ b * x), однако это также дает худшую аппроксимирующую кривую, чем итеративный scipy.optimize.leastsq. - person StacyR; 18.01.2013

arrow_upward
1
arrow_downward

Чтобы немного пояснить точку зрения Хайме, любое нелинейное преобразование данных приведет к другой функции ошибок и, следовательно, к другим решениям. Это приведет к различным доверительным интервалам для параметров подгонки. Таким образом, у вас есть три возможных критерия для принятия решения: какую ошибку вы хотите минимизировать, какие параметры вы хотите больше доверять, и, наконец, если вы используете подгонку для прогнозирования некоторого значения, какой метод дает меньше ошибок в интересующем прогнозируемое значение. Небольшая аналитическая игра в Excel предполагает, что различные виды шума в данных (например, если функция шума масштабирует амплитуду, влияет на постоянную времени или является аддитивной), приводят к различным вариантам решения.

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

person user3117404    schedule 19.12.2013