What is causing "TypeError: only size-1 arrays can be converted to Python scalars least_squares" from implementation of scipy.optimize.least_squares

Question

I am attempting to fit a normal distribution to histogram values from seaborn. I am using the following code:

from scipy.optimize import least_squares
from math import exp


def model_normal_dist(x, u):
    return x[0]*exp(-((x[1]-u)/x[2])**2)

def model_residual(x,u,y):
    print(model_normal_dist(x, u) - y)
    return model_normal_dist(x, u) - y

u=np.array([h.xy[0]+.5*h.get_width() for h in sns.histplot(data=data,x='visitor_prop1').patches])
y=np.array([h.get_height() for h in  sns.histplot(data=data,x='visitor_prop1').patches])
x=np.array([10000,.2,1])


res_1 = least_squares(model_residual, x, args=(u, y), verbose=1)

where data is a dataframe with column, 'visitor_prop1' of float values. When I run this I get TypeError: only size-1 arrays can be converted to Python scalars. I've tried putting print statements in the above functions but they apparently never run.

I've testes that u and y are arrays of floats. An sample of that is:

print(u[0:30])
print(y[0:30])

[0.00109038 0.00327085 0.00545132 0.00763179 0.00981227 0.01199274
 0.01417321 0.01635368 0.01853415 0.02071462 0.02289509 0.02507556
 0.02725603 0.0294365  0.03161697 0.03379744 0.03597791 0.03815838
 0.04033886 0.04251933 0.0446998  0.04688027 0.04906074 0.05124121
 0.05342168 0.05560215 0.05778262 0.05996309 0.06214356 0.06432403]
[ 704.  658.  757.  754.  848.  920.  924.  980. 1048. 1124. 1134. 1300.
 1343. 1420. 1516. 1593. 1689. 1752. 1821. 1970. 1997. 2142. 2162. 2234.
 2409. 2602. 2624. 2715. 2914. 2927.]

I've checked that type(model_residual(x,u[5],y[5])) returns numpy.float64 which appears to be a python scalar What's causing this error?

Answer 1

This error was caused by math.exp() which only works on scalars. scipy.optimize.least_squares needs to be able to operate on arrays. The solution is to swap out math.exp() with numpy.exp():

def model_normal_dist(x, u):
    return x[0]*np.exp(-((x[1]-u)/x[2])**2)

or keep the function as is and change the import:

from numpy import exp

The print statements couldn't run because the error was in the function they had to evaluate to print.