def update(self, measurement, measurement_matrix):
        """
        After each :py:meth:`predict`, this method may be called repeatedly to
        provide additional measurements for each time step.

        Args:
            measurement (MultivariateNormal): Measurement for this
                time step with specified mean and covariance.
            measurement_matrix (array): Measurement matrix for this measurement.

        """
        # Sanitise input arguments
        measurement_matrix = np.atleast_2d(measurement_matrix)
        expected_meas_mat_shape = (measurement.mean.shape[0], self.state_length)
        if measurement_matrix.shape != expected_meas_mat_shape:
            raise ValueError("Measurement matrix is wrong shape ({}). " \
                    "Expected: {}".format(
                        measurement_matrix.shape, expected_meas_mat_shape))

        # Add measurement list
        self.measurements[-1].append(measurement)
        self.measurement_matrices[-1].append(measurement_matrix)

        # "Prior" in this case means "before we've updated with this
        # measurement".
        prior = self.posterior_state_estimates[-1]

        # Compute Kalman gain
        innovation = measurement.mean - measurement_matrix.dot(prior.mean)
        innovation_cov = measurement_matrix.dot(prior.cov).dot(
            measurement_matrix.T)
        innovation_cov += measurement.cov
        kalman_gain = prior.cov.dot(measurement_matrix.T).dot(
            np.linalg.inv(innovation_cov))

        # Update estimates
        post = self.posterior_state_estimates[-1]
        self.posterior_state_estimates[-1] = MultivariateNormal(
            mean=post.mean + kalman_gain.dot(innovation),
            cov=post.cov - kalman_gain.dot(measurement_matrix).dot(prior.cov)
        )