def estimate_fst_freqs_wright(p1, p2):
    """
    calculate mean fst based on called frequencies
    """
    p = (p1 + p2) / 2.0
    is_het = (p > 0) & (p < 1)
    vp = ((p1 - p) ** 2 + (p2 - p) ** 2) / 2.0
    fsts = vp[is_het] / (p * (1 - p))[is_het]
    return fsts.mean()
