Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

I have a task where I want to classify patterns from 2 classes where the samples are drawn from a bivariate Gaussian distribution. I use the 2 discriminant functions ($g_1$ and $g_2$) to classify the patterns, which is no problem. The decision rule would be simply
decide class $\omega_1$ if $g_1()$ > $g_2()$, and $\omega_2$ otherwise.

The code to do that would be very simple (I put my Python implementation at the bottom of this post if you are interested. Now, my problem is to plot the decision boundary, therefore I need to set $g_1()$ = $g_2()$, or $g_1()$ - $g_2() = 0$ all the parameters except $\pmb x$ would be known, but now I have to find $x_2$ for a range of $x_1$ value that yield $g_1()$ - $g_2() = 0$ in order to plot the boundary.

My understanding would be that I therefore have to rearrange the equation $g_1()$ - $g_2() = 0$ so that I have $x_2 = ... $ I tried to solve it analytically, but with no success... It would be awesome if you could help me with that, and I am also very open to alternative approaches. One idea was to minimize the function $g_1()$ - $g_2() = 0$ with an algorithm, but it also proofed to be more difficult than I thought.

I am hoping that you can help me with this problem. Thank you!

$ \Rightarrow g_1(\pmb{x}) = \pmb{x}^{\,t} - \frac{1}{2} \Sigma_1^{-1} \pmb{x} + \bigg( \Sigma_1^{-1} \pmb{\mu}_{\,1}\bigg)^t + \bigg( -\frac{1}{2} \pmb{\mu}_{\,1}^{\,t} \Sigma_{1}^{-1} \pmb{\mu}_{\,1} -\frac{1}{2} ln(|\Sigma_1|)\bigg) \\ \quad g_2(\pmb{x}) = \pmb{x}^{\,t} - \frac{1}{2} \Sigma_2^{-1} \pmb{x} + \bigg( \Sigma_2^{-1} \pmb{\mu}_{\,2}\bigg)^t + \bigg( -\frac{1}{2} \pmb{\mu}_{\,2}^{\,t} \Sigma_{2}^{-1} \pmb{\mu}_{\,2} -\frac{1}{2} ln(|\Sigma_2|)\bigg) $

$ \pmb{x} = \bigg[ \begin{array}{c} x_1 \\ x_2 \\ \end{array} \bigg] $

$\pmb{\mu} = \bigg[ \begin{array}{c} \mu_1 \\ \mu_2 \\ \end{array} \bigg] $

$ \Sigma = \begin{pmatrix} \lambda_{11} \quad \lambda_{12} \\ \lambda_{21} \quad \lambda_{22} \end{pmatrix}$

def discriminant_function(x_vec, 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: 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.

    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)

import operator

def classify_data(x_vec, g, mu_vecs, cov_mats):
    Classifies an input sample into 1 out of x classes determined by
    maximizing the discriminant function g_i().

    Keyword arguments:
        x_vec: A dx1 dimensional numpy array representing the sample.
        g: The discriminant function.
        mu_vecs: A list of mean vectors as input for g.
        cov_mats: A list of covariance matrices as input for g.

    Returns a tuple (g_i()_value, class label).

    assert(len(mu_vecs) == len(cov_mats)), 'Number of mu_vecs and cov_mats must be equal.'

    g_vals = []
    for m,c in zip(mu_vecs, cov_mats): 
        g_vals.append(g(x_vec, mu_vec=m, cov_mat=c))

    max_index, max_value = max(enumerate(g_vals), key=operator.itemgetter(1))
    return (max_value, max_index + 1)

import prettytable

classification_dict, error = empirical_error(all_samples, [1,2], classify_data, [discriminant_function,\
        [mu_est_1, mu_est_2],
        [cov_est_1, cov_est_2]])

labels_predicted = ['w{} (predicted)'.format(i) for i in [1,2]]
labels_predicted.insert(0,'test dataset')

train_conf_mat = prettytable.PrettyTable(labels_predicted)
for i in [1,2]:
    a, b = [classification_dict[i][j] for j in [1,2]]
    # workaround to unpack (since Python does not support just '*a')
    train_conf_mat.add_row(['w{} (actual)'.format(i), a, b])
print('Empirical Error: {:.2f} ({:.2f}%)'.format(error, error * 100))

| test dataset | w1 (predicted) | w2 (predicted) |
|   w1 (actual)    |       49       |       1        |
|   w2 (actual)    |       1        |       49       |
Empirical Error: 0.02 (2.00%)
share|cite|improve this question
You could try Newton's method to solve for $x_2$ numerically. – daw Apr 22 '14 at 11:39

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Browse other questions tagged or ask your own question.