def predict(self, control=None, control_matrix=None,
                process_matrix=None, process_covariance=None):
        """
        Predict the next *a priori* state mean and covariance given the last
        posterior. As a special case the first call to this method will
        initialise the posterior and prior estimates from the
        *initial_state_estimate* and *initial_covariance* arguments passed when
        this object was created. In this case the *process_matrix* and
        *process_covariance* arguments are unused but are still recorded in the
        :py:attr:`.process_matrices` and :py:attr:`.process_covariances`
        attributes.

        Args:
            control (array or None): If specified, the control input for this
                predict step.
            control_matrix (array or None): If specified, the control matrix to
                use for this time step.
            process_matrix (array or None): If specified, the process matrix to
                use for this time step.
            process_covariance (array or None): If specified, the process
                covariance to use for this time step.

        """
        # Sanitise arguments
        if process_matrix is None:
            process_matrix = self._defaults['process_matrix']

        if process_covariance is None:
            process_covariance = self._defaults['process_covariance']

        if control_matrix is None:
            control_matrix = self._defaults['control_matrix']

        if len(self.prior_state_estimates) == 0:
            # Special case: first call
            self.prior_state_estimates.append(self._initial_state_estimate)
        else:
            # Usual case
            process_matrix = as_square_array(process_matrix)
            process_covariance = as_square_array(process_covariance)
            if process_matrix.shape[0] != process_covariance.shape[0]:
                raise ValueError("Process matrix and noise have incompatible " \
                    "shapes: {} vs {}".format(
                        process_matrix.shape, process_covariance.shape))

            if control_matrix is not None:
                control_matrix = np.atleast_2d(control_matrix)
            if control is not None:
                control = np.atleast_1d(control)

            # Update state mean and covariance
            prev_posterior_mean = self.posterior_state_estimates[-1].mean
            prev_posterior_cov = self.posterior_state_estimates[-1].cov

            prior_mean = process_matrix.dot(prev_posterior_mean)
            if control is not None:
                prior_mean += control_matrix.dot(control)

            prior_cov = process_matrix.dot(prev_posterior_cov).dot(
                process_matrix.T) + process_covariance

            self.prior_state_estimates.append(
                MultivariateNormal(mean=prior_mean, cov=prior_cov))

        # Record transition matrix
        self.process_matrices.append(process_matrix)
        self.process_covariances.append(process_covariance)

        # Append empty list to measurements for this time step
        self.measurements.append([])
        self.measurement_matrices.append([])

        # Seed posterior estimates with the prior one.
        self.posterior_state_estimates.append(self.prior_state_estimates[-1])