Runtime error in log function is occurring when using minimize from scipy, how should I fix this?

Question

For context I'm using NumPy. I encountered the following warning message in my code

    RuntimeWarning: invalid value encountered in log
  model = ((-2*r) * np.log(1-s*np.sin(2*np.pi*t/p_orb-phi_orb))) + (-(R**2)/(c*D*9.461E15)*(np.cos(Beta)**2)*np.cos((4*np.pi*t)/(3.154E7)-2*Lambda))

Heres the code for context:

r = 9.85E-6
i = np.pi/2
s = np.sin(i)
N = 157788
t = np.linspace(0 , 9.461E7 , N)
p_orb = 2400
phi_orb = np.pi/3
R = 1.496E11
c = 299752498
D = 1000
Beta = 50*np.pi/180
Lambda = np.pi/4
sigmas = np.ones(N)*0.000001
data = ((-2*r) * np.log(1-s*np.sin(2*np.pi*t/p_orb-phi_orb))) + (-(R**2)/(c*D*9.461E15)*(np.cos(Beta)**2)*np.cos((4*np.pi*t)/(3.154E7)-2*Lambda)) + sigmas

def log_likelihood(theta , t , data , sigmas):
    r , s , D = theta
    model = ((-2*r) * np.log(1-s*np.sin(2*np.pi*t/p_orb-phi_orb))) + (-(R**2)/(c*D*9.461E15)*(np.cos(Beta)**2)*np.cos((4*np.pi*t)/(3.154E7)-2*Lambda))
    return -0.5 * sum(((data - model)**2)/(sigmas**2))

nll = lambda *args: -log_likelihood(*args)
initial = np.array([r , s , D])
soln = minimize(nll , initial , args = (t , data , sigmas))
r_ml, s_ml, D_ml = soln.x
print(r_ml , s_ml , D_ml)

Just running the code in my Jupyter notebook results in the initial values I gave to the minimize command along with the warning message. I checked for zeros and negatives using the following code:

array = 1-s*np.sin(2*np.pi*t/p_orb-phi_orb)
# Check for zeros
zeros_indices = array == 0
zeros_exist = np.any(zeros_indices)

# Check for negative values
negative_indices = array < 0
negatives_exist = np.any(negative_indices)

if zeros_exist:
    print("Array contains zero(s).")
else:
    print("Array does not contain any zero.")

if negatives_exist:
    print("Array contains negative value(s).")
else:
    print("Array does not contain any negative value.")

The array in this case is the argument of the log function. There doesn't seem to be negatives or zeros. I do not know what else to check for.

Answer 1

The problem is that sis one of your optimization variables, so that the minimize() function is free to try whatever value it likes in evaluating things like log(1-s.sin(at-b)). It is almost certain to try values either side of the expected answer of 1.

As soon as a value of s is tried that is above 1 then it is almost guaranteed (depending on the contents of the t[] array) that the argument of the logarithm will go negative (which is invalid).

You can stop the warning during running of your code by using the bounds= option to clamp s to a limited range. Thus the following will "work", in the sense that it will recover your original data parameters (including s=1):

soln = minimize(nll, initial, args=(t,data,sigmas), bounds=((None,None),(0.0,1.0),(None,None)) )

This produces (r,S,D) as

9.93687622297525e-06 1.0 1000.0

However, this is more by luck than judgement. Your discrete set of t values didn't quite produce a value with sin() equal to 1, so log(1-s.sin()) never failed.

I suggest that you DON'T try to minimise like this with model data generated from s=sin(pi/2)=1 and you DO put formal bounds on your optimisation parameters to prevent invalid arguments to functions like log().