Derivation of variance The variance of a random variable $X$ is defined as:
$$
\text{var}(X) = \mathbb{E}[(X - \mathbb{E}[X])^2]
$$
And the definition of expectation for a discrete random variable is:
$$
\mathbb{E}[X] = \sum_{i=1}^{n} x_i \cdot \mathbb{P}(\{X=x_i\})
$$
Combining these two equations and letting $\mathbb{E}[X] = \mu$, I would have said that variance for a discrete random variable $X$ is:
$$
\begin{align}
\text{var}(X) &= \mathbb{E}[(X - \mathbb{E}[X])^2] \\
&= \sum_{i=1}^{n} (x_i - \mu)^2 \cdot \mathbb{P}(\{X = (x_i - \mu)^2\}) \\
\end{align}
$$
But Wikipedia states that the correct formulation is:
$$
= \sum_{i=1}^{n} (x_i - \mu)^2 \cdot \mathbb{P}(\{X = x_i)
$$
What am I missing? Assuming my thinking is correct, this would suggest that:
$$
\mathbb{P}(\{X = (x_i - \mu)^2\}) = \mathbb{P}(\{X = x_i)
$$
My inclination is that the above equation is true since $\mu$ and squaring are both non-random.
 A: Technically, $X$ wouldn't take the value of it's squared deviation or residual. For instance, if $X$ is Bernoulli $\mathbb{P}(X = (x_i - \mu)^2)$ is 0. 
If you define a new variable for the squared deviation/residual, $Y$ so that $y_i = (x_i - \mu)^2$, this has a different distribution and $E[Y]$ is calculated as $\sum_{i=1}^n y_i \mathbb{P}(Y = y_i)$ and is equal to $\text{var}(X)$. A long-winded way of writing that expression is similar to the one you give except that we would say $$\text{var}(X) = \sum_{i=1}^n (x_i - \mu)^2 \mathbb{P}((X-\mu)^2 = (x_i - \mu)^2)$$.
A: Maybe first you should tell us why would you expect $\mathbb{P}(\{X = (x_i - \mu)^2\})$ to see there..? You can easily make up examples of bounded random variables where $X$ would never be equal to $(x_i-\mu)^2$, i.e. the probabilities would be equal to zero and the whole equation would blow up...
This can be easily verified numerically, without going into calculus. The simple example is continuous uniform $\mathcal{U}(-1, 3)$ distribution. It has a variance equal to $(b-a)^2/12 = 1.3333$, if we simulate it and estimate the variance as it is defined and using empirical variance, then both estimates are reasonably close to the correct answer. Same if we numerically integrate the function. However if we numerically integrate your function, it returns a wrong answer.
set.seed(123)
x <- runif(1e6, -1, 3)

mean((x-mean(x))^2)
## [1] 1.334592
var(x)
## [1] 1.334593

integrate(function(x) (x-1)^2 * dunif(x, -1, 3), lower = -1, upper = 3)
## 1.333333 with absolute error < 1.5e-14
integrate(function(x) (x-1)^2 * dunif((x-1)^2, -1, 3), lower = -1, upper = 3)
## 0.8660249 with absolute error < 5.5e-05

More formally, let's go back to the definition of variance
$$
\operatorname{Var}(X) = E[(X-E(X))^2] = E[X^2] - E[X]^2 
$$
as you can see, the expanded definition defines it in terms of exceptions over $X$. You can recall that by LOTUS, the expected value of a function of random variable is
$$
E[g(X)] = \int_{-\infty}^\infty g(x) \, f(x)\, dx 
$$
and this is exactly your case. 
