matplotlib - How to plot classification contours on 2 features out of 4 in python -

i have written neural network classify iris dataset , want plot decision boundary , contours.

when select 2 features (2nd , 4th columns) out of 4 , feed network, can plot , looks this:

the problem when feed 4 features network, can't plot contours.

to resolve this, tried train network on 4 features, replace columns 1 , 3 , plot remaining data, doesn't give meaningful plot.

what should do?

import matplotlib.pyplot plt import numpy np import sklearn.datasets   iris = sklearn.datasets.load_iris()  # extract 2nd , 4th columns/features dataset x = iris.data[:, [1, 3]] y = iris.target  # m = number of test cases # n = number of input features # k = number of output labels # l = number of layers # s[l] = number of units in each layer  # x = [m, n] input matrix # y = [m, k] output matrix  # a[l] = activations each layer # w[l] = weights each layer # b[l] = biases each layer  # d[l] = deltas each layer # dw[l] = weight gradients each layer # db[l] = bias gradients each layer  m = y.size k = np.unique(y).size  y = np.zeros((m, k)) k in range(k):     y[:, k] = (k == y)  l = 4 s = [2, 9, 6, k]  = [none] * l  w = [none] * l  # w[l-1] = none b = [none] * l  # b[l-1] = none  d = [none] * l  # d[0] = none dw = [none] * l  # dw[l-1] = none db = [none] * l  # dw[l-1] = none  eta = 0.001 lambda = 0.001  np.random.seed(0) l in range(l-1):     w[l] = np.random.randn(s[l], s[l + 1]) / np.sqrt(s[l])     b[l] = np.random.randn(1, s[l + 1]) # --------------------------------------------------------------------------   ax = plt.gca()  h = 0.01 x_min, x_max = x[:, 0].min() - .5, x[:, 0].max() + .5 y_min, y_max = x[:, 1].min() - .5, x[:, 1].max() + .5 xx, yy = np.meshgrid(np.arange(x_min, x_max, h), np.arange(y_min, y_max, h)) # --------------------------------------------------------------------------   def classify(x):     a[0] = x      l in range(1, l):         z = np.dot(a[l-1], w[l-1]) + b[l-1]         a[l] = 1 / (1 + np.exp(-z))      return np.argmax(a[l-1], axis=1)   def updateplot():     z = classify(np.c_[xx.ravel(), yy.ravel()])     z = z.reshape(xx.shape)      ax.cla()     ax.contourf(xx, yy, z, c=y, alpha=0.75)     ax.scatter(x[:, 0], x[:, 1], c=y, s=50)      text = 'iteration: %05d\ncost: %.7f' % (i, cost)     ax.text(x_max + 0.2, y_max + 0.2, text, ha='right')     plt.pause(0.000001)   in range(10001):     # ++++++++++++++++++++++++++++++++++++++++++++++++++     # forward feed     a[0] = x      l in range(1, l):         z = np.dot(a[l-1], w[l-1]) + b[l-1]         a[l] = 1 / (1 + np.exp(-z))     # --------------------------------------------------      # ++++++++++++++++++++++++++++++++++++++++++++++++++     # propagation     d[l-1] = np.copy(a[l-1])     d[l-1][range(m), y] -= 1     l in range(l-2, 0, -1):         dw[l] = (a[l].t).dot(d[l + 1]) + lambda * w[l]         db[l] = np.sum(d[l + 1], axis=0, keepdims=true)          d[l] = d[l + 1].dot(w[l].t) * (a[l] * (1 - a[l]))      dw[0] = (a[0].t).dot(d[1]) + lambda * w[0]     db[0] = np.sum(d[1], axis=0, keepdims=true)     # --------------------------------------------------      # ++++++++++++++++++++++++++++++++++++++++++++++++++     # gradient descent     l in range(l-1):         w[l] += -eta * dw[l]         b[l] += -eta * db[l]     # --------------------------------------------------      # ++++++++++++++++++++++++++++++++++++++++++++++++++     # cost calculation , plotting     if % 100 == 0:         loss = -np.sum(y * np.log(a[l-1]) + (1-y) * np.log(1 - a[l-1]))          reg = 0         l in range(l-1):             reg += np.sum(np.square(w[l]))         reg *= lambda / 2         cost = 1 / m * (loss + reg)          updateplot()     # -------------------------------------------------- plt.show()