Ответ 1
Обновлено 4/6/2016
Получение правильных ошибок в параметрах подгонки может быть тонким в большинстве случаев.
Подумайте об установке функции y=f(x), для которой у вас есть набор точек данных (x_i, y_i, yerr_i), где i - это индекс, который проходит по каждой из ваших точек данных.
В большинстве физических измерений ошибка yerr_i представляет собой систематическую неопределенность измерительного устройства или процедуры, и поэтому ее можно рассматривать как константу, которая не зависит от i.
Какую функцию фитинга использовать и как получить ошибки параметра?
Метод optimize.leastsq вернет дробную ковариационную матрицу. Умножая все элементы этой матрицы на остаточную дисперсию (т.е. Приведенную схему хи) и принимая квадратный корень из диагональных элементов, даст вам оценку стандартного отклонения параметров подгонки. Я включил код для этого в одну из следующих функций.
С другой стороны, если вы используете optimize.curvefit, первая часть вышеописанной процедуры (умножение на уменьшенную квадратуру хи) выполняется для вас за кулисами. Затем вам нужно взять квадратный корень из диагональных элементов матрицы ковариации, чтобы получить оценку стандартного отклонения параметров подгонки.
Кроме того, optimize.curvefit предоставляет необязательные параметры для работы с более общими случаями, где значение yerr_i отличается для каждой точки данных. Из документация:
sigma : None or M-length sequence, optional
If not None, the uncertainties in the ydata array. These are used as
weights in the least-squares problem
i.e. minimising ``np.sum( ((f(xdata, *popt) - ydata) / sigma)**2 )``
If None, the uncertainties are assumed to be 1.
absolute_sigma : bool, optional
If False, `sigma` denotes relative weights of the data points.
The returned covariance matrix `pcov` is based on *estimated*
errors in the data, and is not affected by the overall
magnitude of the values in `sigma`. Only the relative
magnitudes of the `sigma` values matter.
Как я могу быть уверен, что мои ошибки верны?
Определение правильной оценки стандартной ошибки в установленных параметрах является сложной статистической задачей. Результаты ковариационной матрицы, реализованные optimize.curvefit и optimize.leastsq, фактически опираются на предположения относительно распределения вероятностей ошибок и взаимодействий между параметрами; взаимодействия, которые могут существовать, в зависимости от вашей конкретной функции соответствия f(x).
На мой взгляд, лучший способ справиться со сложным f(x) заключается в использовании метода начальной загрузки, который описан в этой ссылке.
Посмотрим несколько примеров
Во-первых, некоторый шаблонный код. Пусть задает функцию коротковолновой линии и генерирует некоторые данные со случайными ошибками. Мы создадим набор данных с небольшой случайной ошибкой.
import numpy as np
from scipy import optimize
def f( x, p0, p1, p2):
return p0*x + 0.4*np.sin(p1*x) + p2
def ff(x, p):
return f(x, *p)
# These are the true parameters
p0 = 1.0
p1 = 40
p2 = 2.0
# These are initial guesses for fits:
pstart = [
p0 + random.random(),
p1 + 5.*random.random(),
p2 + random.random()
]
%matplotlib inline
import matplotlib.pyplot as plt
xvals = np.linspace(0., 1, 120)
yvals = f(xvals, p0, p1, p2)
# Generate data with a bit of randomness
# (the noise-less function that underlies the data is shown as a blue line)
xdata = np.array(xvals)
np.random.seed(42)
err_stdev = 0.2
yvals_err = np.random.normal(0., err_stdev, len(xdata))
ydata = f(xdata, p0, p1, p2) + yvals_err
plt.plot(xvals, yvals)
plt.plot(xdata, ydata, 'o', mfc='None')
Теперь установите функцию, используя различные доступные методы:
`optimize.leastsq`
def fit_leastsq(p0, datax, datay, function):
errfunc = lambda p, x, y: function(x,p) - y
pfit, pcov, infodict, errmsg, success = \
optimize.leastsq(errfunc, p0, args=(datax, datay), \
full_output=1, epsfcn=0.0001)
if (len(datay) > len(p0)) and pcov is not None:
s_sq = (errfunc(pfit, datax, datay)**2).sum()/(len(datay)-len(p0))
pcov = pcov * s_sq
else:
pcov = np.inf
error = []
for i in range(len(pfit)):
try:
error.append(np.absolute(pcov[i][i])**0.5)
except:
error.append( 0.00 )
pfit_leastsq = pfit
perr_leastsq = np.array(error)
return pfit_leastsq, perr_leastsq
pfit, perr = fit_leastsq(pstart, xdata, ydata, ff)
print("\n# Fit parameters and parameter errors from lestsq method :")
print("pfit = ", pfit)
print("perr = ", perr)
# Fit parameters and parameter errors from lestsq method :
pfit = [ 1.04951642 39.98832634 1.95947613]
perr = [ 0.0584024 0.10597135 0.03376631]
`optimize.curve_fit`
def fit_curvefit(p0, datax, datay, function, yerr=err_stdev, **kwargs):
pfit, pcov = \
optimize.curve_fit(f,datax,datay,p0=p0,\
sigma=yerr, epsfcn=0.0001, **kwargs)
error = []
for i in range(len(pfit)):
try:
error.append(np.absolute(pcov[i][i])**0.5)
except:
error.append( 0.00 )
pfit_curvefit = pfit
perr_curvefit = np.array(error)
return pfit_curvefit, perr_curvefit
pfit, perr = fit_curvefit(pstart, xdata, ydata, ff)
print("\n# Fit parameters and parameter errors from curve_fit method :")
print("pfit = ", pfit)
print("perr = ", perr)
# Fit parameters and parameter errors from curve_fit method :
pfit = [ 1.04951642 39.98832634 1.95947613]
perr = [ 0.0584024 0.10597135 0.03376631]
`bootstrap`
def fit_bootstrap(p0, datax, datay, function, yerr_systematic=0.0):
errfunc = lambda p, x, y: function(x,p) - y
# Fit first time
pfit, perr = optimize.leastsq(errfunc, p0, args=(datax, datay), full_output=0)
# Get the stdev of the residuals
residuals = errfunc(pfit, datax, datay)
sigma_res = np.std(residuals)
sigma_err_total = np.sqrt(sigma_res**2 + yerr_systematic**2)
# 100 random data sets are generated and fitted
ps = []
for i in range(100):
randomDelta = np.random.normal(0., sigma_err_total, len(datay))
randomdataY = datay + randomDelta
randomfit, randomcov = \
optimize.leastsq(errfunc, p0, args=(datax, randomdataY),\
full_output=0)
ps.append(randomfit)
ps = np.array(ps)
mean_pfit = np.mean(ps,0)
# You can choose the confidence interval that you want for your
# parameter estimates:
Nsigma = 1. # 1sigma gets approximately the same as methods above
# 1sigma corresponds to 68.3% confidence interval
# 2sigma corresponds to 95.44% confidence interval
err_pfit = Nsigma * np.std(ps,0)
pfit_bootstrap = mean_pfit
perr_bootstrap = err_pfit
return pfit_bootstrap, perr_bootstrap
pfit, perr = fit_bootstrap(pstart, xdata, ydata, ff)
print("\n# Fit parameters and parameter errors from bootstrap method :")
print("pfit = ", pfit)
print("perr = ", perr)
# Fit parameters and parameter errors from bootstrap method :
pfit = [ 1.05058465 39.96530055 1.96074046]
perr = [ 0.06462981 0.1118803 0.03544364]
Наблюдения
Мы уже начинаем видеть что-то интересное, параметры и оценки ошибок для всех трех методов почти согласуются. Это хорошо!
Теперь предположим, что мы хотим сказать функциям подгонки, что в наших данных есть какая-то другая неопределенность, возможно, систематическая неопределенность, которая вносит дополнительную ошибку в двадцать раз больше значения err_stdev. На самом деле это большая ошибка, если мы моделируем некоторые данные с такой ошибкой, это будет выглядеть так:
Конечно, нет надежды на то, что мы сможем восстановить подходящие параметры с таким уровнем шума.
Для начала давайте осознаем, что leastsq даже не позволяет нам вводить эту новую систематическую информацию об ошибке. Посмотрим, что делает curve_fit, когда мы рассказываем об ошибке:
pfit, perr = fit_curvefit(pstart, xdata, ydata, ff, yerr=20*err_stdev)
print("\nFit parameters and parameter errors from curve_fit method (20x error) :")
print("pfit = ", pfit)
print("perr = ", perr)
Fit parameters and parameter errors from curve_fit method (20x error) :
pfit = [ 1.04951642 39.98832633 1.95947613]
perr = [ 0.0584024 0.10597135 0.03376631]
Whaat?? Это, безусловно, должно быть неправильно!
Это был конец истории, но в последнее время curve_fit добавил необязательный параметр absolute_sigma:
pfit, perr = fit_curvefit(pstart, xdata, ydata, ff, yerr=20*err_stdev, absolute_sigma=True)
print("\n# Fit parameters and parameter errors from curve_fit method (20x error, absolute_sigma) :")
print("pfit = ", pfit)
print("perr = ", perr)
# Fit parameters and parameter errors from curve_fit method (20x error, absolute_sigma) :
pfit = [ 1.04951642 39.98832633 1.95947613]
perr = [ 1.25570187 2.27847504 0.72600466]
Это несколько лучше, но все же немного подозрительно. curve_fit думает, что мы можем получить соответствие из этого шумного сигнала с уровнем ошибки 10% в параметре p1. Посмотрим, что должен сказать bootstrap:
pfit, perr = fit_bootstrap(pstart, xdata, ydata, ff, yerr_systematic=20.0)
print("\nFit parameters and parameter errors from bootstrap method (20x error):")
print("pfit = ", pfit)
print("perr = ", perr)
Fit parameters and parameter errors from bootstrap method (20x error):
pfit = [ 2.54029171e-02 3.84313695e+01 2.55729825e+00]
perr = [ 6.41602813 13.22283345 3.6629705 ]
Ah, возможно, это лучшая оценка ошибки в нашем параметре fit. bootstrap думает, что знает p1 примерно с 34% неопределенностью.
Резюме
optimize.leastsq и optimize.curvefit предоставляют нам способ оценки ошибок в установленных параметрах, но мы не можем просто использовать эти методы, не ставя под сомнение их немного. bootstrap - статистический метод, который использует грубую силу и, на мой взгляд, имеет тенденцию работать лучше в ситуациях, которые могут быть труднее интерпретировать.
Я настоятельно рекомендую посмотреть конкретную проблему и попробовать curvefit и bootstrap. Если они похожи, то curvefit намного дешевле вычислять, поэтому, вероятно, стоит использовать. Если они значительно отличаются, тогда мои деньги будут на bootstrap.