def find_bias(sf_stack, opt_sf):
    """
    Finds difference between avg(sf) and optimal sf for one 72x46 grid

    Parameters:
        sf_stack (numpy array) : nxMx72x46 array with inverted scale factors
                                 (n=number OSSEs and M=number of months)
        opt_sf   (numpy array) : Mx72x46 optimal scale factor array

    Returns:
        Mx46x76 numpy array of E(sf) - opt_sf

    NOTE:
     - we assume that the 0th index of sf_stack is the OSSE iterations
    """
    assert sf_stack.shape[1] == opt_sf.shape[0]  # same number of months
    assert sf_stack.shape[2] == 72
    assert opt_sf.shape[1] == 72

    # find mean of the given stack of sf draws
    sf_stack_avg = sf_stack.mean(axis=0)

    return sf_stack_avg - opt_sf