I recently wrote some code that I thought very inefficient, but since it included only a few values, I accepted it.  However, I'm still interested at a better algorithm for the following:A list of X objects, each of them are assigned a weightSum up the weightsGenerate a random number from 0 to the sumIterate through the objects, subtracting their weight from the sum until the sum is non-positiveRemove the object from the list, and then add it to the end of the new listItems 2,4, and 5 all take n time, and so it is an O(n^2) algorithm.Can this be improved?As an example of a weighted shuffle, an element has a greater chance at being at the front with a higher weight.Example (I'll generate random numbers to make it real):6 objects with weights 6,5,4,3,2,1; Sum is 21I picked 19: 19-6-5-4-3-2 = -1, thus 2 goes in the first position, weights are now 6,5,4,3,1; Sum is 19I picked 16: 16-6-5-4-3 = -2, thus 3 goes in the second position, weights are now 6,5,4,1; Sum is 16I picked 3: 3-6 = -3, thus 6 goes in the third position, weights are now 5,4,1; Sum is 10I picked 8: 8-5-4 = -1, thus 4 goes in the fourth position, weights are now 5,1; Sum is 6I picked 5: 5-5=0, thus 5 goes in the fifth position, weights are now 1; Sum is 1I picked 1: 1-1=0, thus 1 goes in the last position, I have no more weights, I finish
How to implement a Weighted shuffle
algorithms
Let me add this solution to answer @muskrat better that I could do with yet another comment.I said in my comment that the answer of Winston Ewert should be changed to use exponential distributions to work as expected.More precisely, if ki are independent random samples drawn from an exponential distribution with a rate parameter lamba = 1/wi (assuming wi>0), then keeping the n greatest samples provide a solution to the weighted urn problem with n draws and weights wi.So why does this work? The joint probability of (ki, kj) isp(x,y) = 1/(wi*wj) * exp -(x/wi + y/wj)Now the probability that ki < kj is (after some elementary but boring calculus)p(ki < kj) = ∫(x<y) p(x,y) dx dy = wj/(wi+wj)so that the ratio p(ki < kj) / p(kj < ki) = wj/wiis proportionnal to the ratio of the weights, as desired.EDIT 2As mentionned by @Gordak, the property above does not translate to the global property thatp(ki > kj, ∀j) = wi/∑jwjand therefore may not be appropriate for random shuffling. I modified the code below to display this discrepency. Since the last two variables to be picked would need to have their variables obeying this rule, I am wondering if this means that random shuffling is actually not possible by this technique of sorting random variables.EDITAs requested by @MrTi, here is a small program that not only does the weighted shuffle but actually empirically tests that the resulting order follows the requested weights. (Without loss of generality, it only tests for the relative order of the first two variables).sample_num = 10;% weights are chosen randomly here for test purposes    weights = rand(1, sample_num);round_num = 10000;count = 0;countfirst = zeros(sample_num, 1);for round_idx = 1:round_num  % draw samples from exponential distributions with given weights  scores = arrayfun(@(x) exprnd(x), weights);  [~, I] = sort(scores);  % count the number of time first variable comes ahead of the second  if find(I==1) > find(I==2), count = count + 1; end  countfirst(I(end)) = countfirst(I(end)) + 1;endfprintf('Frequency of the first variable appearing before the second:\n')theoretical = weights(1)/ (weights(1) + weights(2));empirical = count / round_num;fprintf('theoretical: %f\n', theoretical);fprintf('empirical: %f\n', empirical);fprintf('Frequency of variables appearing first:\n')theoretical = weights / sum(weights);empirical = countfirst / round_num;fprintf('theoretical ratio: '); fprintf('%f ', theoretical); fprintf('\n');fprintf('empirical ratio: '); fprintf('%f ', empirical); fprintf('\n');This Matlab/Octave code can be tested online e.g. here.