Edge probability and expected number of edges in the configuration model This question is related to question: Probability that exists at least an edge in the configuration model
There is something I do not understand about the computation of the expected number of edges between $i$ and $j$ nodes in the configurational model, $p_{ij}$. The argument given everywhere I've seen is:


*

*There are $2m$ stubs in the network, with $k_i$ in node $i$ and $k_j$ in node $j$.

*Taking one stub from node $i$, there are $k_j$ possible stubs to connect it to node $j$, so the probability to connect it to node $j$ is $\frac{k_j}{2m-1}$, the $2m-1$ because you can not connect it to the same stub you are coming from.

*There are $k_i$ stubs in node i, so the expected number of edges is just adding up the different probabilities and $p_{ij} = k_i \times \frac{k_j}{2m-1}$.


I do not understand step 3. I would think once there has been an edge between nodes $i$ and $j$, the probability to connect the next stub should change accordingly because there is one less available stub at node $j$: $\frac{k_j-1}{2m-3}$. But also, each new stub considered in node $i$ has two less possible stubs to be connected (because every other edge already connected has two stub ends), so the total available edges in the denominator should decrease as well: $2m-3$, $2m-5$, ..., $2m-2k_i-1$.
Instead, I'd proceed in this way:
$$p_{ij} = 1 - \bar{p}_{ij}, $$
where $\bar{p}_{ij}$ is the probability there isn't any edge between nodes $i$ and $j$. Then,
$$\bar{p}_{ij} = \bar{p}_{{i_1}j} \times \bar{p}_{{i_2}j}\times \dots  \times \bar{p}_{{i_{k_i}}j}, $$
where $\bar{p}_{{i_1}j}$ is the probability there isn't an edge between the first stub in node $i$ to node $j$ and $\bar{p}_{{i_1}j} = \frac{2m-1-k_j}{2m-1}$. Analogously for the other stubs, we get
$$\bar{p}_{ij} = \frac{2m-1-k_j}{2m-1}  \frac{2m-3-k_j}{2m-3} \dots \frac{2m-2k_i-1-k_j}{2m-2k_i-1} = \left( 1 - \frac{k_j}{2m-1} \right)   \left( 1 - \frac{k_j}{2m-3} \right)  \dots \left( 1 - \frac{k_j}{2m-2k_i-1} \right).  $$
So
$$p_{ij} = 1- \left( 1 - \frac{k_j}{2m-1} \right)   \left( 1 - \frac{k_j}{2m-3} \right)  ... \left( 1 - \frac{k_j}{2m-2k_i-1} \right).$$
I can recover from this expression the other one in the large number of edges limit $m \to \infty$, then $2m-2k_i-1 \simeq ... \simeq 2m - 3 \simeq 2m - 1$ and
$$p_{ij} \simeq 1- \left( 1 - \frac{k_j}{2m-1} \right)^{k_i} \simeq 1 - \left( 1 - \frac{k_i k_j}{2m-1} \right) = \frac{k_i k_j}{2m-1},$$
where in the second step I have used the series expansion $(1 - x)^a = 1 - ax + \mathcal{O}(x^2)$ for $x \to 0$.
Question: Does this mean that only the expected number of edges between $i$ and $j$ nodes in the configurational model is $p_{ij} = \frac{k_i k_j}{2m-1}$ in the large number of edges $m$ limit? If that is the case, I find it strange because they don't specify it in any of the sources I've looked. Instead, they seem to say $p_{ij} = \frac{k_i k_j}{2m-1}$ is the general expression which in the large number of edges limit becomes $p_{ij} = \frac{k_i k_j}{2m}$.
 A: The difference between your calculation and the standard one
Actually, $\frac{k_i k_j}{2m-1}$ is the exact expected number of edges between nodes $i$ and $j$.
When you compute
$$
   1 - \bar{p}_{{i_1}j} \times \bar{p}_{{i_2}j}\times \dots  \times \bar{p}_{{i_{k_i}}j}
