How can one replace the less memory consuming way of choosing a random value, or is there an alternative to random.choises()?

Question

Python

There is a method:

def get_job_dist():
    all_dists = [*range(20,20000000)]
    weights = []
    for dist in all_dists:
        weights.append(1/dist)
    dist = choices(all_dists,weights,k=1)
    dist = dist[0]
    return dist

Где приходится объявлять огромный массив с числами от 20 до 20000000 и создавать отдельный массив весов, который вмещает в себе столько же элементов. Веса особые - чем больше число из первого массива, тем меньше его вес. По сути - это гиперболическое распределение случайной величины, но вот не знаю, как это сделать нормально.

Answer 1

Вы знаете, я много чего перепробовал и таки нашёл ответ, который, хоть и математически не совсем соответствует моим запросам, но тем не менее, поставленную задачу решает. Пусть кому-то это покажется костыльным велосипедом, но это работает быстрее и, конечно, не напрягает мою память. Это решение не прямолинейное. А математику решения можно и подкорректировать и теперь, если у меня намертво зависал ноутбук, то сейчас я вызываю функцию в цикле с 1000 итерациями и код у меня выполняется примерно за секунду.
Вот сам метод:

def get_job_dist():
    numbers = '0,1,2,3,4,5,6,7,8,9'
    weights = np.arange(1, 11, 1.0)
    weights = [1/x for x in weights]
    numbers = numbers.split(',')
    steps = np.arange(2,6)
    steps_w = [1/x for x in steps]
    steps_w /= sum(steps_w)
    weights /= sum(weights)
    steps_count = np.random.choice(steps, 1, p= steps_w)
    steps_count = steps_count[0]
    dist_str_arr = np.random.choice(numbers, steps_count, p = weights)
    dist_str = ''
    step = 0
    for sym in dist_str_arr:
        if int(sym) < 2 and step < 2 and steps_count < 3:
            dist_str += str(randint(2,9))
        else:
            dist_str += sym
        step += 1
    dist = int(dist_str)
    if dist < 20:
        dist += randint(20, 99)
    return dist

Answer 2

Как минимум можно передавать аргументы в choices не в виде списков, а в виде генераторов.

r=range(20, 2*10**7);
print(choices(r, (1/x for x in r)))

from guppy import hpy; from random import choices; h=hpy(); print(h.heap()); r=range(20, 2*10**7); print(choices(r, (1/x for x in r))); print(h.heap())                                             
Partition of a set of 257681 objects. Total size = 30586285 bytes.
 Index  Count   %     Size   % Cumulative  % Kind (class / dict of class)
     0  77410  30  8264794  27   8264794  27 str
     1  69723  27  5198776  17  13463570  44 tuple
     2  28709  11  2133716   7  15597286  51 bytes
     3  14447   6  2089272   7  17686558  58 types.CodeType
     4  15188   6  2065568   7  19752126  65 function
     5   2013   1  1932488   6  21684614  71 type
     6   3316   1  1515992   5  23200606  76 dict (no owner)
     7    684   0   937144   3  24137750  79 dict of module
     8   2013   1   932560   3  25070310  82 dict of type
     9   1866   1   633536   2  25703846  84 set
<863 more rows. Type e.g. '_.more' to view.>
[9585326]
Partition of a set of 257636 objects. Total size = 30579517 bytes.
 Index  Count   %     Size   % Cumulative  % Kind (class / dict of class)
     0  77406  30  8264562  27   8264562  27 str
     1  69724  27  5198832  17  13463394  44 tuple
     2  28709  11  2133716   7  15597110  51 bytes
     3  14447   6  2089272   7  17686382  58 types.CodeType
     4  15186   6  2065296   7  19751678  65 function
     5   2013   1  1932488   6  21684166  71 type
     6   3310   1  1514552   5  23198718  76 dict (no owner)
     7    684   0   937144   3  24135862  79 dict of module
     8   2013   1   932560   3  25068422  82 dict of type
     9   1866   1   633536   2  25701958  84 set
<859 more rows. Type e.g. '_.more' to view.>

Но что-то мне подсказывает, что это можно математически привести к операции, работающей за O(1) по времени.
UPD
Залез-таки в исходники и считаю важным отметить, что моё решение не полностью решает вашу проблему с памятью.
Вот что (примерно) делается под капотом:

def choices(self, population, weights=None, *, cum_weights=None, k=1):
        """Return a k sized list of population elements chosen with replacement.

        If the relative weights or cumulative weights are not specified,
        the selections are made with equal probability.

        """
        random = self.random
        if cum_weights is None:
            if weights is None:
                _int = int
                total = len(population)
                return [population[_int(random() * total)] for i in range(k)]
            cum_weights = list(_itertools.accumulate(weights))
        elif weights is not None:
            raise TypeError('Cannot specify both weights and cumulative weights')
        if len(cum_weights) != len(population):
            raise ValueError('The number of weights does not match the population')
        bisect = _bisect.bisect
        total = cum_weights[-1]
        hi = len(cum_weights) - 1
        return [population[bisect(cum_weights, random() * total, 0, hi)]
                for i in range(k)]

Для нашего случая функцию можно упростить до:

def choices(a=20, b=2*10**7, k=1):
    from random import Random
    import itertools
    from bisect import bisect
    random = Random().random

    n = b - a
    population = range(a, b)
    weights = (1/x for x in population)
    cum_weights = list(itertools.accumulate(weights))
    # get_weight = lambda idx: 1 / (a + idx)
    # get_cum_weight = lambda idx: ????

    total = cum_weights[-1]
    hi = len(cum_weights) - 1
    return [
        population[bisect(cum_weights, random() * total, 0, hi)]
        for i in range(k)
    ]

И здесь видно, что происходит материализация генератора в список:
cum_weights = list(itertools.accumulate(weights))
Это уже не три огромных списка, а всего лишь один, но эту проблему тоже можно решить.
Короче. Весь вопрос в том, как вычислить за O(1) следующую сумму ряда:
С[n]=1/(a+0)+1/(a+1)+1/(a+2)+...+1/(a+n)
Если это научиться вычислять за O(1), то bisect(cum_weights, random() * total, 0, hi) будет вычисляться за O(log2(n)). И такая же сложность будет у всей функции в целом.
Однако это не предел. Думаю можно сразу вычислять индекс искомого числа за O(1) без двоичного поиска. Но это уже вы задайте как вопрос с тегом "математика".

Answer 3

Используйте numpy/scipy, например scipy.stats.rv_discrete