FIX: Nonconstant covariances handling

Hello!

Currently if you pass ``transition_covariance`` or ``observation_covariance`` as 3D arrays you'll get wrong results or a script will crash. I'm talking about ``KalmanFilter``, ``CholeskyKalmanFilter`` and ``BiermanKalmanFilter`` classes (I didn't touch unscented filters.)

In ``KalmanFilter`` one can find the following lines:
```Python
transition_covariance = _last_dims(transition_covariance, t - 1) 
observation_covariance = _last_dims(observation_covariance, t)
```
It overwrites references and then you are back to the case of a 2D constant matrix.

Square root filters in turn simply don't account for possible 3D arrays of covariances.

So I tried to fix all that. And also did quite a bit of code reformatting. The changes are not completely subtle (e.g I changed ``transition_covariance`` to ``transition_covariances`` in class intefraces and so on), but since this package is used by rather few people and in research environment, I don't think it's a big deal. The supplied tests are passing all right.

So let's work on it and improve pykalman. Please review my PR.