Multiply Gaussian (mean1, var1) with (mean2, var2) and return the
    results as a tuple (mean, var, scale_factor).

    Strictly speaking the product of two Gaussian PDFs is a Gaussian
    function, not Gaussian PDF. It is, however, proportional to a Gaussian
    PDF. `scale_factor` provides this proportionality constant

    Parameters
    ----------
    mean1 : scalar
         mean of first Gaussian

    var1 : scalar
         variance of first Gaussian

    mean2 : scalar
         mean of second Gaussian

    var2 : scalar
         variance of second Gaussian

    Returns
    -------
    mean : scalar
        mean of product

    var : scalar
        variance of product

    scale_factor : scalar
        proportionality constant


    Examples
    --------
    >>> mul(1, 2, 3, 4)
    (1.6666666666666667, 1.3333333333333333)

    References
    ----------
    Bromily. "Products and Convolutions of Gaussian Probability Functions",
    Tina Memo No. 2003-003.
    http://www.tina-vision.net/docs/memos/2003-003.pdf

def mul_pdf(mean1, var1, mean2, var2):
    """
    Multiply Gaussian (mean1, var1) with (mean2, var2) and return the
    results as a tuple (mean, var, scale_factor).

    Strictly speaking the product of two Gaussian PDFs is a Gaussian
    function, not Gaussian PDF. It is, however, proportional to a Gaussian
    PDF. `scale_factor` provides this proportionality constant

    Parameters
    ----------
    mean1 : scalar
         mean of first Gaussian

    var1 : scalar
         variance of first Gaussian

    mean2 : scalar
         mean of second Gaussian

    var2 : scalar
         variance of second Gaussian

    Returns
    -------
    mean : scalar
        mean of product

    var : scalar
        variance of product

    scale_factor : scalar
        proportionality constant


    Examples
    --------
    >>> mul(1, 2, 3, 4)
    (1.6666666666666667, 1.3333333333333333)

    References
    ----------
    Bromily. "Products and Convolutions of Gaussian Probability Functions",
    Tina Memo No. 2003-003.
    http://www.tina-vision.net/docs/memos/2003-003.pdf
    """

    mean = (var1*mean2 + var2*mean1) / (var1 + var2)
    var = 1. / (1./var1 + 1./var2)

    S = math.exp(-(mean1 - mean2)**2 / (2*(var1 + var2))) / \
                 math.sqrt(2 * math.pi * (var1 + var2))

    return mean, var, S