Я думаю, вы можете использовать функцию curve_fit из scipy.optimize ( документация ). Базовый пример использования:
import numpy as np
from scipy.optimize import curve_fit
def func(x, a, b, c):
return a*x**2 + b*x + c
x = np.linspace(0,4,50)
y = func(x, 5, 3, 4)
yn = y + 0.2*np.random.normal(size=len(x))
popt, pcov = curve_fit(func, x, yn)
Следуя документации, pcov дает:
Расчетная ковариация поп. Диагонали обеспечивают дисперсию
оценки параметра.
Таким образом, вы можете рассчитать оценку погрешности коэффициентов. Чтобы получить стандартное отклонение, вы можете взять квадратный корень из дисперсии.
Теперь у вас есть ошибка в коэффициентах, но она основана только на отклонении между данными и подгонкой. Если вы также хотите учесть ошибку в самих данных, функция curve_fit предоставляет аргумент sigma:
сигма: нет или N-длина последовательности
Если нет, то это стандартное отклонение ydata. это
вектор, если он задан, будет использоваться как вес в наименьших квадратах
проблема.
Полный пример:
import numpy as np
from scipy.optimize import curve_fit
def func(x, a, b, c):
return a*x**2 + b*x + c
x = np.linspace(0,4,20)
y = func(x, 5, 3, 4)
# generate noisy ydata
yn = y + 0.2 * y * np.random.normal(size=len(x))
# generate error on ydata
y_sigma = 0.2 * y * np.random.normal(size=len(x))
popt, pcov = curve_fit(func, x, yn, sigma = y_sigma)
# plot
import matplotlib.pyplot as plt
fig = plt.figure()
ax = fig.add_subplot(111)
ax.errorbar(x, yn, yerr = y_sigma, fmt = 'o')
ax.plot(x, np.polyval(popt, x), '-')
ax.text(0.5, 100, r"a = {0:.3f} +/- {1:.3f}".format(popt[0], pcov[0,0]**0.5))
ax.text(0.5, 90, r"b = {0:.3f} +/- {1:.3f}".format(popt[1], pcov[1,1]**0.5))
ax.text(0.5, 80, r"c = {0:.3f} +/- {1:.3f}".format(popt[2], pcov[2,2]**0.5))
ax.grid()
plt.show()
Тогда еще кое-что об использовании массивов numpy. Одним из основных преимуществ использования numpy является то, что вы можете избежать циклов for, поскольку операции над массивами применяются поэлементно. Таким образом, цикл for в вашем примере также может быть выполнен следующим образом:
trendx = arange(datasetx[0], (datasetx[-1]+1))
trendy = trend[0]*trendx**2 + trend[1]*trendx + trend[2]
Где я использую arange вместо диапазона, так как он возвращает пустой массив вместо списка.
В этом случае вы также можете использовать функцию numpy polyval:
trendy = polyval(trend, trendx)