Генерация случайных чисел при очень специфических ограничениях

Ответ 1

Начните с создания как можно большего количества идентичных кортежей n, если вам нужно:

(A/n, B/n)

Теперь выберите два кортежа в случайном порядке. Сделайте случайное изменение значения a одного и компенсирующего изменения значения a другого, сохраняя все в пределах заданных ограничений. Верните два кортежа.

Теперь выберите другую случайную пару. Это время меняется со значениями b.

Скорее, повторите полоскание.

Ответ 2

Я думаю, что самая простая вещь -

• Используйте свой любимый метод для сброса значений n-1, таких как \sum_i=0,n-1 a_i < A, и установите a_n, чтобы получить нужную сумму. Есть несколько вопросов о том, как это сделать, хотя я никогда не видел ответа, которого я действительно доволен. Может быть, я напишу статью или что-то еще.

• Получите n-1 b, выставив c_i равномерно на допустимый диапазон и установите final b, чтобы получить нужную сумму и проверить окончательный c (я думаю, что это должно быть ОК, но я еще не доказал этого).

Заметим, что, поскольку у нас есть 2 жестких ограничения, мы должны ожидать бросить 2n-2 случайные числа, и этот метод делает именно это (в предположении, что вы можете сделать шаг 1 с n-1 throws.

Ответ 3

Мы ищем кортежи a_i и b_i такие, что

• (a_1,... a_n) и (b_1,... b_n) имеют распределение, которое инвариантно относительно перестановки индексов (что вы бы назвали "несмещенным" )
• отношения a_i/b_i равномерно распределены на [cmin, cmax]
• sum (a_i) = A, sum (b_i) = B

Если c_min и c_max не слишком плохо обучены (т.е. они не очень близки к другому), а n не очень велико, то работает следующее:

• Создать a_i равномерно, чтобы sum a_i = A:
  • Нарисуйте n образцы aa_i (i = 1..n) из некоторого распределения (например, равномерного)
  • Разделите их по их сумме и умножьте на A: a_i = A * aa_i / sum(aa_i) имеет желаемые свойства.
• Создайте b_i таким образом, чтобы sum b_i = B одним и тем же методом.
• Если существует i такое, что a_i / b_i не находится в интервале [cmin, cmax], выбросьте все a_i и b_i и повторите попытку с самого начала.

Он не имеет хорошо с n, потому что набор a_i и b_i, удовлетворяющий ограничениям, становится все более узким по мере увеличения n (и поэтому вы отклоняете больше кандидатов).

Честно говоря, я не вижу другого простого решения. Если n становится большим и cmin ~ cmax, вам нужно будет использовать кувалду (например, MCMC) для генерации выборок из вашего дистрибутива, если не будет какой-то трюк, который мы не видели.

Если вы действительно хотите использовать алгоритмы MCMC, обратите внимание, что вы можете изменить cmin на cmin * B / A (аналогично для cmax) и принять A == B == 1. Тогда задача состоит в том, чтобы равномерно нанести на произведение двух единичных n-симплексов (u_1... u_n, v_1 ​​... v_n) таких, что

u_i / v_i \in [cmin, cmax].

Таким образом, вам нужно использовать алгоритм MCMC (более подходящий для Metropolis-Hastings) для произведения двух единичных n-симплексов с плотностью

f(u_1, ..., u_n, v_1, ..., v_n) = \prod indicator_{u_i/v_i \in [cmin, cmax]}

который определенно выполним (хотя и задействован).

Ответ 4

Забитая выборка Gibbs довольно проста и сходится к правильному распределению (это по строкам того, что предлагает Alexandre).

• Для всех я инициализируйте a _i= A/n и b _i= B/n.
• Выберите я ≠ j равномерно случайным образом. С вероятностью 1/2 обновите _i и _j с равномерными случайными значениями, удовлетворяющими ограничениям. В остальное время сделайте то же самое для b _i и b _j.
• Повторите шаг 2 столько раз, сколько необходимо для вашего приложения. Я понятия не имею, что такое коэффициент конвергенции.

Ответ 5

Итак, вот что я думаю с математической точки зрения. Имеются последовательности a_i и b_i такие, что сумма a_i равна A, а сумма b_i равна B. Более того, A/B находится в (x,y), и поэтому a_i/b_i для каждого i. Кроме того, вы хотите, чтобы a_i/b_i был равномерно распределен в (x,y).

Так сделайте это, начиная с конца. Выберите c_i из (x,y) так, чтобы они были равномерно распределены. Тогда мы хотим иметь следующее равенство a_i/b_i = c_i, поэтому a_i = b_i*c_i.

Поэтому нам нужно найти b_i. Но мы имеем следующую систему линейных уравнений:

A = (sum)b_i*c_i
B = (sum)b_i

где b_i - переменные. Решите его (некоторые причудливые уловки линейной алгебры), и все готово!

Обратите внимание, что для достаточно большой n эта система будет иметь множество решений. Они будут зависеть от некоторых параметров, которые вы можете выбрать произвольно.

Достаточно теоретического подхода, рассмотрим практическое решение.

//РЕДАКТИРОВАТЬ 1: Здесь находится код жесткого ядра Python: D

import random
min = 0.0
max = 10.0
A = 500.0
B = 100.0

def generate(n):
    C = [min + i*(max-min)/(n+1) for i in range(1, n+1)]
    Y = [0]
    for i in range(1,n-1):
        # This line should be changed in order to always get positive numbers
        # It should be relatively easy to figure out some good random generator
        Y.append(random.random())
    val = A - C[0]*B
    for i in range(1, n-1):
        val -= Y[i] * (C[i] - C[0])
    val /= (C[n-1] - C[0])
    Y.append(val)
    val = B
    for i in range(1, n):
        val -= Y[i]
    Y[0] = val
    result = []
    for i in range(0, n):
        result.append([ Y[i]*C[i], Y[i] ])
    return result

Результат представляет собой список пар (x,y), удовлетворяющий вашим условиям, за исключением того, что они могут быть отрицательными (см. строку случайного генератора в коде), то есть первая и последняя пары могут содержать отрицательные числа.

//ИЗМЕНИТЬ 2:

Чтобы убедиться, что они положительные, вы можете попробовать что-то вроде

Y.append(random.random() * B / n)

вместо

Y.append(random.random())

Я не уверен, хотя.

//EDIT 3:

Чтобы получить лучшие результаты, попробуйте что-то вроде этого:

avrg = B / n
ran = avrg / 20
for i in range(1, n-1):
    Y.append(random.gauss(avrg, ran))

вместо

for i in range(1, n-1):
    Y.append(random.random())

Это сделает все b_i ближе к B / n. К сожалению, последний термин по-прежнему будет прыгать высоко. Извините, но избежать этого невозможно (математика), поскольку последние и первые термины зависят от других. При малых n (~ 100) это выглядит неплохо. К сожалению, могут появиться некоторые отрицательные значения.

Выбор правильного генератора не так прост, если вы хотите, чтобы b_i был равномерно распределен.

Ответ 6

Здесь много хороших идей. Благодарю! Идея Rossum казалась самой простой в использовании, поэтому я пошел за ней. Вот код для потомков:

c_min = 0.25
c_max = 0.75
a_sum = 100.0
b_sum = 200.0
n = 1000 

a = [a_sum / n] * n
b = [b_sum / n] * n

while not good_enough(a, b):
    i, j = random.sample(range(n), 2)
    li, ui = c_min * b[i] - a[i], c_max * b[i] - a[i]
    lj, uj = a[j] - c_min * b[j], a[j] - c_max * b[j]
    llim = max((li, uj))
    ulim = min((ui, lj))
    q = random.uniform(llim, ulim)
    a[i] += q
    a[j] -= q

    i, j = random.sample(range(n), 2)
    li, ui = a[i] / c_max - b[i], a[i] / c_min - b[i]
    lj, uj = b[j] - a[j] / c_max, b[j] - a[j] / c_min
    llim = max((li, uj))
    ulim = min((ui, lj))
    q = random.uniform(llim, ulim)
    b[i] += q
    b[j] -= q

Функция good_enough(a, b) может быть очень много. Я пробовал:

• Стандартное отклонение, которое ударяется или промахивается, поскольку вы не знаете, что является достаточно хорошим значением.
• Куртозис, где большое отрицательное значение было бы приятным. Однако относительно медленно вычислять и undefined с начальными значениями (a_sum / n, b_sum / n) (хотя это тривиально для исправления).
• Скос, где значение, близкое к 0, желательно. Но он имеет те же недостатки, что и эксцесс.
• Число итераций, пропорциональных n. 2n иногда не хватало, n ^ 2 - это немного переполняющий и, следовательно, экспоненциальный.

В идеале эвристика, использующая сочетание асимметрии и эксцесса, была бы лучше, но я решил, что каждое значение было изменено с начального (опять же, как rossum, предложенного в комментарии). Хотя теоретической гарантии, что петля не будет завершена, не существует, казалось, что она работает достаточно хорошо для меня.