I think there is a bug in KalmanFilter.compute, more specifically on line 522 of kalmanfilter.py. I noticed that when I tried to perform KalmanFilter.em in a scenario where observations consist of more than 1 variable, the code crashes due to broadcasting errors.

On line 522, result.filtered.gains is initialized as follows:
result.filtered.gains = np.empty((n_vars, n_measurements, n_states, n_states))
while I think it should be:
result.filtered.gains = np.empty((n_vars, n_measurements, n_states, n_obs))

Would it be possible to have a look, and hotfix if needed?
Many thanks in advance!