Drawing decision boundary of two multivariate gaussian in Python
Question: I will borrow the image from the following stack overflow question to help me
describing my problem: [Drawing decision boundary of two multivariate
gaussian](http://stackoverflow.com/questions/19576761/drawing-decision-
boundary-of-two-multivariate-gaussian)

![enter image description here](http://i.stack.imgur.com/KwOPY.png)

I have 2 classes with 2D points, and what I am interested in is the decision
boundary (or discriminant).

I have written the functions that return the result for the discriminant
functions (a float value) that allowed me to classify the samples into those 2
patterns.

if a sample point is e.g., x_i = [x, y]  
I can call the functions and

if `g1(x,y) > g2(x,y)` it's class 1, and vice versa if

`g1(x,y) <= g2(x,y)` it's class 2

So the decision boundary should be at `g1(x,y) == g2(x,y)`

* * *

**EDIT:**

Hope that an example is helpful:  
  
1) Let's assume I take 1 sample `x = [1, 2]` from the dataset  
  
  
2) Then I will call e.g.  
`g1(1,2)` \--> returns `0.345`  
`g2(1,2)` \--> returns `0.453`  
\--> sample x belongs to class 2, since `g2(1,2) > g1(1,2)`  
  
  
3) Now for the decision boundary, I have `g2(x,y) == g1(x,y)`,  
or  
`g1(x,y) - g2(x,y)` == 0  
  
  
4) I generate a range of `x` values, e.g., `1,2,3,4,5`, and now I want to find
the corresponding `y` values that yield `g1(x,y) - g2(x,y) == 0`  
  
  
5) then I can use these `x,y` pairs to plot the decision boundary

* * *

In the StackOverflow post that I linked above, the suggestion would be to

> you can simply plot the contour line of f(x,y) := pdf1(x,y) > pdf2(x,y). So
> you define function f to be 1 iff pdf1(x,y)>pdf2(x,y). This way the only
> contour will be placed along the curve where pdf1(x,y)==pdf2(x,y) which is
> the decision boundary (discriminant). If you wish to define "nice" function
> you can do it simply by setting f(x,y) = sgn( pdf1(x,y) - pdf2(x,y) ), and
> plotting its contour plot will result in exact same discriminant.

But how would I do it in Python and matplotlib, I am really lost setting up
the code to do that. I am grateful for any help!

EDIT:

A little bit more about the function `g()` itself:

    
    
    def discr_func(x, y, cov_mat, mu_vec):
        """
        Calculates the value of the discriminant function for a dx1 dimensional
        sample given covariance matrix and mean vector.
    
        Keyword arguments:
            x_vec: A dx1 dimensional numpy array representing the sample.
            cov_mat: dxd numpy array of the covariance matrix.
            mu_vec: dx1 dimensional numpy array of the sample mean.
    
        Returns a float value as result of the discriminant function.
    
        """
        x_vec = np.array([[x],[y]])
    
        W_i = (-1/2) * np.linalg.inv(cov_mat)
        assert(W_i.shape[0] > 1 and W_i.shape[1] > 1), 'W_i must be a matrix'
    
        w_i = np.linalg.inv(cov_mat).dot(mu_vec)
        assert(w_i.shape[0] > 1 and w_i.shape[1] == 1), 'w_i must be a column vector'
    
        omega_i_p1 = (((-1/2) * (mu_vec).T).dot(np.linalg.inv(cov_mat))).dot(mu_vec)
        omega_i_p2 = (-1/2) * np.log(np.linalg.det(cov_mat))
        omega_i = omega_i_p1 - omega_i_p2
        assert(omega_i.shape == (1, 1)), 'omega_i must be a scalar'
    
        g = ((x_vec.T).dot(W_i)).dot(x_vec) + (w_i.T).dot(x_vec) + omega_i
        return float(g)
    

And when I execute it, it would return a float, e.g.,

`discr_func(1, 2, cov_mat=cov_est_1, mu_vec=mu_est_1)`  
`-3.726426544537969`

if I didn't make a mistake, it should be this equation: ![enter image
description here](http://i.stack.imgur.com/KthRb.png)

Thank you a lot for the suggestion with the contour, however, I have problems
implementing it:

    
    
    import pylab as pl
    
    X, Y = np.mgrid[-6:6:100j, -6:6:100j]
    x = X.ravel()
    y = Y.ravel()
    
    p = (discr_func(x, y, cov_mat=cov_est_1, mu_vec=mu_est_1) -\
         discr_func(x, y, cov_mat=cov_est_2, mu_vec=mu_est_2)).reshape(X.shape)
    
    #pl.scatter(X_train[:, 0], X_train[:, 1])
    pl.contour(X, Y, p, levels=[0])
    
    ---------------------------------------------------------------------------
    ValueError                                Traceback (most recent call last)
    <ipython-input-192-28c1c8787237> in <module>()
          5 y = Y.ravel()
          6 
    ----> 7 p = (discr_func(x, y, cov_mat=cov_est_1, mu_vec=mu_est_1) -     discr_func(x, y, cov_mat=cov_est_2, mu_vec=mu_est_2)).reshape(X.shape)
          8 
          9 #pl.scatter(X_train[:, 0], X_train[:, 1])
    
    <ipython-input-184-fd2f8b7fad82> in discr_func(x, y, cov_mat, mu_vec)
         25     assert(omega_i.shape == (1, 1)), 'omega_i must be a scalar'
         26 
    ---> 27     g = ((x_vec.T).dot(W_i)).dot(x_vec) + (w_i.T).dot(x_vec) + omega_i
         28     return float(g)
    
    ValueError: objects are not aligned
    

My feeling is that the passing of the `.ravel()` list doesn't work well with
how I set up this function... any suggestions?


Answer: calculate the `g1(x, y) - g2(x, y)` on a `mgrid[]` and then draw the line by
`contour(..., levels=[0])`, here is an example. Since you did not post any
sample data and code, I use sklearn to generate the sample data. You only need
the code after `#plot code from here`:

    
    
    import numpy as np
    import pylab as pl
    from sklearn import mixture
    
    np.random.seed(0)
    C1 = np.array([[3, -2.7], [1.5, 2.7]])
    C2 = np.array([[1, 2.0], [-1.5, 1.7]])
    
    X_train = np.r_[
        np.random.multivariate_normal((-5, -5), C1, size=100),
        np.random.multivariate_normal((5, 5), C2, size=100),
    ]
    
    clf = mixture.GMM(n_components=2, covariance_type='full')
    clf.fit(X_train)
    
    #define g1(x, y) and g2(x, y)
    
    def g1(x, y):
        return clf.predict_proba(np.column_stack((x, y)))[:, 0]
    
    def g2(x, y):
        return clf.predict_proba(np.column_stack((x, y)))[:, 1]
    
    #plot code from here
    
    X, Y = np.mgrid[-15:15:100j, -15:15:100j]
    x = X.ravel()
    y = Y.ravel()
    
    p = (g1(x, y) - g2(x, y)).reshape(X.shape)
    
    pl.scatter(X_train[:, 0], X_train[:, 1])
    pl.contour(X, Y, p, levels=[0])
    

here is the output:

![enter image description here](http://i.stack.imgur.com/3Jzmp.png)

