Trazando un límite de decisión que separa 2 clases usando la gráfica de puntos de Matplotlib

Su pregunta es más complicada que una simple ttwig: necesita dibujar el contorno que maximice la distancia entre clases. Afortunadamente, es un campo bien estudiado, especialmente para el aprendizaje automático de SVM.

El método más sencillo es descargar el módulo scikit-learn , que proporciona muchos métodos geniales para dibujar límites: http://scikit-learn.org/stable/modules/svm.html

Código:

 # -*- coding: utf-8 -*- import numpy as np import matplotlib from matplotlib import pyplot as plt import scipy from sklearn import svm mu_vec1 = np.array([0,0]) cov_mat1 = np.array([[2,0],[0,2]]) x1_samples = np.random.multivariate_normal(mu_vec1, cov_mat1, 100) mu_vec1 = mu_vec1.reshape(1,2).T # to 1-col vector mu_vec2 = np.array([1,2]) cov_mat2 = np.array([[1,0],[0,1]]) x2_samples = np.random.multivariate_normal(mu_vec2, cov_mat2, 100) mu_vec2 = mu_vec2.reshape(1,2).T fig = plt.figure() plt.scatter(x1_samples[:,0],x1_samples[:,1], marker='+') plt.scatter(x2_samples[:,0],x2_samples[:,1], c= 'green', marker='o') X = np.concatenate((x1_samples,x2_samples), axis = 0) Y = np.array([0]*100 + [1]*100) C = 1.0 # SVM regularization parameter clf = svm.SVC(kernel = 'linear', gamma=0.7, C=C ) clf.fit(X, Y) 

