Генерация N равномерных случайных чисел, сумма которых равна M

Question

Этот вопрос задавался ранее, но я никогда не видел хорошего ответа.

1. Я хочу сгенерировать 8 случайных чисел, сумма которых равна 0,5.

2. Я хочу, чтобы каждое число выбиралось случайным образом из равномерного распределения (т. Е. Простая функция, представленная ниже, не будет работать, потому что числа не будут распределены равномерно).

   def rand_constrained(n,tot):
       r = [random.random() for i in range(n)]  
       s = sum(r)
       r = [(i/s*tot) for i in r] 
       return r
   

Код должен быть обобщаемым, чтобы вы могли генерировать N равномерных случайных чисел, которые суммируются с M (где M - положительное число с плавающей запятой). Если возможно, не могли бы вы также объяснить (или показать с помощью графика), почему ваше решение генерирует случайные числа равномерно в соответствующем диапазоне?

Смежные вопросы, которые не попадают в цель:

Генерация нескольких случайных чисел, равных значению в python (текущий принятый ответ не является равномерным - другой ответ, который является равномерным, работает только с целыми числами)

Получение N случайных чисел с суммой M (тот же вопрос в Java, в настоящее время принятый ответ просто неверен, также нет ответов с равномерным распределением)

Генерация N случайных целых чисел, которые суммируются с M в R (тот же вопрос, но в R с нормальным, а не равномерным распределением)

Любая помощь с благодарностью.

Answer 1

То, что вы просите, кажется невозможным.

Тем не менее, я буду толковать ваш вопрос так, чтобы он имел больше смысла и мог быть решен. Что вам нужно, это распределение вероятностей на семимерной гиперплоскости x_1 + x_2 + ... + x_8 = 0.5, Поскольку гиперплоскость имеет бесконечную протяженность, равномерное распределение по всей гиперплоскости не будет работать. Что вы, вероятно (?) Хотите, это кусок гиперплоскости, где все x_i>0, Эта область является симплексом, обобщением треугольника, и равномерное распределение на симплексе является частным случаем распределения Дирихле.

Вы можете найти этот раздел статьи Dikichlet Distribution Wikipedia, вырезание струн, особенно освещение.

Реализация

Статья Википедии дает следующую реализацию на Python в разделе Генерация случайных чисел:

params = [a1, a2, ..., ak]
sample = [random.gammavariate(a,1) for a in params]
sample = [v/sum(sample) for v in sample]

То, что вы, вероятно, (?) Хотите, это случай, когда все ai=1 что приводит к равномерному распределению на симплексе. Вот k соответствует номеру N в вашем вопросе. Чтобы получить образцы для суммирования M вместо 1Просто умножьте sample от M,

Обновить

Спасибо Северину Паппадо за то, что он указал, что гаммавариант может возвращать бесконечность при редких обстоятельствах. Это математически "невозможно", но может возникать как артефакт реализации в терминах чисел с плавающей запятой. Мое предложение для обработки этого случая, чтобы проверить его после sample сначала рассчитывается; если какой-либо из компонентов sample равны бесконечности, установите все компоненты, не являющиеся бесконечными, равными 0, и установите все компоненты бесконечности равными 1. Затем, когда xi рассчитываются, результаты как xi=1, all other x's=0, или же xi=1/2, xj=1/2, all other x's=0 в результате все вместе "угловые образцы" и "граничные образцы".

Другая возможность с очень низкой вероятностью заключается в переполнении суммы гамма-вариатов. Я предположил бы, что мы могли бы пройти через всю последовательность псевдослучайных чисел и не увидеть, что это произойдет, но теоретически это возможно (в зависимости от базового генератора псевдослучайных чисел). Ситуация может быть решена путем изменения масштаба sampleнапример, деление всех элементов sample от N, после того, как гамма-вариации были рассчитаны, но до вычисления х. Лично я не стал бы беспокоиться, потому что шансы так низки; сбой программы по другим причинам будет иметь более высокую вероятность.

Answer 2

Вместо выбора "n" чисел из равномерного распределения, которые суммируются с "M", мы можем выбрать "n-1" случайный интервал из равномерного распределения с диапазоном "0-M", затем мы можем извлечь интервалы.

from random import uniform as rand

def randConstrained(n, M):
     splits = [0] + [rand(0, 1) for _ in range(0,n-1)] + [1]
     splits.sort()
     diffs = [x - splits[i - 1] for i, x in enumerate(splits)][1:]
     result = map(lambda x:x*M, diffs)
     return result

res = randConstrained(8,0.5)
print res
print sum(res)

Выход

[0.0004411388173262698,
 0.014832306834761111,
 0.009695872790939863,
 0.04539251424140245,
 0.18791325446494067,
 0.07615024971405443,
 0.07792489571128014,
 0.08764976742529507]

0.5

Answer 3

То же самое, что и решение k4vin, но без ошибки импорта, которую я получаю на random.uniform.

import random

def rand_constrained(n, total):
    # l is a sorted list of n-1 random numbers between 0 and total
    l = sorted([0] + [total * random.random() for i in range(n - 1)] + [total])
    # Return the intervals between each successive element
    return [l[i + 1] - l[i] for i in range(n)]

print(rand_constrained(3, 10))
# [0.33022261483938276, 8.646666440311822, 1.0231109448487956]

Но matplotlib задыхается от установки, поэтому я не могу подготовить ее прямо сейчас.

Answer 4

Это известно как симплексная выборка, которая тесно связана с распределением Дирихле. Sum(x_i) = 1, где каждый x_i равен U(0,1). В вашем случае после симплексной выборки просто умножьте каждый x_i на 0.5.

В любом случае, перевод кода C++ с https://github.com/Iwan-Zotow/SimplexSampling на python (надеюсь, не слишком много ошибок)

И это обращается с бесконечностью как раз

def simplex_sampling(n):
    r = []
    sum = 0.0
    for k in range(0,n):
        x = random.random()
        if x == 0.0:
            return (1.0, make_corner_sample(n, k))

        t = -math.log(x)
        r.append(t)
        sum += t

     return (sum, r)

def make_corner_sample(n, k):
    r = []
    for i in range(0, n):
        if i == k:
            r.append(1.0)
        else:
            r.append(0.0)

    return r

 # main
 sum, r = simplex_sampling(8)

 norm = 0.5 / sum # here is your 0.5 total

 for k in range(0, 8):
     r[k] *= norm

Answer 5

Для полностью обобщенного решения ("Я хочу n числа между low а также high, эта сумма к m):

from random import uniform as rand

def randConstrained(n, m, low, high):
    tot = m
    if not low <= 0 <= high:
        raise ValueError("Cannot guarantee a solution when the input does not allow for 0s")
    answer = []
    for _ in range(n-1):
        answer.append(low + rand(0,tot) * (high-low))
        tot -= answer[-1]
    answer.append(m-sum(answer))
    return answer

Для вашего случая это можно использовать следующим образом:

In [35]: nums = randConstrained(8, 0.5, 0, 1)

In [36]: nums
Out[36]: 
[0.2502590281277123,
 0.082663797709837,
 0.14586995648173873,
 0.011270073049224807,
 0.009328970756471237,
 0.00021993111786291258,
 0.0001831479074098452,
 0.000205094849743237]