Jeppe Gravgaard Jeppe Gravgaard - 27 days ago 9
Python Question

fitting a inverse proportional function

I want to fit the function

f(x) = b+a/x
to my data set. For that I found scipy
from optimize were suitable.

My code is as follows:

x = np.asarray(range(20,401,20))

y is distances that I calculated, but is an array of length 20, here is just random numbers for example

y = np.random.rand(20)

initial guesses of the params a and b

params = np.array([1,1])

function to minimize

def funcinv(x):
return params[0]/x+params[1]

res = least_squares(funinv, params, args=(x, y))

** error given:

return np.atleast_1d(fun(x, *args, **kwargs))

TypeError: funinv() takes 1 positional argument but 3 were given

How can i fit my data?

Answer Source

To make a little of clarity. There are two related problems:

  1. Minimizing a function
  2. Fitting model to data

To fit a model to observed data is to find such parameters of a model which minimize some sort of error between model data and observed data.

least_squares method just minimizes a following function with respect to x (x can be a vector).

F(x) = 0.5 * sum(rho(f_i(x)**2), i = 0, ..., m - 1)

(rho is a loss function and default is rho(x) = x so don't mind it for now)

least_squares(func, x0) expects that call to func(x) will return a vector [a1, a2, a3, ...] for which a sum of squares will be computed: S = 0.5 * (a1^2 + a2^2 + a3^2 + ...).

least_squares will tweak x0 to minimize S.

Thus, in order to use it to fit model to data, one must construct a function of error between a model and actual data - residuals and then minimize that residuals function. In your case you can write it as follows:

import numpy as np
from scipy.optimize import least_squares

x = np.asarray(range(20,401,20))
y = np.random.rand(20)
params = np.array([1,1])

def funcinv(x, a, b):
    return b + a/x

def residuals(params, x, data):
    # evaluates function given vector of params [a, b]
    # and return residuals: (observed_data - model_data)
    a, b = params
    func_eval = funcinv(x, a, b)
    return (data - func_eval)

res = least_squares(residuals, params, args=(x, y))

This gives a result: print(res)

message: '`gtol` termination condition is satisfied.'
    nfev: 4
    njev: 4   optimality: 5.6774618339971994e-10
    status: 1
    success: True
    x: array([ 6.89518618,  0.37118815])

However, as a residuals function pretty much the same all the time (res = observed_data - model_data), there is a shortcut in scipy.optimize called curve_fit: curve_fit(func, xdata, ydata, x0). curve_fit builds residuals function automatically and you can simply write:

import numpy as np
from scipy.optimize import curve_fit

x = np.asarray(range(20,401,20))
y = np.random.rand(20)
params = np.array([1,1])

def funcinv(x, a, b):
    return b + a/x

res = curve_fit(funcinv, x, y, params)
print(res)  #  ... array([ 6.89518618,  0.37118815]), ...