Best way to calculate a long product of numbers between 0 and 1

Hi

I have these two variables:

l_h0 = (x(1,1)^(mb(1,1)+mb(1,5))*x(1,2)^(mb(1,2)+mb(1,6))*x(1,3)^(mb(1,3)+mb(1,7))*x(1,4)^(mb(1,4)+mb(1,8)))*...
       (x(2,1)^(mb(2,1)+mb(2,5))*x(2,2)^(mb(2,2)+mb(2,6))*x(2,3)^(mb(2,3)+mb(2,7))*x(2,4)^(mb(2,4)+mb(2,8)))*...
       (x(3,1)^(mb(3,1)+mb(3,5))*x(3,2)^(mb(3,2)+mb(3,6))*x(3,3)^(mb(3,3)+mb(3,7))*x(3,4)^(mb(3,4)+mb(3,8)))*...
       (x(4,1)^(mb(4,1)+mb(4,5))*x(4,2)^(mb(4,2)+mb(4,6))*x(4,3)^(mb(4,3)+mb(4,7))*x(4,4)^(mb(4,4)+mb(4,8)));
%
l_h1 = (y(1,1)^mb(1,1)*y(1,2)^mb(1,2)*y(1,3)^mb(1,3)*y(1,4)^mb(1,4))*...
       (y(2,1)^mb(2,1)*y(2,2)^mb(2,2)*y(2,3)^mb(2,3)*y(2,4)^mb(2,4))*...
       (y(3,1)^mb(3,1)*y(3,2)^mb(3,2)*y(3,3)^mb(3,3)*y(3,4)^mb(3,4))*...
       (y(4,1)^mb(4,1)*y(4,2)^mb(4,2)*y(4,3)^mb(4,3)*y(4,4)^mb(4,4))*...
       (y(1,5)^mb(1,5)*y(1,6)^mb(1,6)*y(1,7)^mb(1,7)*y(1,8)^mb(1,8))*...
       (y(2,5)^mb(2,5)*y(2,6)^mb(2,6)*y(2,7)^mb(2,7)*y(2,8)^mb(2,8))*...
       (y(3,5)^mb(3,5)*y(3,6)^mb(3,6)*y(3,7)^mb(3,7)*y(3,8)^mb(3,8))*...
       (y(4,5)^mb(4,5)*y(4,6)^mb(4,6)*y(4,7)^mb(4,7)*y(4,8)^mb(4,8));

Where both x, y, and mb are matrices. (I want to do a Neyman-Pearson test, where x is one matrix of likelihoods and y is another, but this is not important)

The thing is, both x and y are probabilities, that is, contain elements between 0 and 1, hence the calculations become really small. I want to divide these variables afterwards, doing l_h0/l_h1.

Which is the best way to adress this problem in terms of the most precision? I have tried using log in each one (inside, like log(x(1,1)..), then using exp() to 'cancel' the logs and get the 'normal' division again. Also I have tried to sum the log of each line, to avoid MATLAB doing the whole product with rounding and then calculating the log, and this two techniques give me different results when evaluating a lot of samples. Moreover, none of them convince me. Any ideas on how to do it in a more suitable way?

Thank you in advance.

Answer 1

[Walter Roberson]

It looks to me like

mb1_4 = mb(:,1:4);
mb5_8 = mb(:,5:8);
mbsum = mb1_4 + mb5_8;

log_l_h0 = sum( log(x(:)) .* mbsum(:) );
log_l_h1 = sum( log(y(:)) .* mb1_4(:) );
h0_over_h1 = exp(log_1_h0 - log_1_h1);

[Juan Gambetta]

Thank you. Seems obvious after seeing it.