# Calculating a Weighted Standard Error of the Fit for Nonlinear Regression

I have a data set of $$N$$ points to which I have fit an equation of $$n$$ parameters $$\theta_{1..n}$$ such that $$y_i \sim f(x_i; \theta_{1..n})$$. These data $$(x_{1..N},y_{1..N})$$ have been provided with standard deviations $$\sigma_{1..N}$$ after having been acquired several times. My nonlinear regression took this into account by minimizing the weighted RSS, where the reciprocal variances $${\sigma_{1..N}}^{-2}$$ were used as weights. I now want to estimate the confidence intervals and covariance for parameters $$\theta_{1..n}$$.

I would have liked to resort to bootstrapping, but the process of resampling is too computationally expensive for the function $$f$$ being used. Instead, I suppose I would like to fall back on estimating the (co)variance matrix using the Jacobian. In the case of an unweighted covariance matrix, I know this can be done (assuming $$J^\mathsf{T}J$$ is invertible) as $$V(\theta_{1..n}) = \sigma^2 \left(J(\theta_{1..n})^\mathsf{T} J(\theta_{1..n})\right)^{-1}$$ for $$n\times n$$ covariance matrix $$V$$, variance $$\sigma^2=\frac{1}{N-n} \sum_{i=1}^N \left(y_i - f(x_i; \theta_{1..n})\right)^2$$, $$N\times n$$ Jacobian matrix $$J$$ (entrywise $$J_{ij}=\frac{\partial}{\partial\theta_j}f(x_i;\theta_{1..n})$$), and $$\theta_{1..n}$$ the optimized parameters. Here, I can estimate the Jacobian with numerical differentiation, allowing evaluation of $$f$$ only a couple of times per parameter $$\theta$$ instead of hundreds/thousands of time as in bootstrapping. How can I adapt this formula to estimate a weighted covariance matrix? Is it as simple as using weighted $$\sigma^2$$ and $$J$$? If so, how does one calculate a weighted'' Jacobian matrix?

Alternatively, is there a more appropriate way to approach such an issue for an expensive function $$f$$? I'd like to make use of the error bars on the data points to their maximum extent to inform the error bars on the fitted parameters.