Aloha :)
Im Folgenden betrachten wir eine Funktion \(f(x_1,x_2,\ldots,x_N)\), die von \(N\) Messwerten \(x_i\) abhängt. Wir tun so als würden wir alle \(N\) Messwerte \(n\)-Mal bestimmen und daraus jeweils einen Wert für \(f\) berechnen. Das heißt, wir haben \(N\) Messwerte und \(n\) Wiederholungen. Bitte das im Folgenden nicht durcheinander werfen.
Aus den \(n\) Wiederholungen erhalten wir die Ergebnisse \(\{f_1,f_2,\ldots,f_n\}\) und können daraus die Varianz berechnen:$$(\delta f)^2=\frac{1}{n-1}\sum\limits_{k=1}^n\left(f_k-\overline f\right)^2$$Gemäß des zentralen Grenzwertsatzes können wir die \(N\) Mittelwerte \(\overline x_i\) aus allen \(n\) Wiederholungen den realen Werten gleichsetzen. Weiter können wir für jede der \(n\) Wiederholungen die Funktion \(f(x_1,x_2,\ldots,x_N)\) durch eine Talyor-Reihe erster Ordnung um die Mittelwerte entwickeln: $$f_k(x_1,x_2,\ldots,x_N)=f_k(\overline x_1, \overline x_2, \ldots, \overline x_N)+\sum\limits_{i=1}^N\frac{\partial f}{\partial x_i}\left(x_{i,k}-\overline x_i\right)$$Hier bitte nicht durcheinander kommen. Der Index \(k\) zählt die \(n\) Wiederholungen, der Index \(i\) zählt die \(N\) Messwerte. Daher ist \(x_{i,k}\) der Messwert \(x_i\) bei der \(k\)-ten Wiederholung.
Das \(f_k(\overline x_1,\overline x_2,\ldots,\overline x_N)\) hat, weil wir für jede Wiederholung \(k\) immer dieselben Mittelwerte einsetzen, einen konstanten Wert, es ist der gesuchte Mittelwert \(\overline f\) der Funktion. Daher ist:$$f_k-\overline f=f_k(x_1,x_2,\ldots,x_N)-f_k(\overline x_1, \overline x_2, \ldots, \overline x_N)=\sum\limits_{i=1}^N\frac{\partial f}{\partial x_i}\left(x_{i,k}-\overline x_i\right)$$Das können wir oben einsetzen:$$(\delta f)^2=\frac{1}{n-1}\sum\limits_{k=1}^n\left(\sum\limits_{i=1}^N\frac{\partial f}{\partial x_i}\left(x_{i,k}-\overline x_i\right)\right)^2$$$$\phantom{(\delta f)^2}=\frac{1}{n-1}\sum\limits_{k=1}^n\left(\sum\limits_{i=1}^N\frac{\partial f}{\partial x_i}\left(x_{i,k}-\overline x_i\right)\sum\limits_{j=1}^N\frac{\partial f}{\partial x_j}\left(x_{j,k}-\overline x_j\right)\right)$$$$\phantom{(\delta f)^2}=\frac{1}{n-1}\sum\limits_{k=1}^n\left(\sum\limits_{i=1}^N\sum\limits_{j=1}^N\frac{\partial f}{\partial x_i}\frac{\partial f}{\partial x_j}\left(x_{i,k}-\overline x_i\right)\left(x_{j,k}-\overline x_j\right)\right)$$$$\phantom{(\delta f)^2}=\sum\limits_{i=1}^N\sum\limits_{j=1}^N\frac{\partial f}{\partial x_i}\frac{\partial f}{\partial x_j}\underbrace{\frac{1}{n-1}\sum\limits_{k=1}^n\left(x_{i,k}-\overline x_i\right)\left(x_{j,k}-\overline x_j\right)}_{=Cov\,(x_i,x_j)}$$$$(\delta f)^2=\sum\limits_{i=1}^N\sum\limits_{j=1}^N\frac{\partial f}{\partial x_i}\frac{\partial f}{\partial x_j}\,\delta(x_i,x_j)$$
Die Formel, die du aus Wikipedia zitiert hast, ist exakt dieselbe. Das siehst du, indem du zunächst aus der Doppelsumme die Terme mit \(i=j\) rausziehst:$$(\delta f)^2=\sum\limits_{i=1}^N\left(\frac{\partial f}{\partial x_i}\,\delta x_i\right)^2+\sum\limits_{i,j=1\;;\;i\ne j}^N\frac{\partial f}{\partial x_i}\frac{\partial f}{\partial x_j}\,\delta(x_i,x_j)$$Wegen der Symmetrie der Summanden in der zweiten Summe, braucht man nur die Terme \(j>i\) zu berechnen und kann danach das Ergebnis verdoppeln:$$(\delta f)^2=\sum\limits_{i=1}^N\left(\frac{\partial f}{\partial x_i}\,\delta x_i\right)^2+2\sum\limits_{i=1}^{N-1}\sum\limits_{j=i+1}^N\frac{\partial f}{\partial x_i}\frac{\partial f}{\partial x_j}\,\delta(x_i,x_j)$$