$$
you are computing something different: the probability that there is at least one edge between $i$ and $j$. (That's because the product $\bar{p}_{{i_1}j} \times \bar{p}_{{i_2}j}\times \dots  \times \bar{p}_{{i_{k_i}}j}$ gives the probability that there are no edges.)
However, in the configuration model, it's possible that there are multiple parallel edges between nodes $i$ and $j$. So the expected number of edges will be bigger than the probability that there is at least one edge.
With typical values (but not all values) of $k_i$, $k_j$, and $m$, it's very unlikely that there are multiple edges between $i$ and $j$: much less likely than having one edge. In that setting, the two values are very close, which is what you're seeing.
The expected value calculation, spelled out
Here's a more detailed justification for the expected value calculation. Number the stubs at node $i$ from $1$ to $k_i$, and number the stubs at node $j$ from $1$ to $k_j$. For $1 \le a \le k_i$ and $1 \le b \le k_j$, define the random variable $X_{i,a}^{j,b}$ to be $1$ if we join the $a^{\text{th}}$ stub at $i$ to the $b^{\text{th}}$ stub at $j$. Let $X_i^j$ be the number of edges between $i$ and $j$. Then
$$
    X_i^j = \sum_{a=1}^{k_i} \sum_{b=1}^{k_j} X_{i,a}^{j,b}
$$
and therefore
$$
   \mathbb E[X_i^j] = \sum_{a=1}^{k_i} \sum_{b=1}^{k_j} \mathbb E[X_{i,a}^{j,b}].
$$
Here we use linearity of expectation, which doesn't care that the random variables $X_{i,a}^{j,b}$ are dependent.
Finally, we have $\mathbb E[X_{i,a}^{j,b}] = \frac1{2m-1}$ for any $a$ and $b$. This doesn't care about what any of the other stubs are doing, because this is a calculation for only one pair of stubs. Therefore $\mathbb E[X_i^j] = \frac{k_i k_j}{2m-1}$ because we add up $k_i k_j$ equal terms.
How to think about these expected values
Regarding the calculation of $\mathbb E[X_{i,a}^{j,b}] = \frac1{2m-1}$: here is how to think about this, and related calculations, painlessly. 
We have a randomizing algorithm for generating a graph from the configuration model: 


*

*Pick one of the $2m$ stubs. Choose another one of the $2m-1$ stubs uniformly at random, and connect them.

*Pick one of the $2m-2$ remaining disconnected stubs. Choose one of the $2m-3$ other stubs uniformly at random, and connect them.

*Repeat until all stubs are connected. Then do the configuration-model-to-graph operation which is irrelevant for now.


This is actually a family of algorithms. In the $i^{\text{th}}$ step, we pick one of the $2m-2i$ remaining stubs, in a way I haven't specified, and then pick one of the $2m-2i-1$ other remaining stubs uniformly at random. We can pick the first stub in many ways: at random, or going in a fixed order, or whatever.
The key fact you should convince yourself of is that no matter how we do that, we end up getting one of the $(2m-1)(2m-3)\dotsm (5)(3)(1)$ matchings of the $2m$ stubs uniformly at random. That means that the way we pick one of $2m-2i$ stubs in the $i^{\text{th}}$ step doesn't matter, and we can do whichever thing is most convenient for us. 
When computing $\mathbb E[X_{i,a}^{j,b}]$, the most convenient rule to use is "In the first step, pick the $a^{\text{th}}$ stub out of node $i$ to connect to a uniformly random stub. In the other steps, do whatever." With this rule, it's clear that $\mathbb E[X_{i,a}^{j,b}] = \frac1{2m-1}$. 
The rule we use shouldn't change the calculation of $\mathbb E[X_{i,a}^{j,b}]$. Therefore it's fine that we use a different rule for every $a$ and for every $b$. If we had to use the same rule for every $a$ and $b$, we'd still get $\frac1{2m-1}$ for all of them, but the calculation would be more painful.
