How can one compute the Shapley value using Monte Carlo?

Question

I am working on an application for which I need to compute the Shapley value (using the first formula from the definition section on the wiki). Since I am dealing with a large number of players (around 300), this clearly cannot be computed exactly, hence the need for some Monte Carlo-like algorithm. Here is what I am doing now:

while(epsilon < 5.0):
    for player in playerSet: 
        random_player_coalition = generateCoalitions(playerSet)

        # scaled, including the factor in front from the definition
       player_marginal_contribution = computeMarginalContribution(player, random_player_coalition)

        # list of marginal contributions for every coalition considered
        # of the form [marg1, marg2, marg3..., margN]
        player_shap_list += [player_marginal_contribution]

        prev_shap = sum(player_shap_list[:-1])
        curr_shap = sum(player_shap_list)

        delta_shap = abs(curr_shap - prev_shap)

        # percentage of change between steps
        epsilon = delta_shap * 100.0 / prev_shap

Is this the "proper" way to perform such an approximation? In reality, I have a threshold that says if epsilon is below a particular percentage for a number of steps, then consider convergence achieved.

My questions are two-fold:

1. Is this algorithm correct, in my assumption that it will reach some kind of convergence?
2. While running it, I noticed that it seems to be converging, as epsilon keeps decreasing with every step, however it's very very slow and takes a large number of samples. What are some tricks I can do to speed it up?

I am new to this area of Monte Carlo simulations, so I can use all the help I can get.