I'm trying to fit a function using curve_fit from the scipy.optimize module. I keep getting the error message:
"TypeError: only size-1 arrays can be converted to Python scalars'.
frequency = []
for i in driven_period:
frequency.append(2*math.pi*i**-1)
amplitude = []
for i in driven_height:
amplitude.append(i-41.3)
def amplitude_func(omega, p, omega0, gamma):
x=((p*omega0**2)/(math.sqrt((omega0**2-omega**2)**2+gamma**2*omega**2)))
return x
parameters = sp.curve_fit(amplitude_func, frequency, amplitude)
p_fit = parameters[0][0]
omega0_fit = parameters[0][1]
gamma_fit = parameters[0][2]
fit=[]
for i in frequency:
fit.append(amplitude_func(i, p_fit, w_fit, y_fit))
plt.scatter(frequency, amplitude)
plt.plot(frequency, fit)
I've tried this same code with a basic linear function and it worked, so I think the problem must be with the function amplitude_func, but I can't figure out what it is. I've tried printing at every stage, but I cant find where I have the mistake. Sorry if this is a silly question, I'm new to programming. Help appreciated.
Edit: more info. Error message:
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
Input In [151], in <cell line: 27>()
24 return (p*(omega)**2)/(math.sqrt((omega0)**2-(omega)**2)**2+(gamma)**2*(omega**2))
26 # lorenzian least-squares fit
---> 27 amplitude_fit, amplitude_covariance = sp.curve_fit(amplitude_func, frequency, amplitude)
28 p_fit = amplitude_fit[0]
29 omega0_fit = amplitude_fit[1]
File ~/opt/anaconda3/lib/python3.9/site-packages/scipy/optimize/minpack.py:789, in curve_fit(f, xdata, ydata, p0, sigma, absolute_sigma, check_finite, bounds, method, jac, **kwargs)
787 # Remove full_output from kwargs, otherwise we're passing it in twice.
788 return_full = kwargs.pop('full_output', False)
--> 789 res = leastsq(func, p0, Dfun=jac, full_output=1, **kwargs)
790 popt, pcov, infodict, errmsg, ier = res
791 ysize = len(infodict['fvec'])
File ~/opt/anaconda3/lib/python3.9/site-packages/scipy/optimize/minpack.py:410, in leastsq(func, x0, args, Dfun, full_output, col_deriv, ftol, xtol, gtol, maxfev, epsfcn, factor, diag)
408 if not isinstance(args, tuple):
409 args = (args,)
--> 410 shape, dtype = _check_func('leastsq', 'func', func, x0, args, n)
411 m = shape[0]
413 if n > m:
File ~/opt/anaconda3/lib/python3.9/site-packages/scipy/optimize/minpack.py:24, in _check_func(checker, argname, thefunc, x0, args, numinputs, output_shape)
22 def _check_func(checker, argname, thefunc, x0, args, numinputs,
23 output_shape=None):
---> 24 res = atleast_1d(thefunc(*((x0[:numinputs],) + args)))
25 if (output_shape is not None) and (shape(res) != output_shape):
26 if (output_shape[0] != 1):
File ~/opt/anaconda3/lib/python3.9/site-packages/scipy/optimize/minpack.py:485, in _wrap_func.<locals>.func_wrapped(params)
484 def func_wrapped(params):
--> 485 return func(xdata, *params) - ydata
Input In [151], in amplitude_func(omega, p, omega0, gamma)
23 def amplitude_func(omega, p, omega0, gamma):
---> 24 return (p*(omega)**2)/(math.sqrt((omega0)**2-(omega)**2)**2+(gamma)**2*(omega**2))
TypeError: only size-1 arrays can be converted to Python scalars
Arrays:
driven_period = [1.027, 0.949, 0.866, 0.806, 0.648, 0.917, 0.935, 1.404, 1.113, 0.887, 0.983]
driven_height = [43.1, 45.35, 43.7, 42.25, 41.3, 45.9, 45.7, 41.5, 42.25, 44.75, 44.0]
I think I figuered it out:
Python seems to have trouble multiplying arrays ( ** operator as well as math library). Switching over to np for calculations works. Sadly I can’t explain in depth why :/
import math
import numpy as np
driven_period = [1.027, 0.949, 0.866, 0.806, 0.648, 0.917, 0.935, 1.404, 1.113, 0.887, 0.983]
driven_height = [43.1, 45.35, 43.7, 42.25, 41.3, 45.9, 45.7, 41.5, 42.25, 44.75, 44.0]
frequency = []
for i in driven_period:
frequency.append(2*math.pi*i**-1)
amplitude = []
for i in driven_height:
amplitude.append(i-41.3)
def amplitude_func(omega, p, omega0, gamma):
#x=((p*omega0**2)/(math.sqrt((omega0**2-omega**2)**2+gamma**2*omega**2)))
p1 = np.power(omega0,2)-np.power(omega,2)
p2 = np.power(gamma*omega,2)
return p*np.power(omega0,2) / (np.sqrt(p1+p2))
parameters = sp.curve_fit(amplitude_func, frequency, amplitude)
p_fit = parameters[0][0]
omega0_fit = parameters[0][1]
gamma_fit = parameters[0][2]
print(parameters)
fit=[]
for i in frequency:
fit.append(amplitude_func(i, p_fit, omega0_fit, gamma_fit))
Right now I’m only having access to Python interpreter via ssh. So I can not ensure whether the plot and fit make sense or not.
Splitting the function to vars p1 and p2 was only for debugging purposes. It could be merged back together.