code: def d2logpdf_dlink2(self, link_f, y, Y_metadata=None):
          
          hess = -special.polygamma(1, self.beta*link_f)*(self.beta**2)
          #old
          #return -self.gp_link.d2transf_df2(gp)*self.beta*np.log(obs) + special.polygamma(1,self.gp_link.transf(gp)*self.beta)*(self.gp_link.dtransf_df(gp)*self.beta)**2 + special.psi(self.gp_link.transf(gp)*self.beta)*self.gp_link.d2transf_df2(gp)*self.beta
          return hess
docs: Hessian at y, given link(f), w.r.t link(f)
          i.e. second derivative logpdf at y given link(f_i) and link(f_j)  w.r.t link(f_i) and link(f_j)
          The hessian will be 0 unless i == j
  
          .. math::
              \\frac{d^{2} \\ln p(y_{i}|\lambda(f_{i}))}{d^{2}\\lambda(f)} = -\\beta^{2}\\frac{d\\Psi(\\alpha_{i})}{d\\alpha_{i}}\\\\
              \\alpha_{i} = \\beta y_{i}
  
          :param link_f: latent variables link(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: Diagonal of hessian matrix (second derivative of likelihood evaluated at points f)
          :rtype: Nx1 array
  
          .. Note::
              Will return diagonal of hessian, since every where else it is 0, as the likelihood factorizes over cases
              (the distribution for y_i depends only on link(f_i) not on link(f_(j!=i))