To distinguish 1000 particles you need 10 bits. For programming convenience, let's say you dedicated two bytes (16 bits) to identify each of 1000 particles. If the temporary list has 1,000,000 entries, we're talking about a grand total of 2MB. That's 0.002GB. Is 2MB an amount that runs the risk of using up all of memory these days? Lookup would be quite brisk using such a table, being O(1). There'd be one multiplication (that's a clock cycle) and one list entry retrieval (what's that, another clock cycle?). What are you aware of that's either faster or easier to understand? This approach should not be dismissed.

As for losing out on some very low probability particles because we may not have chopped the list finely enough... so what? One entry in a list of 1,000,000 entries has a one-in-a-million chance of being selected. Are we supposed to be worrying about particles that have even less of a chance of being selected, too?

If it aches at you that too much RAM is being used, consider making the list length 100,000 instead. Now you're using 200K (0.0002GB) and the process is even faster to set up (but probably runs at the same super-fast speed as before). Hardly any significance is lost when it comes to representing important particles.

@max, following is one elegant way to implement your method #1. Features: 1. No extra memory needed for the new population or temporary array. 2.Preserve / No change to the value of existing arrays and,
3. Minimum calculations (just 2 loops and no multiplication / division)

p3 = []
sumW = sum(w)              # sum of all weights
for i in range(N):         # sample N times
   val = random.uniform(0,sumW)  # generate a random number in range (0,sumW)
   tmp = 0
   for j in range(N):      # chose the index according to the marks
      tmp += w[j]
      if tmp >= val :      # if found the proper index to p array
         p3.append(p[j])   #    chose it as the sampling outcome, and
         break             #    skip to the next round of sampling

This is my approach, that seems simpler and more intuitive than the resampling wheel method

p3 = []
s = sum(w)
while len(p3) != N:
    r = int(random.random()* N)
    p_i = random.random()
    if (p_i < w[r]/s) :
        p3.append(p[r])

Why 1000 x 1000? I did it this way:

s = sum(w)

for i in range(N):
    w[i] = w[i]/s

tmp = []
for i in range(N):
    num = int(w[i]*N)
    for j in range(num):
        tmp.append(p[i])

for i in range(N):
    index = int(random.random() * len(tmp))
    p3.append(tmp[index])

and tmp has ~1000 because sum(w)=1 so sum of all int(w[i]*N) is nearly N (or 1000)

I did the same, first hand. Please note that you can normalize the weights and populate the bin at the same time:

s = sum(w)
# Build a bin according to the (normalized) weights
b = []
N2 = N
for i in range(N):
    n = w[i] / s * N2
    for j in range(int(n)):
        b.append(i)

# Choose N times from the bin randomly
for i in range(N):
    r=int(random.random()*len(b))
    c=b[r]
    p3.append(p[c])

The time complexity seems to be O(2N). Please note that the bin can be made smaller than the sample(without losing much precision), by example by setting N2 = N/10. In that case the time complexity will still be O(2N) but the time complexity of the append loop will only be O(0.1N).

I think that the method is not bad, considering its simplicity.

so you have 1000 particles, and your temporary list would then contain 1000000 particles. You'll run out of memory very soon.

It's relatively easy to code but uses vastly more memory. Additionally, as you note, you sacrifice a bit of accuracy. Additionally, the runtime for k digit accuracy is Nlog((10^k)N) which is meaningfully worse even for k=2.

Rusian Ciurca,
That's essentially what I did, and since sometimes rounding can give less than N, I doubled it. I don't see any need to make it 100 or 1000 times larger, as max did. So memory-wise, it doesn't seem bad. And I don't think your approach loses any accuracy. I'm curious what flaws, if any, there are in this approach.

How about this implementation?

sumw=sum(w)
maxw=max(w)
for i in range(N):
    take=0
    while not take:
        index=int(random.random()*N)
        if w[index]/sumw>random.random()*maxw*N/2.0:
            take=1
            p3.append(p[index])

My naive approach. Can someone tell me if it is correct?

normalizer = sum(w)
w = [wt/normalizer for wt in w]
s = 0
buckets = []
for wt in w:
    s+=wt
    buckets.append(s)
random.seed()
N = len(p)

def inbetween(num, s, e, li):
    """binary search for the particle index"""
    pos = s+(e-s)/2
    if num<=li[pos] and num>li[pos-1]:
        return pos
    if num>li[pos] and num<=li[pos+1]:
        return pos+1
    (s,e) = (s,pos) if li[pos]>num else (pos,e)
    return inbetween(num,s,e,li)

random.seed()

p3 = [p[inbetween(random.random(), 0, N, buckets)] for i in range(N)]