Gradient of marginal likelihood of Gaussian Process w.r.t likelihood parameters with Laplace approximation

The derivation of gradient of the marginal likelihood w.r.t covariance function hyperparameters $$\theta$$ is given in http://www.gaussianprocess.org/gpml/chapters/RW5.pdf, page 125. However, the gradient w.r.t likelihood parameters (let's call them $$\sigma$$) is not derived.

I found the final equations for this gradient w.r.t $$\sigma$$ in the source code of the GPML toolbox by the book's authors and the GPStuff toolbox, as well as in this paper. However, I found no derivation. The equations are similar to those of the gradient w.r.t the covariance function hyperparameters $$\theta$$. But I don't understand how $$\frac{\partial \hat{f}}{\partial \sigma}$$ is derived (supposedly similar to $$\frac{\partial \hat{f}}{\partial \theta}$$ in Eq. (5.24) of the GPML book). The abovementioned paper (Eq. (18)) and the source code of both GPML and GPstuff all use the equation: $$\frac{\partial \hat{f}}{\partial \sigma} = (I + K W)^{-1} K \frac{\partial}{\partial \sigma} \nabla \log p(y | \hat{f}) \text.$$ I tried but could not figure out how to derive this same equation.

Any reference or hint to derive this equation is appreciated. Thanks!

P.S.: I originally posted this question on math stackexchange but has since figured out the stats stachexchange is probably more appropriate for this question.