Gradient of the hessian (d2logpdf_dlink2) 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 hessian evaluated at points f and f_j w.r.t variance parameter
        :rtype: Nx1 array

def d2logpdf_dlink2_dr(self, inv_link_f, y, Y_metadata=None):
        """
        Gradient of the hessian (d2logpdf_dlink2) 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 hessian evaluated at points f and f_j w.r.t variance parameter
        :rtype: Nx1 array
        """
        # c = Y_metadata['censored']

        # c = np.zeros((y.shape[0],))
        c = np.zeros_like(y)
        if Y_metadata is not None and 'censored' in Y_metadata.keys():
            c = Y_metadata['censored']
        link_f = inv_link_f
        y_link_f = y/link_f
        y_link_f_r = y_link_f**self.r
        log_y_link_f = np.log(y) - np.log(link_f)

        #In terms of link_f
        y_link_f_2r = y_link_f**(2*self.r)
        denom2 = (link_f**2)*(1 + y_link_f_r)**2
        denom3 = (link_f**2)*(1 + y_link_f_r)**3

        censored = c*(-((y_link_f_r + self.r + 1)*y_link_f_r)/denom2
                      -(self.r*(y_link_f_r + self.r + 1)*y_link_f_r*log_y_link_f)/denom2
                      -(self.r*y_link_f_r*(y_link_f_r*log_y_link_f + 1))/denom2
                      +(2*self.r*(y_link_f_r + self.r + 1)*y_link_f_2r*log_y_link_f)/denom3
                      )

        uncensored = (1-c)*(-(2*self.r*y_link_f_r + y_link_f_2r - 1)/denom2
                            -(self.r*(2*y_link_f_r + 2*self.r*y_link_f_r*log_y_link_f + 2*y_link_f_2r*log_y_link_f)/denom2)
                            +(2*self.r*(2*self.r*y_link_f_r + y_link_f_2r - 1)*y_link_f_r*log_y_link_f)/denom3
                            )
        d2logpdf_dlink2_dr = censored + uncensored

        return d2logpdf_dlink2_dr