code: def score(self, outcomes, modelparams, expparams, return_L=False):
          r
          
          if len(modelparams.shape) == 1:
              modelparams = modelparams[:, np.newaxis]
          
          # compute likelihood at central point
          L0 = self.likelihood(outcomes, modelparams, expparams)
          
          # allocate space for the score
          q = np.empty([self.n_modelparams, 
                        outcomes.shape[0], 
                        modelparams.shape[0], 
                        expparams.shape[0]])
          h_perturb = np.empty(modelparams.shape)
          
          # just loop over the model parameter as there usually won't be so many
          # of them that vectorizing would be worth the effort.
          for mp_idx in range(self.n_modelparams):
              h_perturb[:] = np.zeros(modelparams.shape)
              h_perturb[:, mp_idx] = self.h[mp_idx]
              # use the chain rule since taking the numerical derivative of a 
              # logarithm is unstable
              q[mp_idx, :] = (
                  self.likelihood(outcomes, modelparams + h_perturb, expparams) - 
                  self.likelihood(outcomes, modelparams - h_perturb, expparams)
              ) / (2 * self.h[mp_idx] * L0)
              
          
          if return_L:
              return q, L0
          else: 
              return q
docs: r"""
          Returns the numerically computed score of the likelihood 
          function, defined as:
          
          .. math::
          
              q(d, \vec{x}; \vec{e}) = \vec{\nabla}_{\vec{x}} \log \Pr(d | \vec{x}; \vec{e}).
              
          Calls are represented as a four-index tensor
          ``score[idx_modelparam, idx_outcome, idx_model, idx_experiment]``.
          The left-most index may be suppressed for single-parameter models.
          
          The numerical gradient is computed using the central difference method, 
          with step size given by the property `~ScoreMixin.h`.
          
          If return_L is True, both `q` and the likelihood `L` are returned as `q, L`.