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.
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().