Nice. Recreating these methods in simple code for yourself is definitely the way to check you understand it. This demo looks nice, clean, and straightforward. (Although I'd rename or comment variables x and y, or give some sort of guidance on what way around the weight matrices are within the code itself.)

It's also worth checking out existing neural net code-bases to see what tricks they have. The fine details usually aren't in papers, and they're not all in the text-books either.

The first potential problem that jumped out at me in this code was the initialization:

    self.weights = [np.array([0])] + [np.random.randn(y, x)
            for y, x in zip(sizes[1:], sizes[:-1])]

If the number of units in a layer is H, the typical size of the input into the layer above will be √(H). For large H, the sigmoid will usually saturate, and the gradients will underflow to zero, making it impossible to learn anything. There are some tricks to avoid the numerical problems, but even if you avoid numerical-underflow, things probably aren't going to work well.

I'd multiply those initial weights by a small constant divided by the square-root of weights going into the same neuron. For multiple layers you might consider layer-by-layer pre-training. For other architectures, like recurrent nets, definitely find a reference on how to do the initialization.

PS I would definitely add a test routine to check that the gradients from back-propagation agree with a finite difference approximation. It's so easy to get gradient code wrong, and it's so easy to test.