python - tutorial - que es numpy

Obtención de errores estándar en parámetros ajustados utilizando el método Optimize.leastsq en python

Question

Tengo un conjunto de datos (desplazamiento frente al tiempo) que he ajustado a un par de ecuaciones usando el método optimise.leastsq. Ahora estoy buscando obtener valores de error en los parámetros ajustados. Mirando a través de la documentación de la matriz de salida es la matriz jacobiana, y debo multiplicar esto por la matriz residual para obtener mis valores. Desafortunadamente, no soy un estadístico, por lo que me estoy ahogando en la terminología.

Por lo que entiendo, todo lo que necesito es la matriz de covarianza que acompaña a mis parámetros ajustados, por lo que puedo cuadrar los elementos diagonales para obtener mi error estándar en los parámetros ajustados. Tengo un vago recuerdo de la lectura de que la matriz de covarianza es lo que se obtiene del método Optimize.leastsq de todos modos. ¿Es esto correcto? Si no, ¿cómo haría para obtener la matriz residual para multiplicar el jacobiano de salida para obtener mi matriz de covarianza?

Cualquier ayuda sería muy apreciada. Soy muy nuevo en Python y, por lo tanto, me disculpo si la pregunta resulta ser básica.

El código de ajuste es el siguiente:

fitfunc = lambda p, t: p[0]+p[1]*np.log(t-p[2])+ p[3]*t # Target function''

errfunc = lambda p, t, y: (fitfunc(p, t) - y)# Distance to the target function

p0 = [ 1,1,1,1] # Initial guess for the parameters


  out = optimize.leastsq(errfunc, p0[:], args=(t, disp,), full_output=1)

Los argumentos t y disp son una matriz de valores de tiempo y descontaminación (básicamente solo 2 columnas de datos). He importado todo lo necesario en la parte superior del código. Los valores ajustados y la matriz proporcionada por la salida son los siguientes:

[  7.53847074e-07   1.84931494e-08   3.25102795e+01  -3.28882437e-11]

[[  3.29326356e-01  -7.43957919e-02   8.02246944e+07   2.64522183e-04]
 [ -7.43957919e-02   1.70872763e-02  -1.76477289e+07  -6.35825520e-05]
 [  8.02246944e+07  -1.76477289e+07   2.51023348e+16   5.87705672e+04]
 [  2.64522183e-04  -6.35825520e-05   5.87705672e+04   2.70249488e-07]]

Sospecho que el ataque es un poco sospechoso de todos modos en este momento. Esto será confirmado cuando pueda sacar los errores.

Answer 1

Actualizado el 4/6/2016

Obtener los errores correctos en los parámetros de ajuste puede ser sutil en la mayoría de los casos.

Pensemos en ajustar una función y=f(x) para la cual tiene un conjunto de puntos de datos (x_i, y_i, yerr_i) , donde i es un índice que se ejecuta sobre cada uno de sus puntos de datos.

En la mayoría de las mediciones físicas, el error yerr_i es una incertidumbre sistemática del dispositivo o procedimiento de medición, por lo que puede considerarse como una constante que no depende de i .

¿Qué función de ajuste usar y cómo obtener los errores de parámetros?

El método Optimize.leastsq devolverá la matriz de covarianza fraccional. Multiplicando todos los elementos de esta matriz por la varianza residual (es decir, el chi cuadrado reducido) y tomando la raíz cuadrada de los elementos diagonales le dará una estimación de la desviación estándar de los parámetros de ajuste. He incluido el código para hacerlo en una de las funciones a continuación.

Por otro lado, si usa optimize.curvefit , la primera parte del procedimiento anterior (multiplicando por el chi cuadrado reducido) se realiza detrás de la escena. Luego, debe tomar la raíz cuadrada de los elementos diagonales de la matriz de covarianza para obtener una estimación de la desviación estándar de los parámetros de ajuste.

Además, optimize.curvefit proporciona parámetros opcionales para tratar casos más generales, donde el valor de yerr_i es diferente para cada punto de datos. De la documentation :

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.

¿Cómo puedo estar seguro de que mis errores son correctos?

Determinar una estimación adecuada del error estándar en los parámetros ajustados es un problema estadístico complicado. Los resultados de la matriz de covarianza, tal como se implementa con optimize.curvefit y optimize.leastsq basan en suposiciones con respecto a la distribución de probabilidad de los errores y las interacciones entre los parámetros; interacciones que pueden existir, dependiendo de su función de ajuste específica f(x) .

En mi opinión, la mejor manera de lidiar con un f(x) complicado es usar el método bootstrap, que se describe en este enlace .

Veamos algunos ejemplos

En primer lugar, algunos código repetitivo. Vamos a definir una función de línea ondulada y generar algunos datos con errores aleatorios. Generaremos un conjunto de datos con un pequeño error aleatorio.

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'')

Ahora, ajustemos la función usando los diversos métodos disponibles:

`optimise.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]

`optimise.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]

Observaciones

Ya comenzamos a ver algo interesante, los parámetros y las estimaciones de error para los tres métodos casi coinciden. ¡Eso es bueno!

Ahora, supongamos que queremos decirle a las funciones de ajuste que existe alguna otra incertidumbre en nuestros datos, tal vez una incertidumbre sistemática que contribuya con un error adicional de veinte veces el valor de err_stdev . Eso es mucho error, de hecho, si simulamos algunos datos con ese tipo de error, se vería así:

Ciertamente, no hay esperanza de que podamos recuperar los parámetros de ajuste con este nivel de ruido.

Para empezar, leastsq que los leastsq no nos permiten ingresar esta nueva información de error sistemático. Veamos qué hace curve_fit cuando le informamos sobre el error:

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]

¿Qué? ¡Esto ciertamente debe estar equivocado!

Esto solía ser el final de la historia, pero recientemente curve_fit agregó el parámetro opcional 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]

Eso es algo mejor, pero aún un poco sospechoso. curve_fit cree que podemos obtener un ajuste de esa señal ruidosa, con un nivel de error del 10% en el parámetro p1 . Veamos lo que tiene que decir 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, tal vez sea una mejor estimación del error en nuestro parámetro de ajuste. bootstrap cree que sabe p1 con un 34% de incertidumbre.

Resumen

optimize.leastsq y optimize.curvefit nos proporcionan una forma de estimar errores en los parámetros ajustados, pero no podemos usar estos métodos sin cuestionarlos un poco. El bootstrap es un método estadístico que usa fuerza bruta, y en mi opinión, tiene una tendencia de trabajar mejor en situaciones que pueden ser más difíciles de interpretar.

Recomiendo encarecidamente ver un problema en particular, y probar curvefit y bootstrap . Si son similares, entonces curvefit es mucho más barato de computar, así que probablemente vale la pena usarlo. Si difieren significativamente, entonces mi dinero estaría en el bootstrap .

Answer 2

Encontré tu pregunta mientras intentaba responder mi question similar. Respuesta corta. El cov_x que se cov_x el mínimo de los resultados debe ser multiplicado por la varianza residual. es decir

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

como en curve_fit.py . Esto se debe a lesssq produce la matriz de covarianza fraccional. Mi gran problema fue que la variación residual aparece como otra cosa cuando se busca en Google.

La variación residual es simplemente una reducción del cuadrado de chi de su ajuste.

Answer 3

Es posible calcular exactamente los errores en el caso de regresión lineal. Y de hecho, la función lesssq da los valores que son diferentes:

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()

salida:

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

Es interesante el resultado que curvefit y bootstrap le darán ...