error computation in multi layered perceptron
How do we get that $\delta_{ok}$?

You calculate the gradient of the network . Have a look at "Tom Mitchel: Machine Learning" if you want to see it in detail. In short, your weight update rule is

$$w \gets w + \Delta w$$
with the $j$-th component of the weight vector update being
$$\Delta w^{(j)} = - \eta \frac{\partial E}{\partial w^{(j)}}$$
where $\eta \in \mathbb{R}_+$ is the learning rate and $E$ is the error of your network.

$\delta_{ok}$ is just $\frac{\partial E}{\partial w^{(o,k)}}$. So I guess $o$ is an output neuron and $k$ a neuron of the last hidden layer.

Where does the chain rule play a role?

The chain rule is applied to compute the gradient. This is where the "backpropagation" comes from. You first calculate the gradient of the last weights, then the weights before, ... This is done by applying the chain rule to the error of the network.

Note

The weight initialization is not only small, but it has to be different for the different weights of one layer. Typically one chooses (pseudo)random weights.

See: Glorotm, Benigo: Understanding the difficulty of training deep feedforward neural networks