code: def dlogpdf_dlink(self, link_f, y, Y_metadata=None):
          
          grad = self.beta*np.log(self.beta*y) - special.psi(self.beta*link_f)*self.beta
          #old
          #return -self.gp_link.dtransf_df(gp)*self.beta*np.log(obs) + special.psi(self.gp_link.transf(gp)*self.beta) * self.gp_link.dtransf_df(gp)*self.beta
          return grad
docs: Gradient of the log likelihood function at y, given link(f) w.r.t link(f)
  
          .. math::
              \\frac{d \\ln p(y_{i}|\\lambda(f_{i}))}{d\\lambda(f)} = \\beta (\\log \\beta y_{i}) - \\Psi(\\alpha_{i})\\beta\\\\
              \\alpha_{i} = \\beta y_{i}
  
          :param link_f: latent variables (f)
          :type link_f: Nx1 array
          :param y: data
          :type y: Nx1 array
          :param Y_metadata: Y_metadata which is not used in gamma distribution
          :returns: gradient of likelihood evaluated at points
          :rtype: Nx1 array