def d2logpdf_dlink2_dvar(self, link_f, y, Y_metadata=None):
        """
        :param link_f: latent variables link(f)
        :type link_f: Nx1 array
        :param y: data
        :type y: Nx1 array
        :param Y_metadata: Y_metadata not used in gaussian
        :returns: derivative of log likelihood evaluated at points link(f) w.r.t variance parameter
        :rtype: Nx1 array
        """
        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 = (1-c)*( 1./(self.variance**2) )
        censored = c*( val*np.exp(-3*(val**2)/(2*self.variance) )/ ((a**3)*np.sqrt(8*np.pi**3)*self.variance**(5/2.))
                       + np.exp(-val**2/self.variance)/((a**2)*4*np.pi*self.variance**2)
                       - np.exp(-val**2/self.variance)*val**2 / ((a**2)*2*np.pi*self.variance**3)
                       + np.exp(-val**2/self.variance)/ ( (a**2)*4*np.pi*self.variance**2)
                       - np.exp(-val**2/ (2*self.variance))*val / ( a*np.sqrt(2*np.pi)*2*self.variance**(5/2.))
                       - np.exp(-val**2/self.variance)*(val**2) / ((a**2)*4*np.pi*self.variance**3)
                       - np.exp(-val**2/ (2*self.variance))*val/ (a*np.sqrt(2*np.pi)*self.variance**(5/2.))
                       + np.exp(-val**2/ (2*self.variance))*(val**3) / (a*np.sqrt(2*np.pi)*2*self.variance**(7/2.)) )
        dlik_hess_dsigma = uncensored + censored
        return dlik_hess_dsigma