Может соответствовать кривой с Scipy минимизировать, но не с Scipy Curve_fit

Question

Я пытаюсь приспособить функцию y= 1-a(1-bx)**n к некоторым экспериментальным данным, используя scipy curve_fit. Модель существует только для y> 0, поэтому я обрезаю вычисленные значения, чтобы применить это. Код показан ниже

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

# Driver function for scipy.minimize

def driver_func(x, xobs, yobs):

    # Evaluate the fit function with the current parameter estimates

    ynew = myfunc(xobs, *x)
    yerr = np.sum((ynew - yobs) ** 2)

    return yerr

# Define function

def myfunc(x, a, b, n):

    y = 1.0 - a * np.power(1.0 - b * x, n) 
    y = np.clip(y, 0.00, None )

    return y

if __name__ == "__main__":

    # Initialise data

    yobs = np.array([0.005, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.004, 
                    0.048, 0.119, 0.199, 0.277, 0.346, 0.395, 0.444, 0.469, 
                    0.502, 0.527, 0.553, 0.582, 0.595, 0.603, 0.612, 0.599])
    xobs = np.array([0.013, 0.088, 0.159, 0.230, 0.292, 0.362, 0.419, 0.471,
                    0.528, 0.585, 0.639, 0.687, 0.726, 0.772, 0.814, 0.854,
                    0.889, 0.924, 0.958, 0.989, 1.015, 1.045, 1.076, 1.078])

    # Initial guess

    p0 = [2.0, 0.5, 2.0]

    # Check fit pre-regression

    yold = myfunc(xobs, *p0)
    plt.plot(xobs, yobs, 'ko', label='data', fillstyle='none')
    plt.plot(xobs, yold, 'g-', label='pre-fit: a=%4.2f, b=%4.2f, n=%4.2f' % tuple(p0))

    # Fit curve using SCIPY CURVE_FIT

    try:
        popt, pcov = scipy.optimize.curve_fit(myfunc, xobs, yobs, p0=p0)
    except:
        print("Could not fit data using SCIPY curve_fit")
    else:
        ynew = myfunc(xobs, *popt)
        plt.plot(xobs, ynew, 'r-', label='post-curve_fit: a=%4.2f, b=%4.2f, n=%4.2f' % tuple(popt))

    # Fit curve using SCIPY MINIMIZE

    res = scipy.optimize.minimize(driver_func, p0, args=(xobs, yobs), method='Nelder-Mead')
    ynw2 = myfunc(xobs, *res.x)
    plt.plot(xobs, ynw2, 'y-', label='post-minimize: a=%4.2f, b=%4.2f, n=%4.2f' % tuple(res.x))

    plt.legend()
    plt.show()

Я также использовал SCIPY MINIMIZE, чтобы добиться того же. Как показано на рисунке ниже, MINIMIZE работает, но CURVE_FIT в основном исчерпывает оценки и сдается, хотя исходное предположение не так уж далеко от решения MINIMIZE (по крайней мере, визуально). Буду признателен за любые мысли по поводу того, почему функция curve_fit здесь не работает.

Спасибо!

Обновление: после комментариев mikuszefski Я сделал следующие корректировки: 1. удалил вырез из функции подгонки следующим образом:

def myfunc_noclip(x, a, b, n):
    y = 1.0 - a * np.power(1.0 - b * x, n) 
    return y

представил обрезанные массивы, удалив данные ниже порога

ymin = 0.01
xclp = xobs[np.where(yobs >= ymin)]
yclp = yobs[np.where(yobs >= ymin)]

улучшил первоначальное предположение (опять же визуально)

p0 = [1.75, 0.5, 2.0]

обновил вызов кривой_fit

popt, pcov = scipy.optimize.curve_fit(myfunc_noclip, xclp, yclp, p0=p0)

Но это, похоже, не помогло, как показывает следующий график:

Другие посты в stackoverflow, кажется, предполагают, что scipy curve_fit имеет проблемы с подгонкой кривых, где один из параметров подгонки является показателем степени, например, SciPy curve_fit не работает, когда один из подгоняемых параметров является степенью поэтому я предполагаю, что у меня та же проблема. Не уверен, как решить это, хотя ...

Answer 1

Эта проблема вызвана ограничением в определении функции. Два метода минимизации работают принципиально по-разному и, следовательно, очень по-разному реагируют на это ограничение. Здесь minimize используется с Nelder-Mead, который является методом без градиента. Следовательно, алгоритм не рассчитывает числовые градиенты и не оценивает якобиан. В отличие от этого, least-squares, который в конечном итоге называется curve_fit, делает именно это. Однако аппроксимация градиента и, следовательно, любого якобиана несколько сомнительна, если функция не является непрерывной. Как упоминалось ранее, этот разрыв вводится np.clip. После удаления можно легко увидеть, что предположение P0 не так хорошо, как кажется с включенным отсечением. curve_fit сходится с увеличением maxfev=5000, тогда как minimize сразу дает сбой при изменении метода на method='CG'. Чтобы увидеть трудности алгоритмов, можно попытаться вручную указать jac.

