def unit_variance(N_tot, N_au, N_alpha, N_beta):
    """Compute the variance of E_{lpha eta}(t)"""
    if N_alpha == 0 or N_beta == 0:
        return float('nan')
    if N_au == 0:
        return 0
    else:
        var = 1.0 / (N_alpha * N_beta) * (N_tot / N_au - 1
            ) ** 2 + 1.0 / N_alpha * (N_tot / N_au - 1) + 1.0 / N_beta * (
            N_tot / N_au - 1)
    return var
