Экспоненциальная кривая распада, подходящая для numpy и scipy

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

В прошлом я делал это с numpy.linalg.lstsq для экспоненциальных функций и scipy.оптимизировать.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()

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

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

2 6
numpy scipy curve-fitting

2 ответа:

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

При использовании 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, I вам не нужно vstack ряд нулей, также работает следующее: 

l2 = np.linalg.lstsq(x[:, None], np.log(y))[0][0]
4

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

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