#!/usr/bin/env python # coding: utf-8 # Nonlinear Support Vector Machine import numpy as np import matplotlib.pyplot as plt import random from scipy.optimize import minimize from skimage import measure import itertools from math import sin, cos, pi d = 3 C = 5. def monomials(x, d): '''Return an array of all products of degree <= d composed from the coordinates of vector x''' # m = [1.] m = [] for l in range(1, d+1): for i in itertools.combinations_with_replacement(x, l): v = np.prod(i) m.append(v) return np.array(m) def f(x): '''Graph of this function separates two classes of points''' return 3.*sin(0.5*x) + 1. # xmin = -10.; xmax = 10.; # ymin = -10.; ymax = 10.; m = 100 class1 = [] class2 = [] for i in range(m): x1 = random.normalvariate(0., 4.) y1 = f(x1) + random.normalvariate(3., 2.) class1.append([x1, y1]) x2 = random.normalvariate(0., 4.) y2 = f(x2) - random.normalvariate(3., 2.) class2.append([x2, y2]) x = np.array(class1 + class2) # print("x:", x) y = np.array([1.]*m + [-1.]*m) # print("y:", y) # %matplotlib inline # get_ipython().run_line_magic('matplotlib', 'inline') ax = plt.gca() ax.set_aspect("equal") # ax.set_xlim(-10., 10.) # ax.set_ylim(-10., 10.) ax.scatter( [xx[0] for xx in x], [xx[1] for xx in x], color = [("red" if yy > 0 else "green") for yy in y] ) xext = np.array([monomials(xx, d) for xx in x]) k = len(xext[0]) print("Polynomial degree =", d) print("Dimension of extended space =", k) # Shuffle the training base perm = np.random.permutation(2*m) x = x[perm] xext = xext[perm] y = y[perm] def hingeLoss(x): return max(1. - x, 0.) def lossFunction(wb): '''Loss function to be minimized in Support Vector Machive''' n = len(xext) w = wb[:-1] b = wb[-1] s = 0. for i in range(n): c = y[i]*(w @ xext[i] - b) l = hingeLoss(c) s += l err = w@w + (C/n)*s # print("w:", w, "b:", b, "err:", err) return err wb = np.array([0.]*k + [0]) res = minimize(lossFunction, wb) # print("res;", res) w = res.x[:-1] b = res.x[-1] print("Parameters of linear classifier in extended space:") print("w =", w) print("b =", b) def classifier(x): '''Linear classifier constructed by SVM method''' xe = monomials(x, d) return w@xe - b # Draw the isoline defined by the equation # classifier(x) = 0 xmin = min([xx[0] for xx in x]) - 1. xmax = max([xx[0] for xx in x]) + 1. ymin = min([xx[1] for xx in x]) - 1. ymax = max([xx[1] for xx in x]) + 1. dx = 0.1; dy = dx xSteps = int((xmax - xmin)/dx) ySteps = int((ymax - ymin)/dy) def xcoord(j): return xmin + j*dx def ycoord(i): return ymin + i*dy # Matrix of function values a = [[0.]*xSteps for i in range(ySteps)] a = np.array(a) for i in range(ySteps): y = ycoord(i) for j in range(xSteps): x = xcoord(j) v = classifier(np.array([x, y])) a[i, j] = v contours = measure.find_contours(a, 0.) for contour in contours: ax.plot( [xcoord(c[1]) for c in contour], [ycoord(c[0]) for c in contour], color="blue" ) plt.show()