(mu, var, with_derivatives=False, **kwargs)
return logP, dlogPdMu, dlogPdSigma, dlogPdMudMu
Computes the probability of every given point to be the minimum
    based on the EPMGP[1] algorithm.

    [1] Cunningham, P. Hennig, and S. Lacoste-Julien.
    Gaussian probabilities and expectation propagation.
    under review. Preprint at arXiv, November 2011.

    Parameters
    ----------
    M: np.ndarray(N,)
        Mean value of each of the N points.

    V: np.ndarray(N, N)
        Covariance matrix for all points

    with_derivatives: bool
        If true than also the gradients are computed

    Returns
    -------
    np.ndarray(N,1)
        pmin distribution
Computes the probability of every given point to be the minimum
    based on the EPMGP[1] algorithm.
Computes the probability of every given point to be the minimum based on the EPMGP [ 1 ] algorithm .
def joint_min(mu, var, with_derivatives=False, **kwargs):
    """
    Computes the probability of every given point to be the minimum
    based on the EPMGP[1] algorithm.

    [1] Cunningham, P. Hennig, and S. Lacoste-Julien.
    Gaussian probabilities and expectation propagation.
    under review. Preprint at arXiv, November 2011.

    Parameters
    ----------
    M: np.ndarray(N,)
        Mean value of each of the N points.

    V: np.ndarray(N, N)
        Covariance matrix for all points

    with_derivatives: bool
        If true than also the gradients are computed

    Returns
    -------
    np.ndarray(N,1)
        pmin distribution
    """

    logP = np.zeros(mu.shape)
    D = mu.shape[0]
    if with_derivatives:
        dlogPdMu = np.zeros((D, D))
        dlogPdSigma = np.zeros((D, int(0.5 * D * (D + 1))))
        dlogPdMudMu = np.zeros((D, D, D))
    for i in range(mu.shape[0]):

        # logP[k] ) self._min_faktor(mu, var, 0)
        a = min_faktor(mu, var, i)

        logP[i] = next(a)
        if with_derivatives:
            dlogPdMu[i, :] = next(a).T
            dlogPdMudMu[i, :, :] = next(a)
            dlogPdSigma[i, :] = next(a).T

    logP[np.isinf(logP)] = -500
    # re-normalize at the end, to smooth out numerical imbalances:
    logPold = logP
    Z = np.sum(np.exp(logPold))
    maxLogP = np.max(logP)
    s = maxLogP + np.log(np.sum(np.exp(logP - maxLogP)))
    s = maxLogP if np.isinf(s) else s

    logP = logP - s
    if not with_derivatives:
        return logP

    dlogPdMuold = dlogPdMu
    dlogPdSigmaold = dlogPdSigma
    dlogPdMudMuold = dlogPdMudMu
    # adjust derivatives, too. This is a bit tedious.
    Zm = sum(np.rot90((np.exp(logPold) * np.rot90(dlogPdMuold, 1)), 3)) / Z
    Zs = sum(np.rot90((np.exp(logPold) * np.rot90(dlogPdSigmaold, 1)), 3)) / Z

    dlogPdMu = dlogPdMuold - Zm
    dlogPdSigma = dlogPdSigmaold - Zs

    ff = np.einsum('ki,kj->kij', dlogPdMuold, dlogPdMuold)
    gg = np.einsum('kij,k->ij', dlogPdMudMuold + ff, np.exp(logPold)) / Z
    Zij = Zm.T * Zm
    adds = np.reshape(-gg + Zij, (1, D, D))
    dlogPdMudMu = dlogPdMudMuold + adds
    return logP, dlogPdMu, dlogPdSigma, dlogPdMudMu
def joint_min ( mu , var , with_derivatives = False , * * kwargs ) : logP = np . zeros ( mu . shape ) D = mu . shape [ 0 ] if with_derivatives : dlogPdMu = np . zeros ( ( D , D ) ) dlogPdSigma = np . zeros ( ( D , int ( 0.5 * D * ( D + 1 ) ) ) ) dlogPdMudMu = np . zeros ( ( D , D , D ) ) for i in range ( mu . shape [ 0 ] ) : # logP[k] ) self._min_faktor(mu, var, 0) a = min_faktor ( mu , var , i ) logP [ i ] = next ( a ) if with_derivatives : dlogPdMu [ i , : ] = next ( a ) . T dlogPdMudMu [ i , : , : ] = next ( a ) dlogPdSigma [ i , : ] = next ( a ) . T logP [ np . isinf ( logP ) ] = - 500 # re-normalize at the end, to smooth out numerical imbalances: logPold = logP Z = np . sum ( np . exp ( logPold ) ) maxLogP = np . max ( logP ) s = maxLogP + np . log ( np . sum ( np . exp ( logP - maxLogP ) ) ) s = maxLogP if np . isinf ( s ) else s logP = logP - s if not with_derivatives : return logP dlogPdMuold = dlogPdMu dlogPdSigmaold = dlogPdSigma dlogPdMudMuold = dlogPdMudMu # adjust derivatives, too. This is a bit tedious. Zm = sum ( np . rot90 ( ( np . exp ( logPold ) * np . rot90 ( dlogPdMuold , 1 ) ) , 3 ) ) / Z Zs = sum ( np . rot90 ( ( np . exp ( logPold ) * np . rot90 ( dlogPdSigmaold , 1 ) ) , 3 ) ) / Z dlogPdMu = dlogPdMuold - Zm dlogPdSigma = dlogPdSigmaold - Zs ff = np . einsum ( 'ki,kj->kij' , dlogPdMuold , dlogPdMuold ) gg = np . einsum ( 'kij,k->ij' , dlogPdMudMuold + ff , np . exp ( logPold ) ) / Z Zij = Zm . T * Zm adds = np . reshape ( - gg + Zij , ( 1 , D , D ) ) dlogPdMudMu = dlogPdMudMuold + adds return logP , dlogPdMu , dlogPdSigma , dlogPdMudMu