fsolve not working properly with random sampling

Question

I am trying to use fsolve in python to find a value of a parameter for which an integral is equal to 1. I'm evaluating the integral through monte carlo sampling of my population. However fsolve never moves away from the initial guess.

In detail here is my code:

#Setting our Np number of pairs for our population, as well as other parameters Gamma and m
Np=1000
Gamma=20/100
m=5

#Keep the same initialization Step
Normals=[]
Hats=[]


for i in range(Np):
    #Creating our Omegas that are complex and uniformly distributed
    x_Omega=np.random.uniform(0,100)
    y_Omega=np.random.uniform(0,100)
    Omega=complex(x_Omega,y_Omega)
    
    #Creating our hs that are complex and uniformly distributed
    x_h=np.random.uniform(0,200)
    y_h=np.random.uniform(0,100)
    H=complex(x_h,y_h)
    
    #Creating our pairs of Omegas and Hs
    pair=(Omega,H)
    Normals.append(pair)
    
    #Creating our Omega hats that are complex and uniformly distributed
    x_Omega_Hat=np.random.uniform(0,100)
    y_Omega_Hat=np.random.uniform(0,100)
    Omega_Hat=complex(x_Omega_Hat,y_Omega_Hat)
    
    #Creating our h hats that are complex and uniformly distributed
    x_h_hat=np.random.uniform(0,200)
    y_h_hat=np.random.uniform(0,100)
    H_hat=complex(x_h_hat,y_h_hat)
    
    #Creating our Hat pairs
    pair=(Omega_Hat, H_hat)
    Hats.append(pair)
    
#Setting the s distributions:
s_dist=[0]*500
s_options=[0]*500
for s in range(500):
    s_dist[s]=poisson.pmf(s, m/Gamma)
    s_options[s]=s
    
#Defining the function
def F(Lambda, Hats, n):
    Integral=0
    for i in range(n):
        s=np.random.choice(s_options, p=s_dist)
        Selected=random.sample(Hats, s)
        Numerator=sum(item[1] for item in Selected)
        Denominator=1j*Lambda-sum(item[0] for item in Selected)
        Integral+=(Numerator/Denominator)**2
    return (Integral/n) -1 

The key part is the function at the end, in which I use Monte Carlo simulation to evaluate the Integral.

Then I use fsolve as such:

# Define the function that returns a system of equations, where the imaginary part is set to 0
def equations(lambda_real_imag):
    Lambda = lambda_real_imag[0] + 1j * lambda_real_imag[1]
    Criteria=F(Lambda, Hats, n)
    result = [Criteria.real, Criteria.imag]

    # Print progress information
    print("Equations: {}, Lambda: {}".format(result, Lambda))
        
    return result

# Initial guess for Lambda_real and Lambda_imag
initial_guess = [1.0, 1.0]

# Use fsolve to find the root of the system of equations
solution = fsolve(equations, initial_guess)

However, the results never converge to 0, and the value of the parameter Lambda that fsolve uses never strays from the initial guess. Results with random sampling

Interestingly if I take the random sampling out of the definition of the function and define the 'selected' sample at the start then it works perfectly, and a value for Lambda is found. I.e. if I write instead:

s=np.random.choice(s_options, p=s_dist)
Selected=random.sample(Hats, s)
def F(Lambda, Hats, n):
    Integral=0
    for i in range(n):
        Numerator=sum(item[1] for item in Selected)
        Denominator=1j*Lambda-sum(item[0] for item in Selected)
        Integral+=(Numerator/Denominator)**2
    return (Integral/n) -1 

However this removes the random sampling and is totally invalid.