Directly discarding big term in the proof of error propagation formula of variance from random variable $x$ to $f(x)$? I read the error propagation formula scanario said that, the connection between the variance of a random variable $x$ and $f(x)$ is $\frac{var(f(x))}{|\partial_xf(x)|_{x=\bar{x}}|^2}=var(x)$. While I'm confused in the prove of it in the process of Taylor series:
The proof of it states that, if we define the expected value of random variable $x$ and variance of it as follows:
$$
\begin{aligned}
\bar{x} &=\sum_{i=1}^{N} p_{i} x_{i} \\
(\Delta x)^{2} &=\sum_{i=1}^{N} p_{i}\left(x_{i}-\bar{x}\right)^{2} \\
&=\overline{x^{2}}-\bar{x}^{2}
\end{aligned}
$$
And the same for $f(x)$:
$$
\begin{aligned}
\overline{f(x)} &=\sum_{i=1}^{N} p_{i} f\left(x_{i}\right) \\
(\Delta f)^{2} &=\sum_{i=1}^{N} p_{i}\left(f\left(x_{i}\right)-\bar{f}\right)^{2}
\end{aligned}\tag{1}
$$
Then we can Taylor expanse it as:
$$
\begin{aligned}
f\left(x_{i}\right) &=\left.\sum_{n !}^{\infty} \frac{1}{n !}\left(x_{i}-\bar{x}\right)^{n} f^{(n)}(x)\right|_{x=\bar{x}} \\
&\approx f(\bar{x})+\left(x_{i}-\bar{x}\right) f^{(1)}(\bar{x})
\end{aligned}\tag{2}
$$
ignoring the square and higher-order term. Then $\overline{f(x)}$ can be written as
$$
\begin{aligned}
\overline{f(x)} &\approx f(\bar{x})+\overline{\left(x_{i}-\bar{x}\right)} f^{(1)}(\bar{x}) \\
&=f(\bar{x})
\end{aligned}\tag{3}
$$
Then replace eq (3) and (2) into (1)， we can get $\frac{var(f(x))}{|\partial_xf(x)|_{x=\bar{x}}|^2}=var(x)$.
My question is, $x_i-\bar{x}$ is not always a small term!!! Why should we just discard it? Or, are there more formal references of this error propagation formula, because all I search online is some other kind of form.

The computation like it is shown above discards terms which have higher order than 1 when the function is expanded. Therefore, the proposed formula only applies to functions f which are ("well enough") representable as a linear function. This is the level of error propagation like you will typically find it.
If you want to do error propagation beyond linear taylor expansion because your function f is highly non-linear (or to check the assumptions made during propagation of error) I suggest to go back to the idea of error propagation itself. The following image shows the limitation of using a linearization (resulting in $+/- \Delta y$) instead of all terms of the taylor expansion (resulting in $+\Delta y_+$ and $- \Delta y_-$):

Biggerj1, CC BY-SA 4.0 https://creativecommons.org/licenses/by-sa/4.0, via Wikimedia Commons
Picture source: https://commons.wikimedia.org/wiki/File:Propagation_of_error.svg