. Некоторые примечания: 1) Что касается отсечения, может быть хорошей идеей удалить обрезанные данные, чтобы избежать такой соответствующей проблемы. 2) Глядя на ковариационную матрицу, ошибка n и корреляция с другими значениями чрезвычайно высоки.

Итак, вот что я получаю от

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

# Driver function for scipy.minimize
def driver_func( x, xobs, yobs ):
    # Evaluate the fit function with the current parameter estimates
    ynew = myfunc( xobs, *x)
    yerr = np.sum( ( ynew - yobs ) ** 2 )
    return yerr

# Define functions
def myfunc( x, a, b, n ):
    y = 1.0 - a * np.power( 1.0 - b * x, n ) 
    y = np.clip( y, 0.00, None )
    return y

def myfunc_noclip( x, a, b, n ):
    y = 1.0 - a * np.power( 1.0 - b * x, n ) 
    return y

if __name__ == "__main__":

    # Initialise data
    yobs = np.array([
        0.005, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.004, 
        0.048, 0.119, 0.199, 0.277, 0.346, 0.395, 0.444, 0.469, 
        0.502, 0.527, 0.553, 0.582, 0.595, 0.603, 0.612, 0.599
    ])
    xobs = np.array([
        0.013, 0.088, 0.159, 0.230, 0.292, 0.362, 0.419, 0.471,
        0.528, 0.585, 0.639, 0.687, 0.726, 0.772, 0.814, 0.854,
        0.889, 0.924, 0.958, 0.989, 1.015, 1.045, 1.076, 1.078
    ])

    # Clipped data
    ymin = 0.01
    xclp = xobs[ np.where( yobs >= ymin ) ]
    yclp = yobs[ np.where( yobs >= ymin ) ]

    # Initial guess
    p0 = [ 2.0, 0.5, 2.0 ]

    # Check fit pre-regression
    yold = myfunc( xobs, *p0 )
    plt.plot( xobs, yobs, 'ko', label='data', fillstyle='none' )
    plt.plot( xobs, yold, 'g-', label='pre-fit: a=%4.2f, b=%4.2f, n=%4.2f' % tuple( p0 ) )

    # Fit curve using SCIPY CURVE_FIT
    try:
        popt, pcov = scipy.optimize.curve_fit( myfunc, xobs, yobs, p0=p0, maxfev=5000 )
    except:
        print("Could not fit data using SCIPY curve_fit")
    else:
        ynew = myfunc( xobs, *popt )
        plt.plot( xobs, ynew, 'r-', label="curve-fit: a=%4.2f, b=%4.2e, n=%4.2f" % tuple( popt ) )

    # Fit curve using SCIPY CURVE_FIT on clipped data
    p0 = [ 1.75, 1e-4, 1e3 ]
    try:
        popt, pcov = scipy.optimize.curve_fit( myfunc_noclip, xclp, yclp, p0=p0 )
    except:
        print("Could not fit data using SCIPY curve_fit")
    else:
        ynew = myfunc_noclip( xobs, *popt )
        plt.plot( xobs, ynew, 'k-', label="curve-fit clipped data: a=%4.2f, b=%4.2e, n=%4.2f" % tuple( popt ) )

    # Fit curve using SCIPY MINIMIZE
    p0 = [ 2.0, 0.5, 2.0 ]
    res = scipy.optimize.minimize( driver_func, p0, args=( xobs, yobs ), method='Nelder-Mead' )
    # ~res = scipy.optimize.minimize(driver_func, p0, args=(xobs, yobs), method='CG')
    ynw2 = myfunc( xobs, *res.x )
    plt.plot( xobs, ynw2, 'y--', label='Nelder-Mead 1: a=%4.2f, b=%4.2f, n=%4.2f' % tuple( res.x ) )
    p0 = [ 2.4, 3.6e-4, 5.6e3 ]
    res = scipy.optimize.minimize( driver_func, p0, args=( xobs, yobs ), method='Nelder-Mead' )
    # ~res = scipy.optimize.minimize(driver_func, p0, args=(xobs, yobs), method='CG')
    ynw2 = myfunc( xobs, *res.x )
    plt.plot( xobs, ynw2, 'b:', label='Nelder-Mead 2: a=%4.2f, b=%4.2e, n=%4.2e' % tuple( res.x ) )

    plt.legend( loc=2 )
    plt.ylim( -0.05, 0.7 )
    plt.grid()
    plt.show()

Так что я бы сказал, что это работает okeyi sh. Однако я получаю предупреждение о переполнении.