def pu_score(y, yhat):
    """
    Returns a score used for PU learning as introduced in [LEE2003]_.

    :param y: true function values
    :param yhat: predicted function values
    :returns:
        .. math:: \\frac{P(\hat{y}=1 | y=1)^2}{P(\hat{y}=1)}

    y and yhat must be boolean vectors.

    Higher is better.

    .. [LEE2003] Wee Sun Lee and Bing Liu. Learning with positive and unlabeled examples
        using weighted logistic regression. In Proceedings of the Twentieth
        International Conference on Machine Learning (ICML), 2003.
    """
    num_pos = sum(y)
    p_pred_pos = float(sum(yhat)) / len(y)
    if p_pred_pos == 0:
        return 0.0
    tp = sum([all(x) for x in zip(y, yhat)])
    return tp * tp / (num_pos * num_pos * p_pred_pos)