Search code examples
pythonscipycurve-fittingbounds

How to set bounds for only one parameter


I'm using curve_fit from scipy.optimize to fit my data. I have a function that fits three parameters (Z1, Z2, Z3). I wannt to provide bounds. However, I'd like to only provide a bound to Z2 (Z2 shall be below 40). I do not want to set bounds for Z1 and Z3. Is that possible?

            popt, pcov = curve_fit(func, xdata, ydata, p0 = [Z1, Z2, Z3],
                           bounds = ((10, 20, 5), (100, 50, 100,)))

            # This way I provide bounds to Z1, Z2 and Z3
            # I, however, only want to say that Z2 < 40
            # Also interesting would be to say Z2 < Z1, with no bounds for Z1 or Z3

Solution

  • Here I provide a pseudo-code solution that uses just remapping of parameters instead of actual boundary conditions

    originally we would do something like:

    bestA, bestB, bestC = fit( function( x, a, b, c, guesses=( guessA, guessB, guessC) ) )
    

    However, we want the restriction that b < c. In such a case we fit the following instead (omitting the guesses)

    wrapper_function( x, a, d, c ) = function( x, a, c - d**2, c )   
    bestA, bestD, bestC = fit( wrapper_function( x, a, d, c ) )
    

    like this the value b send to function() will always be less than c. Assuming the fit converges we just calculate b = c - d**2. If we are also interested in the error of b we have to do error-propagation including correlation. ( see e.g. Joel Tellinghuisen, J. Phys. Chem. A 2001, 105, 3917-3921)

    So s_b**2 = gT V g. where V is the variance-covariance matrix and g = df/du and u in [a, d, c]. That is gT=(0, -2 * bestD, 1 ).

    Now we want b < min( c, 40 ). That is, as mentioned in my comment, a bit more complicated, but also possible. We rewrite the wrapper function and have

    wrapper_function( x, a, d, c ) = function( x, a, 0.5 * ( c + 40 - abs( c - 40 ) ) - d**2, c )   
    bestA, bestD, bestC = fit( wrapper_function( x, a, d, c ) )
    

    It might not be obvious that this does the trick, but if one plots 0.5 * ( c + 40 - abs( c - 40 ) ) it becomes clear. It is again straight forward to calculate b. For the error we have to be careful with calculating g, though. We get

    g = ( 0, -2 * bestD, 1 - numpy.heaviside( bestC - 40, 0.5 ) )
    

    Note: if the value and error of c are such that the discontinuity of the remapping is within error margins, this needs to be reconsidered, though.