Relation between Covariance matrix and Jacobian in Nonlinear Least Squares

I saw that CovB = inv(J'*J)*MSE in a MATLAB documentation here at http://www.mathworks.com/help/stats/nlinfit.html However, I cant find any sources for the relationship . I believe in it but I need to find a source to refer it. I have looked but I cant find. It would be great if anyone can help me a bit.

Consider the nonlinear least squares problem: minimize $1/2r(x)^Tr(x)$. Let $J$ = Jacobian of r(x). The Hessian of the objective = $J^TJ +$ higher order terms. The Gauss-Newton or Levenberg-Marquardt approximation is to ignore the higher order terms, and approximate the Hessian as $J^TJ$. This approximation for the Hessian is what is used in the formula CovB = inv(J'*J)*MSE in MATLAB's nlinfit.