Compute the "gradient" of the model for the current parameters

        The default implementation computes the gradients numerically using
        a first order forward scheme. For better performance, this method
        should be overloaded by subclasses. The output of this function
        should be an array where the first dimension is ``full_size``.

def compute_gradient(self, *args, **kwargs):
        """
        Compute the "gradient" of the model for the current parameters

        The default implementation computes the gradients numerically using
        a first order forward scheme. For better performance, this method
        should be overloaded by subclasses. The output of this function
        should be an array where the first dimension is ``full_size``.

        """
        _EPS = 1.254e-5
        vector = self.get_parameter_vector()
        value0 = self.get_value(*args, **kwargs)
        grad = np.empty([len(vector)] + list(value0.shape), dtype=np.float64)
        for i, v in enumerate(vector):
            vector[i] = v + _EPS
            self.set_parameter_vector(vector)
            value = self.get_value(*args, **kwargs)
            vector[i] = v
            self.set_parameter_vector(vector)
            grad[i] = (value - value0) / _EPS
        return grad