I know variants of this question have been asked a million times, but rather than just asking "how do I derive the covariance matrix" I ask you to check the error in my calculations, because I am stumped...
I have the least squares problem
$$\min_p \frac{1}{2}\lVert r(p)\rVert^2,$$
where $r(p) = W(y-f(p))$ is my residual vector. Here $y\in \mathbb{R}^{N_y}$ is the vector of observations, $f(p)$ is the vector-valued model function and $p \in \mathbb{R}^{N_p}$ are the parameters. The weight matrix is $W=\text{diag}(1/\sigma_1,...,1/\sigma_{N_y})$, with $\sigma_j$ the standard deviation of data at index $j$.
What I am trying to do is to derive this expression for the covariance matrix of the best-fit parameters $p^\dagger$:
$$C_p = \sigma^2 (J^T J)^{-1}$$
where $J$ is the Jacobian of $r$ and $\sigma^2$ is
$$\sigma^2 = \frac{\lVert r(p^\dagger)\rVert^2}{N_y-N_p}$$
What I tried
I tried taking a Bayesian approach, where least squares minimization is derived as the max-a posteriori estimate of the parameters given the observations under uniform priors. I'll gloss over this part a bit, because this is all well known. I'm happy to provide clarification. Also I'm happy to hear if I glossed over mistakes here.
Assuming statistical independence of all $y_i$, $y_j$ I should be able to write the posterior as:
$$P(p|y) = K \exp\left(-\frac{1}{2} \sum_j \frac{(y_j-f_j(p))^2}{\sigma_j^2} \right) = K \exp\left( -\frac{1}{2} \lVert r(p) \rVert^2 \right) = K \exp(-g(p)),$$
where $K$ is a constant of integration and for simplicity I wrote:
$$g(p) := \frac{1}{2} \lVert r(p) \rVert^2$$
Then I use I Taylor expansion for $g(p)$ around the best-fit parameter $p^\dagger$, which is
$$g(p) \approx g(p^\dagger) + \nabla g(p^\dagger)(p-p^\dagger) + (p-p^\dagger)^T H(p^\dagger) (p-p^\dagger) = g(p^\dagger) + \frac{1}{2}(p-p^\dagger)^T H(p^\dagger) (p-p^\dagger),$$
where $H(p^\dagger)$ is the Hessian of $r(p)$ evaluated at $p^\dagger$ and I have used the fact that the gradient vanishes at $p^\dagger$ (see here p.3 for the formulae).
Now I can approximate the posterior as
$$P(p|y) = K^\prime \exp\left( -\frac{1}{2}(p-p^\dagger)^T H(p^\dagger) (p-p^\dagger) \right),$$
where I have absorbed the further constant factors into $K^\prime$. To my mind this is a multivariate normal distribution with covariance matrix $H(p^\dagger)^{-1}$. Now I know that we often approximate the Hessian in least squares approximation as $H=J^TJ$, so that I get for the covariance matrix:
$$C_p = (J^T J)^{-1}$$
Compared to the expression that I wanted to derive, this is missing the $\sigma^2$ and I have no idea where my logic is flawed. I appreciate all help.
UPDATE: the Univariate Case
@whuber suggested to work through the univariate case to spot errors. So for maximum simplicity I have only one parameter $p\in \mathbb{R}$ and one data point $y\in\mathbb{R}$.
For the univariate case we can write our likelihood as
$$P(y|p) = \frac{1}{\sqrt{2\pi\sigma^2}} \exp \left( -\frac{1}{2} \frac{(y-f(p))^2}{\sigma^2}\right)$$
Now by the same arguments as above the posterior is
$$P(p|y) = K \exp(g(p))$$
where I again define the "sum of squares" as $g(x)$, where
$$g(p) = \frac{1}{2} \frac{(y-f(p))^2}{\sigma^2}$$
Now we Taylor expand $g$ again around the best fit, such that
$$g(p) \approx g(p^\dagger) + g'(p^\dagger)(p-p^\dagger) + 1/2 g''(p^\dagger) (p-p^\dagger)^2$$
where $g'$ and $g''$ are the first and second derivative of $g$ with respect to $p$. Again, the first derivative vanishes and I end up with
$$P(p|y) \approx K' \exp(-\frac{1}{2}g''(p^\dagger) (p-p^\dagger))$$
which again is a Gaussian with variance $\sigma^2=(g''(p^\dagger))^{-1}$, which is the same result as above... I can even approximate $g'(p) \approx r'(p)\cdot r'(p)$... Maybe this example is too simplified because I have only one data point and one degree of freedom??