function se = get_BF_SE(RAB_B,RBA_A,BF)
    % Use delta method to compute approximate standard error of BFhat
    n = size(RAB_B,1);
    a = RAB_B./(RAB_B + BF);
    b = RBA_A./(1+BF*RBA_A);
    abar = mean(a);
    bbar = mean(b);
    aprime = -RAB_B./((RAB_B + BF).^2);
    bprime = -b.^2;
    aprimebar = mean(aprime);
    bprimebar = mean(bprime);
    var_a = std(a)^2;
    var_b = std(b)^2;
    cov_ab = 0;
    var_abar = var_a/n;
    var_bbar = var_b/n;
    cov_abar = cov_ab/n;

    var_bhat = ((BF^2)*var_abar - 2*BF*cov_abar + var_bbar);
    tmp = (abar+BF*aprimebar-bprimebar)^2;
    var_bhat = var_bhat/tmp;
    se = sqrt(var_bhat);

end