code: def d3logpdf_dlink3(self, link_f, y, Y_metadata=None):
          
          # d3lik_dlink3 = (1. - self.beta) / (y - link_f)**3
  
          c = np.zeros_like(y)
          if Y_metadata is not None and 'censored' in Y_metadata.keys():
              c = Y_metadata['censored']
          # uncensored = (1-c)* ((y ** self.r) * np.exp(-link_f))
          # censored = c*np.exp(-link_f)*y**self.r
          uncensored = (1-c)*(-2/link_f**3+ 6*y**self.r/link_f**4)
          censored = c*6*y**self.r/link_f**4
  
          d3lik_dlink3 = uncensored + censored
          # d3lik_dlink3 = (y ** self.r) * np.exp(-link_f)
          return d3lik_dlink3
docs: Third order derivative log-likelihood function at y given link(f) w.r.t link(f)
  
          .. math::
              \\frac{d^{3} \\ln p(y_{i}|\lambda(f_{i}))}{d^{3}\\lambda(f)} = -\\beta^{3}\\frac{d^{2}\\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: third derivative of likelihood evaluated at points f
          :rtype: Nx1 array