Working out error on fit parameters for nonlinear fit I am struggling to find a concrete formula for the Hessian or Jacobian in respects to fitting parameters. 
I have implemented some fitting in Java using the Apache Common Maths package for the function $ y = a * tan^{-1}(b x)$ where $a$ and $b$ are fitting parameters. The fitting package does not seem to provide the covariance matrix. 
The derivatives w.r.t to the parameters are straight-forwardly
$$ \frac{\partial y}{\partial a} = \frac{\partial y}{\partial x_1} = tan^{-1}(bx) \ , \ \ \frac{\partial y}{\partial b} = \frac{\partial y}{\partial x_2} =\frac{ax}{b^2 x^2 + 1} $$
What is the correct of calculating the Jacobina or Hessian or Covariance matrix?
The generic formula for the Jacobian is
$$J_{ij} = \frac{\partial f_i}{\partial x_j}.$$
In the case of fitting, $\vec{x}$ are the parameters we are fitting, i.e. $a$ and $b$. What I cannot work out is what $f_i$ should be.
Once I have the Jacobian, I can estimate the Hessian as $H = J^{T} J$ and then invert it to get the covariance matrix, 
$$cov = H^{-1}$$
scale it 
$$cov_s = cov * \frac{RSS}{dof}$$
where $RSS$ is the residual of the sum of squares and $dof$ is the degrees of freedoms. 
The uncertainty of the $i$th parameter then is 
$$\sqrt{cov_{s, ii}}$$
How to calculate the Jacobian? Or is there a better way of estimating the uncertainty of the fitting parameters?
 A: The $f_i$ is the residual of the model fit for the $i$th data point.
$f_i = y_i - a * \tan^{-1}(b x_i)$ 
$J^TJ$ may be inaccurate estimate of Hessian:
Note that the Hessian equals $J^TJ$ plus a (higher order) matrix. This additional matrix is close to the zero matrix if  the residuals are small enough.  But if the optimal residuals are not small, $J^TJ$ can severely underestimate the true Hessian, resulting in underestimation of the true uncertainty of the parameters.  You can eliminate this source of inaccuracy by calculating the actual Hessian, rather than using $J^TJ$.
Bootstrap provides most accurate estimate of uncertainty of parameters: Furthermore, the inherent "linearization" involved in using (RSS/(degrees of freedom) * inverse of Hessian as estimate of the covariance matrix limits its accuracy. Presuming there is no systematic structure in your data, you should be able to get the most accurate estimate of the covariance, and more generally the uncertainty in the parameter estimates, by using the bootstrap.  Generate some number $M$, say 100 or 1000, of bootstrap samples of size $n$, with replacement, where $n$ is the number of data points you have. For each bootstrap sample, calculate the least squares estimate of the parameters.  At the end of this procedure, you will have $M$ samples of parameter pairs ($\hat{a},\hat{b})$.  This collection of $M$ samples represents your best estimate of the distribution of $(a,b)$.  You can calculate the sample covariance of these $M$ samples and/or make assessments directly from these $M$ pairs.
Edit: Note that to make $J^TJ$ match up as an approximation to the Hessian, and for use of scaled inverse Hessian for covariance, the function of which the Hessian is approximated or computed is $\frac{1}{2}\Sigma_{i=1}^nf_i^2$
