All:
I am optimizing a likelihood function which contains really small probabilities (on the order of 10^-35) and want to use Octave to help me find optimal parameters that maximize the likelihood function. Using sqp on a function with such small values appears not to do anything and returns the guess that I give it. After getting some advice on where to seek help in the IRC channel I posted this apparent issue to the
bug tracker. It was closed as not being a bug but I was advised to post the issue to the mailing list in an effort to improve Octave. The exact inputs I can use to recreate the situation from a fresh start of Octave is in the post script. I should mention that I'm running Octave 4.2.2 on Windows 10 home and when I installed Octave, I was warned that it might behave incorrectly on Windows 10. I have some background with mathematics but almost none with programing.
The suggested change of taking the log of the likelihood function works. As I can change the tolerance for the sqp function to something that should still work with the unmodified function, this led me to believe that sqp would work. If sqp cannot work on functions that take on such small values then maybe that should be noted in the
documentation. In any event, it should be tested to see if indeed this is an issue and if it is an issue on other operating systems before changing the documentation. If desired, I could test sqp on Windows 10 to figure out what the true minimum tolerance is for sqp. I don't have time at the moment to do so.
I hope this is useful to someone else. I don't know what more details could be requested but I would be happy to provide them or if I need to clarify my thoughts for you'all.
In Gratitude,
Benjamin A Schwab
Post Script:
>> ttt = @(x) x(1)^2+x(2)^2
ttt =
@(x) x (1) ^ 2 + x (2) ^ 2
>> sqp([70,40],ttt,[],[],[-500,-500],[500,500],200,1e-40)
ans =
0
0
>> ggg=@(x) -factorial(27)*prod(1/sqrt(2*pi*x(2)^2).*exp(-([1,2,3,7,20,26,35,52,58,60,64,69,72,77,84,90,90,101,104,111,121,121,132,139,141,155,15
9].-x(1)).^2./(2*x(2)^2)))
ggg =
@(x) -factorial (27) * prod (1 / sqrt (2 * pi * x (2) ^ 2) .* exp (-([1, 2, 3, 7, 20, 26, 35, 52, 58, 60, 64, 69, 72, 77, 84, 90, 90, 101, 104, 11
1, 121, 121, 132, 139, 141, 155, 159] - x (1)) .^ 2 ./ (2 * x (2) ^ 2)))
>> ggg([70,40])
ans = -2.9010e-035
>> ggg([77,47])
ans = -1.1955e-034
>> ggg([77.543,47.707])
ans = -1.2047e-034
>> sqp([70,40],ggg,[],[],[-500,0],[500,500],200,1e-40)
ans =
70
40
>> bbb=@(x) -log(factorial(27)*prod(1/sqrt(2*pi*x(2)^2).*exp(-([1,2,3,7,20,26,35,52,58,60,64,69,72,77,84,90,90,101,104,111,121,121,132,139,141,15
5,159].-x(1)).^2./(2*x(2)^2))))
bbb =
@(x) -log (factorial (27) * prod (1 / sqrt (2 * pi * x (2) ^ 2) .* exp (-([1, 2, 3, 7, 20, 26, 35, 52, 58, 60, 64, 69, 72, 77, 84, 90, 90, 101, 10
4, 111, 121, 121, 132, 139, 141, 155, 159] - x (1)) .^ 2 ./ (2 * x (2) ^ 2))))
>> bbb([70,40])
ans = 79.525
>> bbb([77,47])
ans = 78.109
>> bbb([77.543,47.707])
ans = 78.102
>> sqp([70,40],bbb,[],[],[-500,0],[500,500])
ans =
77.556
47.691
>>