@staticmethod
def _split(x, n_factors, n_variables, n_lower_diag):
    """
        Given a one-dimensional array, split it
        into a set of multi-dimensional loading,
        error variance, factor variance, and factor
        covariance matrices.

        Parameters
        ----------
        x : array-like
            The combined (X, 1) array to split.
        n_factors : int
            The number of factors.
        n_variables : int
            The number of variables.
        n_lower_diag : int
            The number of elements in the `factor_covs`
            array, which is equal to the lower.

        Returns
        -------
        loadings : array-like
            The loadings matrix (n_factors * n_variables)
        error_vars : array-like (n_variables * 1)
            The error variance array.
        factor_vars : array-like (n_factors * 1)
            The factor variance array.
        factor_covs : array-like (n_lower_diag * 1)
            The factor covariance array
        """
    loadings_ix = int(n_factors * n_variables)
    error_vars_ix = n_variables + loadings_ix
    factor_vars_ix = n_factors + error_vars_ix
    factor_covs_ix = n_lower_diag + factor_vars_ix
    return x[:loadings_ix].reshape((n_variables, n_factors), order='F'), x[
        loadings_ix:error_vars_ix].reshape((n_variables, 1), order='F'), x[
        error_vars_ix:factor_vars_ix].reshape((n_factors, 1), order='F'), x[
        factor_vars_ix:factor_covs_ix].reshape((n_lower_diag, 1), order='F')