Diagtwig lineal (tomado de http://scikit-learn.org/stable/auto_examples/svm/plot_svm_margin.html )

 w = clf.coef_[0] a = -w[0] / w[1] xx = np.linspace(-5, 5) yy = a * xx - (clf.intercept_[0]) / w[1] plt.plot(xx, yy, 'k-') 

Gráfica multilínea (tomada de http://scikit-learn.org/stable/auto_examples/svm/plot_iris.html )

 C = 1.0 # SVM regularization parameter clf = svm.SVC(kernel = 'rbf', gamma=0.7, C=C ) clf.fit(X, Y) h = .02 # step size in the mesh # create a mesh to plot in x_min, x_max = X[:, 0].min() - 1, X[:, 0].max() + 1 y_min, y_max = X[:, 1].min() - 1, X[:, 1].max() + 1 xx, yy = np.meshgrid(np.arange(x_min, x_max, h), np.arange(y_min, y_max, h)) # Plot the decision boundary. For that, we will assign a color to each # point in the mesh [x_min, m_max]x[y_min, y_max]. Z = clf.predict(np.c_[xx.ravel(), yy.ravel()]) # Put the result into a color plot Z = Z.reshape(xx.shape) plt.contour(xx, yy, Z, cmap=plt.cm.Paird) 

Implementación

Si desea implementarlo usted mismo, necesita resolver la siguiente ecuación cuadrática:

El artículo de wikipedia.

Desafortunadamente, para los límites no lineales como el que se dibuja, es un problema difícil que se basa en un truco del núcleo, pero no existe una solución clara.

En función de la forma en que haya escrito decision_boundary , querrá usar la función de contour , como Joe señaló anteriormente. Si solo desea la línea de límite, puede dibujar un contorno único en el nivel 0:

 f, ax = plt.subplots(figsize=(7, 7)) c1, c2 = "#3366AA", "#AA3333" ax.scatter(*x1_samples.T, c=c1, s=40) ax.scatter(*x2_samples.T, c=c2, marker="D", s=40) x_vec = np.linspace(*ax.get_xlim()) ax.contour(x_vec, x_vec, decision_boundary(x_vec, mu_vec1, mu_vec2), levels=[0], cmap="Greys_r") 

Que hace:

Puedes crear tu propia ecuación para el límite:

donde tienes que encontrar las posiciones x0 y y0 , así como las constantes ai y bi para la ecuación del radio. Entonces, tienes 2*(n+1)+2 variables. Usar scipy.optimize.leastsq es sencillo para este tipo de problema.

El código que se adjunta a continuación genera el residuo para el leastsq penalización de los puntos que superan el límite del tamaño. El resultado para su problema, obtenido con:

 x, y = find_boundary(x2_samples[:,0], x2_samples[:,1], n) ax.plot(x, y, '-k', lw=2.) x, y = find_boundary(x1_samples[:,0], x1_samples[:,1], n) ax.plot(x, y, '--k', lw=2.) 

utilizando n=1 :

utilizando n=2 :

usng n=5 :

utilizando n=7 :

 import numpy as np from numpy import sin, cos, pi from scipy.optimize import leastsq def find_boundary(x, y, n, plot_pts=1000): def sines(theta): ans = np.array([sin(i*theta) for i in range(n+1)]) return ans def cosines(theta): ans = np.array([cos(i*theta) for i in range(n+1)]) return ans def residual(params, x, y): x0 = params[0] y0 = params[1] c = params[2:] r_pts = ((x-x0)**2 + (y-y0)**2)**0.5 thetas = np.arctan2((y-y0), (x-x0)) m = np.vstack((sines(thetas), cosines(thetas))).T r_bound = m.dot(c) delta = r_pts - r_bound delta[delta>0] *= 10 return delta # initial guess for x0 and y0 x0 = x.mean() y0 = y.mean() params = np.zeros(2 + 2*(n+1)) params[0] = x0 params[1] = y0 params[2:] += 1000 popt, pcov = leastsq(residual, x0=params, args=(x, y), ftol=1.e-12, xtol=1.e-12) thetas = np.linspace(0, 2*pi, plot_pts) m = np.vstack((sines(thetas), cosines(thetas))).T c = np.array(popt[2:]) r_bound = m.dot(c) x_bound = x0 + r_bound*cos(thetas) y_bound = y0 + r_bound*sin(thetas) return x_bound, y_bound 

Esas fueron algunas sugerencias, muchas gracias por su ayuda! Terminé resolviendo la ecuación analíticamente y esta es la solución con la que terminé (solo quiero publicarlo para futuras referencias)

Y el código se puede encontrar aquí.

EDITAR:

También tengo una función de conveniencia para trazar regiones de decisión para clasificadores que implementan un método de fit y predict , por ejemplo, los clasificadores en scikit-learn, que es útil si la solución no se puede encontrar analíticamente. Una descripción más detallada de cómo funciona se puede encontrar aquí .

Solo resolví un problema muy similar con un enfoque diferente (búsqueda de raíz) y quería publicar esta alternativa como respuesta aquí para referencia futura:

  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: 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) #g1 = discr_func(x, y, cov_mat=cov_mat1, mu_vec=mu_vec_1) #g2 = discr_func(x, y, cov_mat=cov_mat2, mu_vec=mu_vec_2) x_est50 = list(np.arange(-6, 6, 0.1)) y_est50 = [] for i in x_est50: y_est50.append(scipy.optimize.bisect(lambda y: discr_func(i, y, cov_mat=cov_est_1, mu_vec=mu_est_1) - \ discr_func(i, y, cov_mat=cov_est_2, mu_vec=mu_est_2), -10,10)) y_est50 = [float(i) for i in y_est50] 

Aquí está el resultado: (azul el caso cuadrático, rojo el caso lineal (varianzas iguales)

Sé que esta pregunta ha sido respondida de una manera muy analítica. Solo quería compartir un posible ‘hack’ al problema. Es difícil de manejar pero hace el trabajo.

Comience por construir una cuadrícula de malla del área 2d y luego, basándose en el clasificador, simplemente construya un mapa de clase de todo el espacio. Posteriormente, detecte los cambios en la decisión tomada en la fila y almacene los puntos de los bordes en una lista y distribuya los puntos.

 def disc(x): # returns the class of the point based on location x = [x,y] temp = 0.5 + 0.5*np.sign(disc0(x)-disc1(x)) # disc0() and disc1() are the discriminant functions of the respective classes return 0*temp + 1*(1-temp) num = 200 a = np.linspace(-4,4,num) b = np.linspace(-6,6,num) X,Y = np.meshgrid(a,b) def decColor(x,y): temp = np.zeros((num,num)) print x.shape, np.size(x,axis=0) for l in range(num): for m in range(num): p = np.array([x[l,m],y[l,m]]) #print p temp[l,m] = disc(p) return temp boundColorMap = decColor(X,Y) group = 0 boundary = [] for x in range(num): group = boundColorMap[x,0] for y in range(num): if boundColorMap[x,y]!=group: boundary.append([X[x,y],Y[x,y]]) group = boundColorMap[x,y] boundary = np.array(boundary) 

Ejemplo de límite de decisión para un clasificador gaussiano bivariado simple

Me gusta la biblioteca mglearn para dibujar límites de decisión. Aquí hay un ejemplo del libro “Introducción al aprendizaje automático con Python” de A. Mueller:

 fig, axes = plt.subplots(1, 3, figsize=(10, 3)) for n_neighbors, ax in zip([1, 3, 9], axes): clf = KNeighborsClassifier(n_neighbors=n_neighbors).fit(X, y) mglearn.plots.plot_2d_separator(clf, X, fill=True, eps=0.5, ax=ax, alpha=.4) mglearn.discrete_scatter(X[:, 0], X[:, 1], y, ax=ax) ax.set_title("{} neighbor(s)".format(n_neighbors)) ax.set_xlabel("feature 0") ax.set_ylabel("feature 1") axes[0].legend(loc=3) 

Si desea utilizar scikit learn, puede escribir su código así:

 import numpy as np import pandas as pd import matplotlib.pyplot as plt from sklearn.linear_model import LogisticRegression # read data data = pd.read_csv('ex2data1.txt', header=None) X = data[[0,1]].values y = data[2] # use LogisticRegression log_reg = LogisticRegression() log_reg.fit(X, y) # Coefficient of the features in the decision function. (from theta 1 to theta n) parameters = log_reg.coef_[0] # Intercept (aka bias) added to the decision function. (theta 0) parameter0 = log_reg.intercept_ # Plotting the decision boundary fig = plt.figure(figsize=(10,7)) x_values = [np.min(X[:, 1] -5 ), np.max(X[:, 1] +5 )] # calcul y values y_values = np.dot((-1./parameters[1]), (np.dot(parameters[0],x_values) + parameter0)) colors=['red' if l==0 else 'blue' for l in y] plt.scatter(X[:, 0], X[:, 1], label='Logistics regression', color=colors) plt.plot(x_values, y_values, label='Decision Boundary') plt.show() 

ver: Construir una logística-regresión-con-Scikit-aprender