u"""
        Modifies the factor with marginalized values.

        Let C(X,Y ; K, h, g) be some canonical form over X,Y where,

        k = [[K_XX, K_XY],        ;       h = [[h_X],
             [K_YX, K_YY]]                     [h_Y]]

        In this case, the result of the integration operation is a canonical
        from C (K', h', g') given by,

        .. math:: K' = K_{XX} - K_{XY} * {K^{-1}}_{YY} * K_YX
        .. math:: h' = h_X - K_{XY} * {K^{-1}}_{YY} * h_Y
        .. math:: g' = g + 0.5 * (|Y| * log(2*pi) - log(|K_{YY}|) + {h^T}_Y * K_{YY} * h_Y)

        Parameters
        ----------

        variables: list or array-like
                List of variables over which to marginalize.

        inplace: boolean
                If inplace=True it will modify the distribution itself,
                else would return a new distribution.

        Returns
        -------
        CanonicalDistribution or None :
                if inplace=True (default) returns None
                if inplace=False return a new CanonicalDistribution instance

        Examples
        --------
        >>> import numpy as np
        >>> from pgmpy.factors.continuous import CanonicalDistribution
        >>> phi = CanonicalDistribution(['X1', 'X2', 'X3'],
        ...                       np.array([[1, -1, 0], [-1, 4, -2], [0, -2, 4]]),
        ...                       np.array([[1], [4], [-1]]), -2)
        >>> phi.K
        array([[ 1, -1,  0],
                [-1,  4, -2],
                [ 0, -2,  4]])

        >>> phi.h
        array([[ 1],
                [ 4],
                [-1]])

        >>> phi.g
        -2

        >>> phi.marginalize(['X3'])

        >>> phi.K
        array([[ 1., -1.],
                [-1.,  3.]])

        >>> phi.h
        array([[ 1. ],
                [ 3.5]])

        >>> phi.g
        0.22579135

def marginalize(self, variables, inplace=True):
        u"""
        Modifies the factor with marginalized values.

        Let C(X,Y ; K, h, g) be some canonical form over X,Y where,

        k = [[K_XX, K_XY],        ;       h = [[h_X],
             [K_YX, K_YY]]                     [h_Y]]

        In this case, the result of the integration operation is a canonical
        from C (K', h', g') given by,

        .. math:: K' = K_{XX} - K_{XY} * {K^{-1}}_{YY} * K_YX
        .. math:: h' = h_X - K_{XY} * {K^{-1}}_{YY} * h_Y
        .. math:: g' = g + 0.5 * (|Y| * log(2*pi) - log(|K_{YY}|) + {h^T}_Y * K_{YY} * h_Y)

        Parameters
        ----------

        variables: list or array-like
                List of variables over which to marginalize.

        inplace: boolean
                If inplace=True it will modify the distribution itself,
                else would return a new distribution.

        Returns
        -------
        CanonicalDistribution or None :
                if inplace=True (default) returns None
                if inplace=False return a new CanonicalDistribution instance

        Examples
        --------
        >>> import numpy as np
        >>> from pgmpy.factors.continuous import CanonicalDistribution
        >>> phi = CanonicalDistribution(['X1', 'X2', 'X3'],
        ...                       np.array([[1, -1, 0], [-1, 4, -2], [0, -2, 4]]),
        ...                       np.array([[1], [4], [-1]]), -2)
        >>> phi.K
        array([[ 1, -1,  0],
                [-1,  4, -2],
                [ 0, -2,  4]])

        >>> phi.h
        array([[ 1],
                [ 4],
                [-1]])

        >>> phi.g
        -2

        >>> phi.marginalize(['X3'])

        >>> phi.K
        array([[ 1., -1.],
                [-1.,  3.]])

        >>> phi.h
        array([[ 1. ],
                [ 3.5]])

        >>> phi.g
        0.22579135
        """
        if not isinstance(variables, (list, tuple, np.ndarray)):
            raise TypeError("variables: Expected type list or array-like, "
                            "got type {var_type}".format(var_type=type(variables)))

        if not all([var in self.variables for var in variables]):
            raise ValueError("Variable not in scope.")

        phi = self if inplace else self.copy()

        # index_to_keep -> i vector
        index_to_keep = [self.variables.index(var) for var in self.variables
                         if var not in variables]
        # index_to_marginalize -> j vector
        index_to_marginalize = [self.variables.index(var) for var in variables]

        K_i_i = self.K[np.ix_(index_to_keep, index_to_keep)]
        K_i_j = self.K[np.ix_(index_to_keep, index_to_marginalize)]
        K_j_i = self.K[np.ix_(index_to_marginalize, index_to_keep)]
        K_j_j = self.K[np.ix_(index_to_marginalize, index_to_marginalize)]
        K_j_j_inv = np.linalg.inv(K_j_j)
        h_i = self.h[index_to_keep]
        h_j = self.h[index_to_marginalize]

        phi.variables = [self.variables[index] for index in index_to_keep]

        phi.K = K_i_i - np.dot(np.dot(K_i_j, K_j_j_inv), K_j_i)
        phi.h = h_i - np.dot(np.dot(K_i_j, K_j_j_inv), h_j)
        phi.g = self.g + 0.5 * (len(variables) * np.log(2 * np.pi) -
                                np.log(abs(np.linalg.det(K_j_j))) + np.dot(np.dot(h_j.T, K_j_j), h_j))[0][0]

        if not inplace:
            return phi