Составные операторы присваивания в библиотеке Numpy Python

«Векторизация» причудливой индексации библиотекой numpy Python иногда дает неожиданные результаты. Например:

import numpy
a = numpy.zeros((1000,4), dtype='uint32')
b = numpy.zeros((1000,4), dtype='uint32')
i = numpy.random.random_integers(0,999,1000)
j = numpy.random.random_integers(0,3,1000)

a[i,j] += 1
for k in xrange(1000):
    b[i[k],j[k]] += 1

Дает разные результаты в массивах 'a' и 'b' (т.е. появление кортежа (i,j) отображается как 1 в 'a' независимо от повторов, тогда как повторы подсчитываются в 'b'). Это легко проверить следующим образом:

numpy.sum(a)
883
numpy.sum(b)
1000

Примечательно также, что причудливая версия индексации почти на два порядка быстрее, чем цикл for. Мой вопрос: «Существует ли эффективный способ для numpy вычислять счетчики повторов, реализованные с помощью цикла for в приведенном примере?»