Sympy nsolve returns lambdify TypeErrors

I am trying to solve an overdefined system of equations nummerically with Sympy nsolve, but I keep getting lambdify TypeErrors.

import sympy as sp

f25 = 61900.
f50 = 66200.
f75 = 70700.
f90 = 77700.

#mu = sp.Symbol('mu')
#sigma = sp.Symbol('sigma')

mu = sp.var('mu')
sigma = sp.var('sigma')

def f(x):
    f = (1 / sp.sqrt(2 * sp.pi * sp.Pow(sigma,2))) * sp.euler(-1 * (sp.Pow(x - mu,2)) / (2 * sp.Pow(sigma,2)))
    return f

f1 = f(f25) - 0.25
f2 = f(f50) - 0.50
f3 = f(f75) - 0.75
f4 = f(f90) - 0.90

sp.nsolve((f1, f2, f3, f4), (mu, sigma), (66000, 10000))

returns

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
File ~\Anaconda3\lib\site-packages\mpmath\calculus\optimization.py:937, in findroot(ctx, f, x0, solver, tol, verbose, verify, **kwargs)
    936 try:
--> 937     fx = f(*x0)
    938     multidimensional = isinstance(fx, (list, tuple, ctx.matrix))

File <lambdifygenerated-5>:2, in _lambdifygenerated(mu, sigma)
      1 def _lambdifygenerated(mu, sigma):
----> 2     return ImmutableDenseMatrix([[mpf((1, 1, -2, 1)) + (mpf(1)/mpf(2))*sqrt(2)*euler(-mpf((0, 239475625, 3, 28))*(1 - mpf((0, 2384070316472963, -67, 52))*mu)**2/sigma**2)/(sqrt(pi)*sqrt(sigma**2))], [mpf((1, 1, -1, 1)) + (mpf(1)/mpf(2))*sqrt(2)*euler(-mpf((0, 68475625, 5, 27))*(1 - mpf((0, 139325861583909, -63, 47))*mu)**2/sigma**2)/(sqrt(pi)*sqrt(sigma**2))], [mpf((1, 3, -2, 2)) + (mpf(1)/mpf(2))*sqrt(2)*euler(-mpf((0, 312405625, 3, 29))*(1 - mpf((0, 8349304248355101, -69, 53))*mu)**2/sigma**2)/(sqrt(pi)*sqrt(sigma**2))], [mpf((1, 8106479329266893, -53, 53)) + (mpf(1)/mpf(2))*sqrt(2)*euler(-mpf((0, 377330625, 3, 29))*(1 - mpf((0, 3798557338215609, -68, 52))*mu)**2/sigma**2)/(sqrt(pi)*sqrt(sigma**2))]])

File ~\Anaconda3\lib\site-packages\mpmath\ctx_mp_python.py:346, in _constant.__call__(self, prec, dps, rounding)
    345 if dps: prec = dps_to_prec(dps)
--> 346 return self.context.make_mpf(self.func(prec, rounding))

File ~\Anaconda3\lib\site-packages\mpmath\libmp\libelefun.py:116, in def_mpf_constant.<locals>.f(prec, rnd)
    115 wp = prec + 20
--> 116 v = fixed(wp)
    117 if rnd in (round_up, round_ceiling):

File ~\Anaconda3\lib\site-packages\mpmath\libmp\libelefun.py:97, in constant_memo.<locals>.g(prec, **kwargs)
     96 if prec <= memo_prec:
---> 97     return f.memo_val >> (memo_prec-prec)
     98 newprec = int(prec*1.05+10)

TypeError: unsupported operand type(s) for >>: 'int' and 'mpf'

During handling of the above exception, another exception occurred:

TypeError                                 Traceback (most recent call last)
Cell In[3], line 24
     21 f3 = f(f75) - 0.75
     22 f4 = f(f90) - 0.90
---> 24 sp.nsolve((f1, f2, f3, f4), (mu, sigma), (66000, 10000))

File ~\Anaconda3\lib\site-packages\sympy\utilities\decorator.py:88, in conserve_mpmath_dps.<locals>.func_wrapper(*args, **kwargs)
     86 dps = mpmath.mp.dps
     87 try:
---> 88     return func(*args, **kwargs)
     89 finally:
     90     mpmath.mp.dps = dps

File ~\Anaconda3\lib\site-packages\sympy\solvers\solvers.py:3000, in nsolve(dict, *args, **kwargs)
   2998 J = lambdify(fargs, J, modules)
   2999 # solve the system numerically
-> 3000 x = findroot(f, x0, J=J, **kwargs)
   3001 if as_dict:
   3002     return [dict(zip(fargs, [sympify(xi) for xi in x]))]

File ~\Anaconda3\lib\site-packages\mpmath\calculus\optimization.py:940, in findroot(ctx, f, x0, solver, tol, verbose, verify, **kwargs)
    938     multidimensional = isinstance(fx, (list, tuple, ctx.matrix))
    939 except TypeError:
--> 940     fx = f(x0[0])
    941     multidimensional = False
    942 if 'multidimensional' in kwargs:

TypeError: _lambdifygenerated() missing 1 required positional argument: 'sigma'

I have tried to define my variables as both sp.Symbol and sp.var.

I have tried to bracket the variables to solve for: sp.nsolve((f1, f2, f3, f4), [(mu, sigma)], (66000, 10000)). But that just gives me a new TypeError: TypeError: nsolve expected exactly 1 guess vectors, got 2

Bracketing my input guess values yields yet another TypeError: sp.nsolve((f1, f2, f3, f4), [(mu, sigma)], [(66000, 10000)]). TypeError: cannot create mpf from (66000, 10000).

As I am sure you can tell I have all but given up, and are going at it brute force "monkeys and typewriters" style.

Can anyone offer up any advice?

>Solution :

Looking at your f(x), they appear to be some kind of normal statistical distribution. If so, you might want to replace sp.euler with sp.E ** , where sp.E is Euler’s number.