The sampling approach rigorously would look like this. For each player $i$, we want to estimate the expected marginal contribution, where the expectation is taken over a the subset of players that precede $i$ in the permutation ordering. So for each $i$, we do the following.
Let $X_i$ be a random variable equal to the marginal contribution of $i$, when we draw the permutation of players randomly. Let $\mu_i$ be the true Shapley value of $i$ (which we do not know). Then $\mathbb{E} X_i = \mu_i$. Now, if we sample many independent copies of $X_i$ and average them, this average, call it $\bar{X}_i$, should be very close to $\mu_i$, and the closeness is given by e.g. Hoeffding's inequality, which says
Let $K$ be an upper bound on any marginal contribution (so they are always between $0$ and $K$). Let $m$ be the number of independent copies of $X_i$ that we draw. Then
$$ \Pr[ | \bar{X}_i - \mu_i | \geq \epsilon ] \leq 2e^{-2m \epsilon^2 / K^2} . $$
So for each $i$, we sample $m$ permutations randomly and calculate the average marginal contribution $\bar{X}_i$.
For example, if $K=100$ and we want an accuracy of $\epsilon = 0.01$ and a probability of failure of $2e^{-50}$, then we need $2m\epsilon^2/K^2 = 50$, so $m = 25K^2/\epsilon^2 = 2.5$ billion. So we need to sample $2.5$ billion permutations to "guarantee" accuracy of $0.01$, except with a miniscule $2e^{-50}$ chance of failure.
As comments have mentioned, there are big improvements for many special cases.