Получение стандартных ошибок для установленных параметров с использованием метода optimize.leastsq в python

Ответ 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.

Ответ 2

Нашел ваш вопрос, пытаясь ответить на мой собственный question. Короткий ответ. cov_x, что наименьшие выходы должны быть умножены на остаточную дисперсию. то есть.

s_sq = (func(popt, args)**2).sum()/(len(ydata)-len(p0))
pcov = pcov * s_sq

как в curve_fit.py. Это связано с тем, что lesssq выводит дробную ковариационную матрицу. Моя большая проблема заключалась в том, что остаточная дисперсия проявляется как что-то еще при поиске в ней.

Остаточная дисперсия просто уменьшена на квадрат от вашей подгонки.

Ответ 3

Можно точно вычислить ошибки в случае линейной регрессии. И действительно, функция lesssq дает разные значения:

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

A = 1.353
B = 2.145

yerr = 0.25

xs = np.linspace( 2, 8, 1448 )
ys = A * xs + B + np.random.normal( 0, yerr, len( xs ) )

def linearRegression( _xs, _ys ):
    if _xs.shape[0] != _ys.shape[0]:
        raise ValueError( 'xs and ys must be of the same length' )
    xSum = ySum = xxSum = yySum = xySum = 0.0
    numPts = 0

    for i in range( len( _xs ) ):
        xSum += _xs[i]
        ySum += _ys[i]
        xxSum += _xs[i] * _xs[i]
        yySum += _ys[i] * _ys[i]
        xySum += _xs[i] * _ys[i]
        numPts += 1

    k = ( xySum - xSum * ySum / numPts ) / ( xxSum - xSum * xSum / numPts )
    b = ( ySum - k * xSum ) / numPts
    sk = np.sqrt( ( ( yySum - ySum * ySum / numPts ) / ( xxSum - xSum * xSum / numPts ) - k**2 ) / numPts )
    sb = np.sqrt( ( yySum - ySum * ySum / numPts ) - k**2 * ( xxSum - xSum * xSum / numPts ) ) / numPts

    return [ k, b, sk, sb ]

def linearRegressionSP( _xs, _ys ):
    defPars = [ 0, 0 ]
    pars, pcov, infodict, errmsg, success = \
    leastsq( lambda _pars, x, y: y - ( _pars[0] * x + _pars[1] ), defPars, args = ( _xs, _ys ), full_output=1 )

    errs = []
    if pcov is not None:
        if( len( _xs ) > len(defPars) ) and pcov is not None:
            s_sq = ( ( ys - ( pars[0] * _xs + pars[1] ) )**2 ).sum() / ( len( _xs ) - len( defPars ) )
            pcov *= s_sq

        for j in range( len( pcov ) ):
            errs.append( pcov[j][j]**0.5 )
    else:
        errs = [ np.inf, np.inf ]
    return np.append( pars, np.array( errs ) )

regr1 = linearRegression( xs, ys )
regr2 = linearRegressionSP( xs, ys )

print( regr1 )
print( 'Calculated sigma = %f' %  ( regr1[3] * np.sqrt( xs.shape[0] ) ) )
print( regr2 )
#print( 'B = %f must be in ( %f,\t%f )' % ( B, regr1[1] - regr1[3], regr1[1] + regr1[3] ) )


plt.plot( xs, ys, 'bo' )
plt.plot( xs, regr1[0] * xs + regr1[1] )
plt.show()

выход:

[1.3481681543925064, 2.1729338701374137, 0.0036028493647274687, 0.0062446292528624348]
Calculated sigma = 0.237624 # quite close to yerr
[ 1.34816815  2.17293387  0.00360534  0.01907908]

Интересно, какие результаты curvefit и bootstrap дадут...