def fit(self, X, design, nuisance=None, scan_onsets=None, coords=None,
            inten=None):
        """Compute the Bayesian RSA

        Parameters
        ----------
        X: numpy array, shape=[time_points, voxels]
            If you have multiple scans of the same participants that you
            want to analyze together, you should concatenate them along
            the time dimension after proper preprocessing (e.g. spatial
            alignment), and specify the onsets of each scan in scan_onsets.
        design: numpy array, shape=[time_points, conditions]
            This is the design matrix. It should only include the hypothetic
            response for task conditions. You should not include
            regressors for a DC component or motion parameters, unless you
            want to estimate their pattern similarity with response patterns
            to your task conditions. If you want to model head motion,
            you should include them in nuisance regressors.
            If you have multiple run, the design matrix
            of all runs should be concatenated along the time dimension,
            with every column for one condition across runs.
            For example, if you have 3 runs of experiment of one participant,
            with each run lasting 200 TR. And you have 4 conditions,
            then design should be a 600 x 4 numpy array.
        nuisance: optional, numpy array, shape=[time_points, nuisance_factors]
            The responses to these regressors will be marginalized out from
            each voxel, which means they are considered, but won't be assumed
            to share the same pseudo-SNR map with the design matrix.
            Therefore, the pseudo-SNR map will only reflect the
            relative contribution of design matrix to each voxel.
            You can provide time courses such as those for head motion
            to this parameter.
            Note that if auto_nuisance is set to True, the first
            n_nureg principal components of residual (excluding the response
            to the design matrix and the user-provided nuisance regressors
            and a constant baseline)
            will be included as additional nuisance regressor after the
            first round of fitting.
            If auto_nuisance is set to False, the nuisance regressors supplied
            by the users together with DC components will be used as
            nuisance time series.
            Please do not include time course of constant baseline in nuisance.
        scan_onsets: optional, numpy array, shape=[runs,]
            This specifies the indices of X which correspond to the onset
            of each scanning run. For example, if you have two experimental
            runs of the same subject, each with 100 TRs, then scan_onsets
            should be [0,100].
            If you do not provide the argument, the program will
            assume all data are from the same run.
            The effect of them is to make the inverse matrix
            of the temporal covariance matrix of noise block-diagonal.
        coords: optional, numpy array, shape=[voxels,3]
            This is the coordinate of each voxel,
            used for implementing Gaussian Process prior.
        inten: optional, numpy array, shape=[voxel,]
            This is the average fMRI intensity in each voxel.
            It should be calculated from your data without any preprocessing
            such as z-scoring. Because it should reflect
            whether a voxel is bright (grey matter) or dark (white matter).
            A Gaussian Process kernel defined on both coordinate and intensity
            imposes a smoothness prior on adjcent voxels
            but with the same tissue type. The Gaussian Process
            is experimental and has shown good performance on
            some visual datasets.
        """

        logger.info('Running Bayesian RSA')
        self.random_state_ = check_random_state(self.random_state)
        # setting random seed
        logger.debug('RandState set to {}'.format(self.random_state_))

        assert not self.GP_inten or (self.GP_inten and self.GP_space),\
            'You must speficiy GP_space to True'\
            'if you want to use GP_inten'

        # Check input data
        assert_all_finite(X)
        assert X.ndim == 2, 'The data should be 2-dimensional ndarray'

        assert np.all(np.std(X, axis=0) > 0),\
            'The time courses of some voxels do not change at all.'\
            ' Please make sure all voxels are within the brain'

        # check design matrix
        assert_all_finite(design)
        assert design.ndim == 2,\
            'The design matrix should be 2-dimensional ndarray'
        assert np.linalg.matrix_rank(design) == design.shape[1], \
            'Your design matrix has rank smaller than the number of'\
            ' columns. Some columns can be explained by linear '\
            'combination of other columns. Please check your design matrix.'
        assert np.size(design, axis=0) == np.size(X, axis=0),\
            'Design matrix and data do not '\
            'have the same number of time points.'
        assert self.rank is None or self.rank <= design.shape[1],\
            'Your design matrix has fewer columns than the rank you set'

        # Check the nuisance regressors.
        if nuisance is not None:
            assert_all_finite(nuisance)
            assert nuisance.ndim == 2,\
                'The nuisance regressor should be 2-dimensional ndarray'
            assert np.linalg.matrix_rank(nuisance) == nuisance.shape[1], \
                'The nuisance regressor has rank smaller than the number of'\
                'columns. Some columns can be explained by linear '\
                'combination of other columns. Please check your nuisance' \
                'regressors.'
            assert np.size(nuisance, axis=0) == np.size(X, axis=0), \
                'Nuisance regressor and data do not have the same '\
                'number of time points.'
        # check scan_onsets validity
        assert scan_onsets is None or\
            (np.max(scan_onsets) <= X.shape[0] and np.min(scan_onsets) >= 0),\
            'Some scan onsets provided are out of the range of time points.'

        # check the size of coords and inten
        if self.GP_space:
            logger.info('Fitting with Gaussian Process prior on log(SNR)')
            assert coords is not None and coords.shape[0] == X.shape[1],\
                'Spatial smoothness was requested by setting GP_space. '\
                'But the voxel number of coords does not match that of '\
                'data X, or voxel coordinates are not provided. '\
                'Please make sure that coords is in the shape of '\
                '[n_voxel x 3].'
            assert coords.ndim == 2,\
                'The coordinate matrix should be a 2-d array'
            if self.GP_inten:
                assert inten is not None and inten.shape[0] == X.shape[1],\
                    'The voxel number of intensity does not '\
                    'match that of data X, or intensity not provided.'
                assert np.var(inten) > 0,\
                    'All voxels have the same intensity.'
        if (not self.GP_space and coords is not None) or\
                (not self.GP_inten and inten is not None):
            logger.warning('Coordinates or image intensity provided'
                           ' but GP_space or GP_inten is not set '
                           'to True. The coordinates or intensity are'
                           ' ignored.')
        # Estimate the number of necessary nuisance regressors
        if self.auto_nuisance:
            if self.n_nureg is None:
                logger.info('number of nuisance regressors is determined '
                            'automatically.')
                run_TRs, n_runs = self._run_TR_from_scan_onsets(
                    X.shape[0], scan_onsets)
                ts_dc = self._gen_legendre(run_TRs, [0])
                _, ts_base, _ = self._merge_DC_to_base(
                    ts_dc, nuisance, False)
                ts_reg = np.concatenate((ts_base, design), axis=1)
                beta_hat = np.linalg.lstsq(ts_reg, X)[0]
                residuals = X - np.dot(ts_reg, beta_hat)
                self.n_nureg_ = np.max(
                    [1, Ncomp_SVHT_MG_DLD_approx(residuals,
                                                 self.nureg_zscore)])
                logger.info('Use {} nuisance regressors to model the spatial '
                            'correlation in noise.'.format(self.n_nureg_))
                self.n_nureg_ = np.int32(self.n_nureg_)
            else:
                self.n_nureg_ = self.n_nureg
            self.n_nureg_ = np.int32(self.n_nureg_)
        # Run Bayesian RSA
        # Note that we have a change of notation here. Within _fit_RSA_UV,
        # design matrix is named X and data is named Y, to reflect the
        # generative model that data Y is generated by mixing the response
        # X to experiment conditions and other neural activity.
        # However, in fit(), we keep the tradition of scikit-learn that
        # X is the input data to fit and y, a reserved name not used, is
        # the label to map to from X.

        if not self.GP_space:
            # If GP_space is not requested, then the model is fitted
            # without imposing any Gaussian Process prior on log(SNR^2)
            self.U_, self.L_, self.nSNR_, self.beta_, self.beta0_,\
                self._beta_latent_, self.sigma_, self.rho_, _, _, _,\
                self.X0_ = self._fit_RSA_UV(X=design, Y=X, X_base=nuisance,
                                            scan_onsets=scan_onsets)
        elif not self.GP_inten:
            # If GP_space is requested, but GP_inten is not, a GP prior
            # based on spatial locations of voxels will be imposed.
            self.U_, self.L_, self.nSNR_, self.beta_, self.beta0_,\
                self._beta_latent_, self.sigma_, self.rho_, \
                self.lGPspace_, self.bGP_, _, \
                self.X0_ = self._fit_RSA_UV(
                    X=design, Y=X, X_base=nuisance,
                    scan_onsets=scan_onsets, coords=coords)
        else:
            # If both self.GP_space and self.GP_inten are True,
            # a GP prior based on both location and intensity is imposed.
            self.U_, self.L_, self.nSNR_, self.beta_, self.beta0_,\
                self._beta_latent_, self.sigma_, self.rho_, \
                self.lGPspace_, self.bGP_, self.lGPinten_, self.X0_ = \
                self._fit_RSA_UV(X=design, Y=X, X_base=nuisance,
                                 scan_onsets=scan_onsets,
                                 coords=coords, inten=inten)

        self.C_ = utils.cov2corr(self.U_)
        self.design_ = design.copy()
        self._rho_design_, self._sigma2_design_ = \
            self._est_AR1(self.design_, same_para=True)
        self._rho_X0_, self._sigma2_X0_ = self._est_AR1(self.X0_)
        # AR(1) parameters of the design matrix and nuisance regressors,
        # which will be used in transform or score.

        # Finally, we fit a null model with the same setting except
        # that there is no response to X
        self.beta0_null_, self.sigma_null_, self.rho_null_, \
            self.X0_null_ = self._fit_null(Y=X, X_base=nuisance,
                                           scan_onsets=scan_onsets)
        self._rho_X0_null_, self._sigma2_X0_null_ =\
            self._est_AR1(self.X0_null_)
        return self
