function  [BF,log_BF,se_BF,se_log_BF] = get_BF(L_A_A,L_B_A,L_A_B,L_B_B)
    % Get BF and its standard error
     % Form ratios
     RAB_B = exp(L_A_B - L_B_B);
     RBA_A = exp(L_B_A - L_A_A);
    
     BF = 1;
     jj = 0;
     iter = 0;
     while jj == 0
         iter = iter+1;
         BF_old = BF;
         num = RAB_B./(RAB_B + BF);
         den = RBA_A./(1+BF*RBA_A);
         BF = mean(num)/mean(den);
         dif = abs(BF-BF_old);
         % [iter BF]
         if dif < 0.00001
             jj = 1;
         end
         if iter > 1000
             error('no convergence')
         end
     end
     % Approximate SE for BF
     se_BF = get_BF_SE(RAB_B,RBA_A,BF);
     log_BF = log(BF);
     se_log_BF = se_BF/BF;
     
end