Создать пуассоновоподобное распределение с N случайными числами, сумма которых является константой (C)

Question

Я хочу сгенерировать случайное распределение чисел Пуассона, где сумма сгенерированных чисел равна 1000, а нижняя верхняя граница распределения равна (3-30).

Я могу использовать numpy для генерации случайного числа:

In [2]: np.random.poisson(5, 150)
array([ 4,  4,  6,  4,  8,  6,  4,  2,  6,  8,  8,  8,  1,  4,  3,  4,  1,
        3,  7,  6,  7,  4,  5,  5,  7,  6,  5,  3,  3,  5,  4,  6,  2,  0,
        3,  5,  6,  2,  5,  2,  4,  7,  4,  7,  8,  5,  6,  1,  4,  4,  7,
        4,  7,  2,  7,  4,  3,  8, 10,  2,  5,  7,  6,  3,  5,  7,  8,  5,
        4,  7,  8,  8,  2,  2, 10,  6,  3,  5,  2,  5,  5,  6,  4,  6,  4,
        0,  4,  3,  5,  8,  6,  7,  4,  4,  4,  3,  3,  4,  4,  6,  7,  6,
        3,  9,  7,  7,  4,  5,  2,  4,  3,  6,  5,  6,  3,  6,  8,  9,  6,
        3,  4,  4,  7,  3,  9, 12,  4,  5,  5,  7,  6,  5,  2, 10,  1,  3,
        4,  4,  6,  5,  4,  4,  7,  5,  6,  5,  7,  2,  5,  5])

Но я хочу добавить к этому что-то еще:

- The random number should be minimal of 3 and max of 30 
- The sum of the generated random number should be 1000.

Я знаю, что, возможно, я не создаю точное распределение Пуассона, если буду манипулировать.Но я хочу что-то вроде Пуассона, но с предлагаемыми элементами управления.

Answer 1

Позвольте мне написать что-то, что могло бы работать или нет, мы увидим

Свойство распределения Пуассона состоит в том, что один параметр - λ - это мера среднего значения и дисперсии одновременно.Давайте попробуем другой дистрибутив, который действительно суммирует до 1000 и достаточно близок к Пуассону.

Я бы попробовал Полиномиальное распределение .Давайте рассмотрим выборку 200 чисел из полинома.Мы будем сдвигать каждое выбранное число на 3, поэтому выполняется минимальное граничное условие.Это означает, что для выборочной полиномиальной суммы (параметр n) равен 1000 - 3 * 200 = 400. Вероятности p _i будут установлены на 1/200.

Таким образом, для полиномиального среднего E [x _i ] = np _i = 400/200 = 2. Отклонение от полинома будет равно = np _i (1 - p _i ), а поскольку p _i очень мало, член (1 - p _i ) будет очень близок к 1,таким образом делая выборочные целые числа, напоминающие Пуассона, со средним, равным дисперсии.Проблема в том, что после сдвига среднее значение будет равно 5, но дисперсия останется на уровне ~ 2.

В любом случае, некоторый код.

import numpy as np

N = 200
shift = 3
n = 1000 - N*shift
p = [1.0 / float(N)] * N

q = np.random.multinomial(n, p, size=1)
print(np.sum(q))
print(np.mean(q))
print(np.var(q))

result = q + shift
print(np.sum(result))
print(np.mean(result))
print(np.var(result))

Answer 2

Вот еще одна альтернатива, основанная на предварительном распределении минимума для каждой ячейки, вычислении количества оставшихся наблюдений и наборе скорости Пуассона для каждой оставшейся ячейки, определяемой количеством наблюдений и количеством оставшихся ячейок, но подлежащих принятию /отклонение на основе верхней границы для бина.

Поскольку пуассон представляет собой подсчет количества наблюдений, попавших в интервал, если не все были выделены на начальном этапе, они случайным образом распределяются один за другимЕмкости с оставшейся емкостью.

Вот оно:

import numpy as np

def make_poissonish(n, num_bins):
    if n > 30 * num_bins:
        print("requested n exceeds 30 / bin")
        exit(-1)
    if n < 3 * num_bins:
        print("requested n cannot fill 3 / bin")
        exit(-1)

    # Disperse minimum quantity per bin in all bins, then determine remainder
    lst = [3 for _ in range(num_bins)]
    number_remaining = n - num_bins * 3

    # Allocate counts to all bins using a truncated Poisson
    for i in range(num_bins):
        # dial the rate up or down depending on whether we're falling
        # behind or getting ahead in allocating observations to bins
        rate = number_remaining / float(num_bins - i)  # avg per remaining bin

        # keep generating until we meet the constraint requirement (acceptance/rejection)
        while True:
            x = np.random.poisson(rate)
            if x <= 27 and x <= number_remaining: break
        # Found an acceptable count, put it in this bin and move on
        lst[i] += x
        number_remaining -= x

    # If there are still observations remaining, disperse them
    # randomly across bins that have remaining capacity
    while number_remaining > 0:
        i = np.random.randint(0, num_bins)
        if lst[i] >= 30:    # not this one, it's already full!
            continue
        lst[i] += 1
        number_remaining -= 1
    return lst

Пример вывода:

result = make_poissonish(150, 10)
print(result)                    # => [16, 19, 11, 16, 21, 18, 12, 17, 8, 12]
print(sum(result))               # => 150

result = make_poissonish(50, 10)
print(result)                    # => [3, 5, 5, 4, 3, 3, 15, 3, 6, 3]
print(sum(result))               # => 50

Answer 3

Вы можете легко сделать это, используя цикл while и случайный модуль, и он сделает эту работу:

from random import randint
nums_sum = 0
nums_lst = list()
while nums_sum < 1000:
    n = randint(3, 31)
    nums_sum += n
    nums_lst.append(str(n))
    print(nums_sum)
    if 1000-nums_sum > 30: # means if the sum is more than 30 then complete ..
        continue
    else:
        nums_sum += 1000-nums_sum
print(nums_sum)
print(nums_lst)

так просто.