Compute the Bayesian RSA

        Parameters
        ----------
        X: numpy array, shape=[time_points, voxels]
            If you have multiple scans of the same participants that you
            want to analyze together, you should concatenate them along
            the time dimension after proper preprocessing (e.g. spatial
            alignment), and specify the onsets of each scan in scan_onsets.
        design: numpy array, shape=[time_points, conditions]
            This is the design matrix. It should only include the hypothetic
            response for task conditions. You should not include
            regressors for a DC component or motion parameters, unless you
            want to estimate their pattern similarity with response patterns
            to your task conditions. If you want to model head motion,
            you should include them in nuisance regressors.
            If you have multiple run, the design matrix
            of all runs should be concatenated along the time dimension,
            with every column for one condition across runs.
            For example, if you have 3 runs of experiment of one participant,
            with each run lasting 200 TR. And you have 4 conditions,
            then design should be a 600 x 4 numpy array.
        nuisance: optional, numpy array, shape=[time_points, nuisance_factors]
            The responses to these regressors will be marginalized out from
            each voxel, which means they are considered, but won't be assumed
            to share the same pseudo-SNR map with the design matrix.
            Therefore, the pseudo-SNR map will only reflect the
            relative contribution of design matrix to each voxel.
            You can provide time courses such as those for head motion
            to this parameter.
            Note that if auto_nuisance is set to True, the first
            n_nureg principal components of residual (excluding the response
            to the design matrix and the user-provided nuisance regressors
            and a constant baseline)
            will be included as additional nuisance regressor after the
            first round of fitting.
            If auto_nuisance is set to False, the nuisance regressors supplied
            by the users together with DC components will be used as
            nuisance time series.
            Please do not include time course of constant baseline in nuisance.
        scan_onsets: optional, numpy array, shape=[runs,]
            This specifies the indices of X which correspond to the onset
            of each scanning run. For example, if you have two experimental
            runs of the same subject, each with 100 TRs, then scan_onsets
            should be [0,100].
            If you do not provide the argument, the program will
            assume all data are from the same run.
            The effect of them is to make the inverse matrix
            of the temporal covariance matrix of noise block-diagonal.
        coords: optional, numpy array, shape=[voxels,3]
            This is the coordinate of each voxel,
            used for implementing Gaussian Process prior.
        inten: optional, numpy array, shape=[voxel,]
            This is the average fMRI intensity in each voxel.
            It should be calculated from your data without any preprocessing
            such as z-scoring. Because it should reflect
            whether a voxel is bright (grey matter) or dark (white matter).
            A Gaussian Process kernel defined on both coordinate and intensity
            imposes a smoothness prior on adjcent voxels
            but with the same tissue type. The Gaussian Process
            is experimental and has shown good performance on
            some visual datasets.