Why do you think there is a difference between those two statements of the expectation value? Notice that the cost function only contains the operators $Z$ and $I$, so we only need to measure in the $z$-basis. The classical cost function is
$$
\begin{align}
C &= \sum_{\langle jk \rangle \in \text{edges}} C_{\langle jk \rangle}\\
C_{\langle jk \rangle} &= \frac{1}{2}(1 - g(j)g(k))
\end{align}
$$
where $g(\cdot)$ assigns $+1$ or $-1$ to a vertex, depending on whether it is assigned to the left or the right of the cut. The cost is maximized when for as many edges as possible, we have assigned the vertices of the edge to different sides of the cut. This is directly mapped to a quantum Hamiltonian when we replace $g(\cdot)$ with $Z$ and interpret computational bitstrings as cut assignments.
$$
\begin{align}
C_{\langle jk \rangle} \rightarrow \hat{C}_{\langle jk \rangle} = \frac{1}{2}(I - Z_j \otimes Z_k)
\end{align}
$$
Notice that the Qiskit notebook says that "the problem Hamiltonian specific to the Max-Cut problem up to a constant is..." because we don't care about the identity term in each $\hat{C}_{\langle jk \rangle}$. When we sample the state $|\vec{\gamma},\vec{\beta}\rangle$ and obtain a bitstring, evaluating the Hamiltonian for that string, by construction, yields the classical cost of that cut assignment, so it's not an abuse of notation to evaluate the classical cost function on a classical bitstring. The probabilistic nature of the sum is implicitly included by the fact that the $i$-th bitstring is obtained with some probability that depends on the particular distribution of $|\vec{\gamma},\vec{\beta}\rangle$. The final equality you mention holds in the limit of infinite samples, where for a finite number of samples we would expect an error that scales like $1/\sqrt{N}$$\sigma/\sqrt{N}$ (where $\sigma$ is the true standard deviation).
Please let me know if this helps or if I should include more details somewhere :)