Почему scipy.optimize панели ошибок для параметров выходят за пределы - PullRequest
0 голосов
/ 11 июля 2019

Я просматривал Получение стандартных ошибок для подогнанных параметров с помощью метода optimize.leastsq в python , поскольку я хочу использовать optimize.least_squares для подгонки некоторых данных, однако у меня есть ограничения (по физическим причинам) на моемпараметры, например, параметр A должен быть [0,1].Я вижу, что метод с использованием ковариационной матрицы и умножения на невязки приводит к появлению полос ошибок, выходящих за пределы параметров.Что нужно сделать статистически правильно и строго?Искусственно отрезать мои панели ошибок на допустимых границах или есть какой-то другой способ позаботиться о том, чтобы панели ошибок были в границах?

import numpy as np
from scipy import optimize
import random

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 #In the fit we will bound p0 from 0 to 1 for e.g., some physical/modeling reason
p1 = 40
p2 = 2.0

# These are initial guesses for fits:
pstart = [0.8,40,2]

# Generate data with a bit of randomness
# (the noise-less function that underlies the data is shown as a blue line)

xdata = np.array(np.linspace(0., 1, 120))
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

def fit_leastsq(p0, datax, datay, function, bounds1):

    errfunc = lambda p, x, y: function(x,p) - y

    all_results= optimize.least_squares(errfunc, p0, args=(datax, datay),bounds=bounds1 )
    pfit=all_results.x
    J=final_results.jac #the jacobian matrix
    pcov=np.linalg.inv(2*J.T.dot(J)) #This calculates the inverse Hessian which is used to approximate the covariance

    s_sq = (errfunc(pfit, datax, datay)**2).sum()/(len(datay)-len(p0))
    pcov = pcov * s_sq

    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 

#Set Bounds on Parameters 
bounds=([0,0,0],[1,40,2.5]) #this means first parameter is restricted 0-1, second is 0-40, and this is 0-2.5

pfit, perr = fit_leastsq(pstart, xdata, ydata, ff, bounds)

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.         39.97612883  1.98430673]
perr =  [0.07545234 0.07611141 0.01546065]

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

1 Ответ

0 голосов
/ 12 июля 2019

Существуют разные способы интерпретации ошибки такого соответствия. Вы пытаетесь думать об этом так: если у меня x_0 с ошибкой x_e, есть 68% -ная вероятность, что реальное значение будет между x_0 - x_e и x_0 + x_e. Алгоритм нахождения вашего x_0 не знает об этом a, поэтому вы можете увидеть ошибку следующим образом: Линейным образом ошибка оценивается по кривизне вашей функции хи-квадрат как минимум. Таким образом, ошибка является лишь мерой «качества соответствия». Более того, он линеаризован. При больших ошибках могут вступать в действие нелинейные свойства ваших функций, которые не фиксируются ни этим анзацем, ни стандартным распространением ошибок.

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

...