code: def normal_equation(x, y):
      
      num_samples = y.size
      x = np.insert(x, 0, np.ones(num_samples), axis=1)
  
      return pinv(x)*y
docs: X: matrix of features, one sample per row (without bias unit)
      y: values (continuous) corresponding to rows (samples) in X