code: def d2logpdf_dlink2(self, inv_link_f, y, Y_metadata=None):
          
          #d2logpdf_dlink2 = -y/(inv_link_f**2) - (1-y)/((1-inv_link_f)**2)
          #d2logpdf_dlink2 = np.where(y, -1./np.square(inv_link_f), -1./np.square(1.-inv_link_f))
          arg = np.where(y==1, inv_link_f, 1.-inv_link_f)
          ret =  -1./np.square(np.clip(arg, 1e-9, 1e9))
          if np.any(np.isinf(ret)):
              stop
          return ret
docs: Hessian at y, given inv_link_f, w.r.t inv_link_f the hessian will be 0 unless i == j
          i.e. second derivative logpdf at y given inverse link of f_i and inverse link of f_j  w.r.t inverse link of f_i and inverse link of f_j.
  
  
          .. math::
              \\frac{d^{2}\\ln p(y_{i}|\\lambda(f_{i}))}{d\\lambda(f)^{2}} = \\frac{-y_{i}}{\\lambda(f)^{2}} - \\frac{(1-y_{i})}{(1-\\lambda(f))^{2}}
  
          :param inv_link_f: latent variables inverse link of f.
          :type inv_link_f: Nx1 array
          :param y: data
          :type y: Nx1 array
          :param Y_metadata: Y_metadata not used in bernoulli
          :returns: Diagonal of log hessian matrix (second derivative of log likelihood evaluated at points inverse link of 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 inverse link of f_i not on inverse link of f_(j!=i)