Tengo un conjunto de datos (desplazamiento vs tiempo) que he instalado un par de ecuaciones mediante el optimizar.leastsq método. Ahora estoy mirando para obtener los valores de error en el conjunto de los parámetros. Mirando a través de la documentación de la matriz mostrada es la matriz jacobiana, y debo multiplicarlo por el valor residual de la matriz para obtener mis valores. Por desgracia yo no soy un estadístico así que me estoy ahogando un poco en la terminología.

Por lo que puedo entender todo lo que necesito es la matriz de covarianza que va con mi equipada parámetros, por lo que puedo raíz cuadrada de los elementos de la diagonal para obtener mi error estándar en el conjunto de los parámetros. Tengo un vago recuerdo de la lectura que la matriz de covarianza es ¿cuál es la salida de la a optimizar.leastsq método de todos modos. Es esto correcto? Si no ¿cómo ir sobre cómo obtener el valor residual de la matriz de multiplicar la emitida Jacobiana por conseguir mi matriz de covarianza?

Cualquier ayuda sería muy apreciada. Soy muy nuevo en python y por lo tanto disculpas si la pregunta resulta ser uno básico.

el correspondiente código 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)

La args t y disp es y la matriz de tiempo y displcement valores (básicamente, sólo 2 columnas de datos). Me ha importado todo lo necesario en el tope del código. Los valores ajustados y la matriz proporcionada por el resultado es como sigue:

