code: def d3logpdf_dlink3(self, link_f, y, Y_metadata=None):
          
          c = np.zeros_like(y)
          if Y_metadata is not None and 'censored' in Y_metadata.keys():
              c = Y_metadata['censored']
  
          val = np.log(y) - link_f
          val_scaled = val/np.sqrt(self.variance)
          val_scaled2 = val/self.variance
          a = (1 - stats.norm.cdf(val_scaled))
          uncensored = 0
          censored = c *( 2*np.exp(-3*(val**2)/(2*self.variance)) / ((a**3)*(2*np.pi*self.variance)**(3/2.))
                          - val*np.exp(-(val**2)/self.variance)/ ( (a**2)*np.pi*self.variance**2)
                          - val*np.exp(-(val**2)/self.variance)/ ( (a**2)*2*np.pi*self.variance**2)
                          - np.exp(-(val**2)/(2*self.variance))/ ( a*(self.variance**(1.50))*np.sqrt(2*np.pi))
                          + (val**2)*np.exp(-(val**2)/(2*self.variance))/ ( a*np.sqrt(2*np.pi*self.variance)*self.variance**2 ) )
          d3pdf_dlink3 = uncensored + censored
          return d3pdf_dlink3
docs: Gradient of the log-likelihood function at y given f, w.r.t shape parameter
  
          .. math::
  
          :param inv_link_f: latent variables link(f)
          :type inv_link_f: Nx1 array
          :param y: data
          :type y: Nx1 array
          :param Y_metadata: includes censoring information in dictionary key 'censored'
          :returns: derivative of likelihood evaluated at points f w.r.t variance parameter
          :rtype: float