Using a combination of all, logical_andand operator <=, one can express the basic idea in 1 line.
import random
import numpy as np
from matplotlib import pyplot
points = [(random.random(), random.random()) for i in range(100)]
bx1, bx2 = sorted([random.random(), random.random()])
by1, by2 = sorted([random.random(), random.random()])
pts = np.array(points)
ll = np.array([bx1, by1])
ur = np.array([bx2, by2])
inidx = np.all(np.logical_and(ll <= pts, pts <= ur), axis=1)
inbox = pts[inidx]
outbox = pts[np.logical_not(inidx)]
rect = np.array([[bx1, by1], [bx1, by2], [bx2, by2], [bx2, by1], [bx1, by1]])
pyplot.plot(inbox[:, 0], inbox[:, 1], 'rx',
outbox[:, 0], outbox[:, 1], 'bo',
rect[:, 0], rect[:, 1], 'g-')
pyplot.show()