[  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 ajuste es un poco sospechoso de todos modos por el momento. Esto será confirmado cuando puedo conseguir los errores a cabo.

  • docs.scipy.org/doc/scipy/reference/generated/…
  • No del todo seguro de qué preguntar, tal vez algunos ejemplos podrían ser de utilidad.
  • Esa es la principal docs he estado mirando. El código que he utilizado ha sido añadido a la descripción anterior
  • Creo que he encontrado una manera alrededor de ella (aunque sea un poco incómodo, en términos de reescritura de código), yo lo de la ‘optimizar.curve_fit’ salidas de la covarience de la matriz, de la que usted puede conseguir sus errores, y utiliza la misma mínimos cuadrados el método de regresión como el ‘optimizar.leastsq’. ¿Alguien puede confirmar que esto es correcto?
  • Sí, curve_fit devuelve la matriz de covarianza para la estimación del parámetro (la incertidumbre). Si desea utilizar leastsq directamente, también se puede comprobar la fuente de curve_fit.
InformationsquelleAutor Phil | 2013-01-29

3 Comentarios

  1. 74

    Actualizado en 4/6/2016

    Conseguir la corrección de errores en el ajuste de los parámetros pueden ser sutiles en la mayoría de los casos.

    Vamos a pensar en el ajuste de una función y=f(x) para los que 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 los puntos de datos.

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

    Que ajuste la función a utilizar, y cómo obtener el parámetro de errores?

    La optimize.leastsq método devuelve la fracción de una matriz de covarianza. Multiplicando todos los elementos de esta matriz por la varianza residual (es decir, la reducción de la chi al cuadrado) y tomando la raíz cuadrada de los elementos de la diagonal le dará una estimación de la desviación estándar del ajuste de los parámetros. He incluido el código para hacer que en una de las siguientes funciones.

    Por otro lado, si utiliza optimize.curvefit, la primera parte del procedimiento anterior (multiplicando por la reducción de la chi al cuadrado) está hecho para usted detrás de las escenas. A continuación, deberá tomar la raíz cuadrada de los elementos de la diagonal de la matriz de covarianza para obtener una estimación de la desviación estándar del ajuste de parámetros.

    Además, optimize.curvefit proporciona parámetros opcionales para lidiar con los casos más generales, donde el yerr_i valor es diferente para cada punto de datos. Desde el documentación:

    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 correctas?

    La determinación de una adecuada estimación del error estándar en el conjunto de los parámetros es un complicado problema estadístico. Los resultados de la matriz de covarianza, como el implementado por optimize.curvefit y optimize.leastsq en realidad se basan en supuestos acerca de la distribución de probabilidad de los errores y de las interacciones entre los parámetros; las interacciones que pueden existir, dependiendo del ajuste de la función de f(x).

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

    Vamos a ver algunos ejemplos

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

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

    Obtener los errores estándar en armarios parámetros utilizando el optimizar.leastsq método en python

    Ahora, vamos a ajuste de la función utilizando los diferentes métodos disponibles:

    optimizar.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]

    optimizar.curve_fit`

    def fit_curvefit(p0, datax, datay, function, yerr=err_stdev, **kwargs):
        """
        Note: As per the current documentation (Scipy V1.1.0), sigma (yerr) must be:
            None or M-length sequence or MxM array, optional
        Therefore, replace:
            err_stdev = 0.2
        With:
            err_stdev = [0.2 for item in xdata]
        Or similar, to create an M-length sequence for this example.
        """
        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 hemos empezado a ver algo interesante, los parámetros y estimaciones de error para los tres métodos casi de acuerdo. Eso es bueno!

    Ahora, supongamos que queremos decir a las funciones de ajuste que hay algunos otros incertidumbre en nuestros datos, tal vez una sistemática de incertidumbre que podrían contribuir a un error adicional de veinte veces el valor de err_stdev. Que es mucho de error, de hecho, si queremos simular algunos de los datos con que tipo de error se vería así:

    Obtener los errores estándar en armarios parámetros utilizando el optimizar.leastsq método en python

    Ciertamente no hay esperanza de que podemos recuperar el ajuste de los parámetros de este nivel de ruido.

    Para empezar, vamos a darnos cuenta de que leastsq no nos permiten la entrada de este nuevo error sistemático de la información. Vamos a ver lo que curve_fit ¿ cuando nos dicen que 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]

    Whaat?? Esto sin duda debe ser un error!

    Este solía ser el final de la historia, pero recientemente curve_fit añadido la absolute_sigma parámetro opcional:

    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]

    Que es algo mejor, pero todavía un poco de pescado. curve_fit piensa que puede obtener un ajuste de esa señal con ruido, con un nivel de error del 10% en el p1 parámetro. Vamos a ver lo que bootstrap tiene que decir:

    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, que es tal vez una mejor estimación del error en el ajuste de parámetros. bootstrap piensa que sabe p1 con cerca de un 34% en la incertidumbre.

    Resumen

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

    Yo recomiendo mirar un problema en particular, y tratando de curvefit y bootstrap. Si son similares, entonces curvefit es mucho más barato para calcular, por lo que probablemente vale la pena utilizar. Si difieren significativamente, luego de mi dinero estaría en la bootstrap.

    • Así que, déjame ver si entiendo lo que usted hizo. Generar valores aleatorios de todo el suministrado Y basado en el adaptador de errores (o residuo) de hacerlo 100 veces, obtener el ajuste de los parámetros en cada uno de los pasos. A continuación, tomar la SEXUAL, como el error? ¿Cómo es en comparación con statsmodels’ JUEGO?
    • Hola, Yotam, no estoy familiarizado con el PREMIO de la función en statsmodels. Lo que yo busco es tratar con las mediciones que tienen individualmente un gran error/uncertainy pero cuyos medios se ajustan a un modelo razonablemente bien. El objetivo es proporcionar un valor razonable de que el error del ajuste de los parámetros. Incluso si la media de los valores se encuentran perfectamente en la línea de ajuste, si el error de cada medida es gigantesco, a continuación, el ajuste del parámetro de error debe reflejar eso.
    • Una forma de hacerlo es la ‘fuerza bruta’ método de arranque (que he mostrado aquí), de otra manera sería estudiar la dependencia de la función de error (chi cuadrado, $\chi^{2}$) en cada uno equipado parámetro y el establecimiento de un intervalo de confianza para cada parámetro. Algunas personas utilizan el conocido derivada de la función de error con respecto a la equipada con parámetro ( $\mathrm{d}\chi^{2} / \mathrm{d}p$ ) para acceder rápidamente a este intervalo de confianza, o que numéricamente variar el parámetro con el resto de los parámetros mantiene fijo y definir un intervalo de confianza basado en el cambio observado de $\chi^{2}$
    • Hola, tu respuesta es muy instructivo y útil. Usted ha hablado de «Un beneficio adicional de optmize.curvefit es que no se restringe para el caso de la IGUALDAD de los ERRORES», ahora tengo un ajuste de la función con la DESIGUALDAD de los Errores, pero no sé cómo implementar el uso de Scipy curve_fit. Me podrían ayudar con eso? Muchas gracias. Mi pregunta es en el siguiente link: stackoverflow.com/questions/36417506/…
    • usted tiene que usar el sigma parámetro opcional en la curva de ajuste, y de paso es una matriz que tiene los errores para cada una de sus datos y puntos.
    • también he actualizado esta respuesta tal vez para hacerla más clara. Doy la bienvenida a cualquier comentario que pueda tener, que puede ayudar a mejorar la respuesta.
    • Muchas gracias por su instructiva explicación.
    • la aplicación de su código de inicio a mi curve_fit devuelve TypeError: func() takes exactly 4 arguments (2 given) que apunta en la línea donde se errfunc = lambda p, x, y: function(x,p) - y. Cómo resolverlo?
    • Tras llamar a bootstrap, he dado los argumentos de derecho. Mi ajuste de la función es a * numpy.exp(-c*x)+d. Es este el problema?
    • el p0 argumento en fit_bootstrap es la matriz-como. Por lo que su función debe ser definida como: def fit_function(x, p): p[0] * numpy.exp(-p[1]*x) + p[2]. Déjeme saber si funciona para usted.
    • Parámetros indicados para el ejemplo se puede extraer con precisión por transformada de fourier (frecuencia y parámetros dc) así que no es tan «sospechoso» que los métodos de «precisión» determinar los parámetros. Intente rehacer con un simple polinomio o una oscilación de la onda sinusoidal y ver cómo pescado es entonces.
    • lo que me parece sospechoso no es el hecho de que los parámetros pueden ser obtenidos, pero que los errores de parámetros que se obtienen son el 10% del valor. Como los datos es más ruidoso que me lo esperaba más grande de los errores en los valores ajustados de los parámetros.
    • Mi punto es que esa suposición es malo para el problema de ejemplo. Esto está relacionado con el bloqueo en la detección en la que una de oscilación de la señal puede ser sacado de una terrible ruido de fondo (un LED intermitente puede ser detectado en un fotodiodo metros inundado por la luz de la sala con todos los 50-60 Hertz ruido, etc). Los algoritmos que dan resultados diferentes debido a inadecuado escalas (como usted bien señala), pero el algoritmo bootstrap en este caso es una mala elección, porque se debería haber hecho un trabajo mucho mejor – me hincapié «en este caso».
    • Yo diría que la comparación de este bloqueo en la detección no es una comparación justa. Para hacer algo como bloqueo en la detección usted necesita saber la frecuencia de la fuente de antemano. En este ejemplo, la frecuencia de la onda sinusoidal es uno de los parámetros que estamos estimando.
    • Excelente explicación! Sólo una pequeña confirmación. Así que si tengo errores individuales para cada uno de mis y valores, luego de obtener el corregir los errores en el ajuste de los parámetros, necesito usar absolute_sigma = True, estoy en lo correcto? Si no la uso, puedo ver que los errores obtenidos son muy pequeñas, ¿por qué es esto así? ¿Qué absolute_sigma = False (el caso predeterminado) ¿cuando cada y valor tiene su propio error ?
    • Por defecto, los valores de sigma debe ser entendida sólo como pesos a los puntos de datos individuales. Los pesos serán utilizados por curvefit y tendrá un efecto sobre los valores ajustados pero no en la matriz de covarianza. Cuando dices absolute_sigma=True que están diciendo curvefit que su sigmas no son sólo acaba de pesos, pero también real incertidumbres en las mismas unidades que $y$. En ese caso curvefit va por delante y devuelve la onu-normalizado de la matriz de covarianza en pcov. De esa manera usted puede tomar la raíz cuadrada de la diagonal entradas en pcov obtener directamente el parámetro de errores.
    • para optimize.curve_fit puedo obtener ValueError: sigma` ha incorrecto forma`. Otros se encontraron con este problema?
    • Me acaba de editar en un docstring a explicar cómo solucionar ese problema en el ejemplo 🙂
    • sí, funciona ahora, de hecho,

  2. 10

    Encontrado a su pregunta, mientras que tratando de responder a mis propias similar pregunta.
    Respuesta corta. El cov_x que leastsq salidas 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 es debido a que leastsq salidas de las fracciones de una matriz de covarianza. Mi gran problema fue que la varianza residual se muestra como algo más que buscar en google es.

    Varianza Residual se reduce simplemente chi cuadrado de su ajuste.

  3. 2

    Es posible calcular exactamente los errores en el caso de la regresión lineal. Y, de hecho, la leastsq función 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()

    de 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 que los resultados de curvefit y bootstrap dará…

Dejar respuesta

Please enter your comment!
Please enter your name here