def score(self, X, design, scan_onsets=None):
        """ After fit() is applied to the data of a group of participants,
            use the parameters estimated by fit() function to evaluate
            from some data of a set of participants to evaluate
            the log likelihood of some new data of the same participants
            given these estimated parameters.
            Design matrices of the same set of experimental
            conditions in the testing data should be provided, with each
            column corresponding to the same condition as that column
            in the design matrix of the training data.
            Unknown nuisance time series will be marginalized, assuming
            they follow the same spatial pattern as in the training
            data. The hypothetical response captured by the design matrix
            will be subtracted from data before the marginalization
            when evaluating the log likelihood. For null model,
            nothing will be subtracted before marginalization.

            There is a difference between the form of likelihood function
            used in fit() and score(). In fit(), the response amplitude
            beta to design matrix X and the modulation beta0 by nuisance
            regressor X0 are both marginalized, with X provided and X0
            estimated from data. In score(), posterior estimation of
            beta and beta0 from the fitting step are assumed unchanged
            in testing data; X is assumed given by the user,
            and X0 is marginalized.
            The logic underlying score() is to transfer
            as much as what we can learn from training data when
            calculating a likelihood score for testing data. This is done
            at the cost of using point estimation for beta and beta0.

            If you z-scored your data during fit step, you should
            z-score them for score function as well. If you did not
            z-score in fitting, you should not z-score here either.

        Parameters
        ----------
        X : List of 2-D arrays. For each item, shape=[time_points, voxels]
            fMRI data of new data of the same participants.
            The voxels of each participants should
            match those used in the fit() function. If data are z-scored
            (recommended) when fitting the model, data should be z-scored
            as well when calling transform()
        design : List of 2-D arrays. shape=[time_points, conditions] for each
            Each corresponds to one participant.
            Design matrices expressing the hypothetical response of
            the task conditions in data X.
        scan_onsets : List of 2-D arrays, shape=[#fMRI runs] for each
            Each array corresponds to one participant.
            Lists of indices corresponding to the onsets of
            scans in the data X.
            If not provided, data will be assumed
            to be acquired in a continuous scan.

        Returns
        -------
        ll: list, shape=[number of participants]
            The log likelihoods of the new data based on the model and its
            parameters fit to the training data.
            If data of some participants are not provided, the corresponding
            entry will be None.
        ll_null: list, shape=[number of participants]
            The log likelihood of the new data based on a null model
            which assumes the same as the full model for everything
            except for that there is no response to any of the
            task conditions.
        """
        X = self._check_data_GBRSA(X, for_fit=False)
        scan_onsets = self._check_scan_onsets_GBRSA(scan_onsets, X)
        design = self._check_design_GBRSA(design, X)
        assert len(X) == self.n_subj_
        ll = [None] * self.n_subj_
        ll_null = [None] * self.n_subj_
        for subj in np.arange(self.n_subj_):
            if X[subj] is not None:
                ll[subj] = self._score(
                    Y=X[subj], design=design[subj], beta=self.beta_[subj],
                    scan_onsets=scan_onsets[subj], beta0=self.beta0_[subj],
                    rho_e=self.rho_[subj], sigma_e=self.sigma_[subj],
                    rho_X0=self._rho_X0_[subj],
                    sigma2_X0=self._sigma2_X0_[subj])
                ll_null[subj] = self._score(
                    Y=X[subj], design=None, beta=None,
                    scan_onsets=scan_onsets[subj], beta0=self.beta0_[subj],
                    rho_e=self.rho_[subj], sigma_e=self.sigma_[subj],
                    rho_X0=self._rho_X0_[subj],
                    sigma2_X0=self._sigma2_X0_[subj])
        return ll, ll_null