Ответ 1

Если вы готовы приблизиться к p на основе 256 возможных значений, и у вас есть PRNG, который может генерировать однородные значения, в которых отдельные биты независимы друг от друга, тогда вы можете использовать векторизованное сравнение для создания нескольких смещенных бит из одного случайного числа.

Это стоит делать, если (1) вы беспокоитесь о качестве случайных чисел и (2) вам, вероятно, потребуется большое количество бит с одинаковым смещением. Второе требование, по-видимому, вытекает из первоначального вопроса, который критикует предлагаемое решение следующим образом: "Недостатком этого решения является то, что он может генерировать только 8 бит одновременно, даже при большой работе, в то время как непредвзятый PRNG может генерировать 64 одновременно с помощью всего лишь нескольких арифметических инструкций". Здесь подразумевается, что полезно создать большой блок смещенных битов за один вызов.

Качество случайных чисел - сложный вопрос. Трудно, если не невозможно измерить, и поэтому разные люди будут предлагать разные показатели, которые подчеркивают и/или девальвируют различные аспекты "случайности". Как правило, можно обменять скорость генерации случайных чисел на более низкое "качество"; стоит ли этого делать, зависит от вашего точного применения.

Простейшие возможные тесты качества случайных чисел связаны с распределением отдельных значений и продолжительностью цикла генератора. Стандартные реализации функций библиотеки C rand и Posix random, как правило, передают тест распределения, но длины циклов не подходят для долгосрочных приложений.

Эти генераторы, как правило, очень быстрые, однако: для реализации glibc random требуется всего несколько циклов, тогда как классический линейный конгруэнтный генератор (LCG) требует умножения и добавления. (Или, в случае реализации glibc, три из вышеперечисленных, чтобы генерировать 31 бит.) Если этого достаточно для ваших требований к качеству, то нет смысла пытаться оптимизировать, особенно если вероятность смещения часто изменяется.

Имейте в виду, что длина цикла должна быть намного больше, чем количество ожидаемых образцов; в идеале он должен быть больше квадрата этого числа, поэтому линейно-конгруэнтный генератор (LCG) с длиной цикла 2 ³¹ не подходит, если вы ожидаете генерировать гигабайты случайных данных. Даже триномиальный нелинейный генератор обратной связи Gnu, длина цикла которого составляет примерно 2 ³⁵ не должна использоваться в приложениях, для которых потребуются миллионы образцов.

Другая проблема качества, которая намного сложнее проверить, связана с независимостью в последовательных выборках. Короткие длины цикла полностью терпят неудачу по этой метрике, потому что как только начинается повторение, генерируемые случайные числа точно коррелируют с историческими значениями. Триномиальный алгоритм Gnu, хотя его цикл длиннее, имеет четкую корреляцию в результате того, что генерируемое случайное число я ^th r _i всегда является одним из два значения r _{я & minus; 3} & plus; r _{я & minus; 31} или r _{я & minus; 3} & plus; r _{я & minus; 31суб > & плюс; 1. У этого могут быть удивительные или, по крайней мере, загадочные последствия, особенно с экспериментами Бернулли.}

Здесь реализация, использующая Agner Fog, полезную библиотеку векторных классов, которая абстрагирует многие досадные детали в SSE-intrinsics, а также помогает поставляется с быстрым векторизованным генератором случайных чисел (находится в special.zip внутри архива vectorclass.zip), что позволяет нам генерировать 256 бит из восьми вызывает 256-битный PRNG. Вы можете прочитать д-р Фог, объяснив, почему он считает, что даже Мерседенский твистер имеет проблемы с качеством и его предлагаемое решение; Я не умею комментировать, правда, но он, по крайней мере, кажется, дает ожидаемые результаты в экспериментах Бернулли, которые я пробовал с ним.

#include "vectorclass/vectorclass.h"
#include "vectorclass/ranvec1.h"

class BiasedBits {
  public:
    // Default constructor, seeded with fixed values
    BiasedBits() : BiasedBits(1)  {}
    // Seed with a single seed; other possibilities exist.
    BiasedBits(int seed) : rng(3) { rng.init(seed); }

    // Generate 256 random bits, each with probability `p/256` of being 1.
    Vec8ui random256(unsigned p) {
      if (p >= 256) return Vec8ui{ 0xFFFFFFFF };
      Vec32c output{ 0 };
      Vec32c threshold{ 127 - p };
      for (int i = 0; i < 8; ++i) {
        output += output;
        output -= Vec32c(Vec32c(rng.uniform256()) > threshold);
      }
      return Vec8ui(output);
    }

  private:
    Ranvec1 rng;
};

В моем тесте это произвело и подсчитало 268435456 бит в 260 мс или один бит на наносекунду. Тест-машина - i5, поэтому у нее нет AVX2; YMMV.

В фактическом варианте использования, с 201 возможными значениями для p, вычисление 8-битных пороговых значений будет досадно неточным. Если эта неточность нежелательна, вы можете адаптировать приведенное выше значение для использования 16-битных пороговых значений за счет генерации в два раза больше случайных чисел.

В качестве альтернативы вы можете вручную нарисовать векторизацию на основе 10-битных пороговых значений, что даст вам очень хорошее приближение к шагом 0,5%, используя стандартную манипуляцию с использованием бит-манипуляций для выполнения векторизованного порогового сравнения путем проверки на заимствование на каждый 10-й бит вычитания вектора значений и повторного порога. В сочетании с, скажем, std::mt19937_64, это даст вам в среднем по шесть бит каждого 64-битного случайного числа.

Ответ 2

Одна вещь, которую вы можете сделать, состоит в том, чтобы сделать выборку из базового несмещенного генератора несколько раз, получить несколько 32-битных или 64-битных слов, а затем выполнить побитовую булевскую арифметику. Например, для 4 слов b1,b2,b3,b4 вы можете получить следующие дистрибутивы:

    expression             | p(bit is 1)
    -----------------------+-------------
    b1 & b2 & b3 & b4      |  6.25%
    b1 & b2 & b3           | 12.50%
    b1 & b2 & (b3 | b4)    | 18.75%
    b1 & b2                | 25.00%
    b1 | (b2 & (b3 | b4))  | 31.25%
    b1 & (b2 | b3)         | 37.50%
    b1 & (b2 | b3 | b4))   | 43.75%
    b1                     | 50.00%

Аналогичные конструкции могут быть сделаны для более точных разрешений. Он становится немного утомительным и требует больше вызовов генератора, но, по крайней мере, не один бит. Это похоже на ответ a3f, но, вероятно, проще реализовать и, я подозреваю, быстрее, чем сканировать слова для 0xF nybbles.

Обратите внимание, что для желаемого 0,5% -ного разрешения вам понадобится 8 несмещенных слов для одного смещенного слова, что даст вам разрешение (0.5 ^ 8) = 0.390625%.

Ответ 3

С теоретико-информационной точки зрения смещенный поток бит (с p != 0.5) имеет в нем меньше информации, чем несмещенный поток, поэтому теоретически он должен принимать (в среднем) менее 1 бит непредвзятого ввод для создания одного бита смещенного выходного потока. Например, entropy случайной переменной Bernoulli с p = 0.1 является битами -0.1 * log2(0.1) - 0.9 * log2(0.9), который находится вокруг бит 0.469. Это говорит о том, что для случая p = 0.1 мы должны иметь возможность создавать чуть более двух битов выходного потока на несмещенный входной бит.

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

Метод 1: арифметическое (де) кодирование

Практический метод состоит в том, чтобы декодировать ваш непредвзятый входной поток, используя арифметическое (де) кодирование, как уже описано в ответ от alexis. Для этого простого случая не сложно что-то подделать. Вот какой-то неоптимизированный псевдокод (кашель, Python), который делает это:

import random

def random_bits():
    """
    Infinite generator generating a stream of random bits,
    with 0 and 1 having equal probability.
    """
    global bit_count  # keep track of how many bits were produced
    while True:
        bit_count += 1
        yield random.choice([0, 1])

def bernoulli(p):
    """
    Infinite generator generating 1-bits with probability p
    and 0-bits with probability 1 - p.
    """
    bits = random_bits()

    low, high = 0.0, 1.0
    while True:
        if high <= p:
            # Generate 1, rescale to map [0, p) to [0, 1)
            yield 1
            low, high = low / p, high / p
        elif low >= p:
            # Generate 0, rescale to map [p, 1) to [0, 1)
            yield 0
            low, high = (low - p) / (1 - p), (high - p) / (1 - p)
        else:
            # Use the next random bit to halve the current interval.
            mid = 0.5 * (low + high)
            if next(bits):
                low = mid
            else:
                high = mid

Здесь пример использования:

import itertools
bit_count = 0

# Generate a million deviates.
results = list(itertools.islice(bernoulli(0.1), 10**6))

print("First 50:", ''.join(map(str, results[:50])))
print("Biased bits generated:", len(results))
print("Unbiased bits used:", bit_count)
print("mean:", sum(results) / len(results))

Вышеприведенный пример дает следующий пример:

First 50: 00000000000001000000000110010000001000000100010000
Biased bits generated: 1000000
Unbiased bits used: 469036
mean: 0.100012

Как и было обещано, мы генерировали 1 миллион бит нашего выходного смещенного потока, используя менее пятисот тысяч из исходного несмещенного потока.

В целях оптимизации при переводе на C/С++ может возникнуть смысл кодировать это, используя арифметику с фиксированной точкой на основе целого числа, а не с плавающей точкой.

Метод 2: алгоритм на основе целых чисел

Вместо того, чтобы пытаться преобразовать метод арифметического декодирования для непосредственного использования целых чисел, здесь используется более простой подход. Это не совсем арифметическое декодирование, но оно не полностью несвязано, и оно достигает близкого к тому же соотношению выходного смещения-бит/вход-несмещенное-бит, как версия с плавающей запятой выше. Он организован так, чтобы все величины вписывались в 32-битное целое без знака, поэтому его легко перевести на C/С++. Код специализирован для случая, когда p является точным кратным 1/200, но этот подход будет работать для любого p, который может быть выражен как рациональное число с достаточно малым знаменателем.

def bernoulli_int(p):
    """
    Infinite generator generating 1-bits with probability p
    and 0-bits with probability 1 - p.

    p should be an integer multiple of 1/200.
    """
    bits = random_bits()
    # Assuming that p has a resolution of 0.05, find p / 0.05.
    p_int = int(round(200*p))

    value, high = 0, 1
    while True:
        if high < 2**31:
            high = 2 * high
            value = 2 * value + next(bits)
        else:
            # Throw out everything beyond the last multiple of 200, to
            # avoid introducing a bias.
            discard = high - high % 200
            split = high // 200 * p_int
            if value >= discard:  # rarer than 1 time in 10 million
                value -= discard
                high -= discard
            elif value >= split:
                yield 0
                value -= split
                high = discard - split
            else:
                yield 1
                high = split

Главное наблюдение заключается в том, что каждый раз, когда мы достигаем начала цикла while, value равномерно распределяется среди всех целых чисел в [0, high) и не зависит от всех ранее выведенных битов. Если вы заботитесь о скорости больше, чем совершенная правильность, вы можете избавиться от ветки discard и value >= discard: именно там, чтобы гарантировать, что мы выводим 0 и 1 точно с правильными вероятностями. Оставьте это осложнение, и вместо этого вы получите почти правильные вероятности. Кроме того, если вы сделаете разрешение для p равным 1/256, а не 1/200, то потенциально трудоемкие операции деления и по модулю могут быть заменены битовыми операциями.

С тем же тестовым кодом, что и раньше, но используя bernoulli_int вместо bernoulli, я получаю следующие результаты для p=0.1:

First 50: 00000010000000000100000000000000000000000110000100
Biased bits generated: 1000000
Unbiased bits used: 467997
mean: 0.099675

Ответ 4

Скажем, вероятность появления 1 равна 6,25% (1/16). Существует 16 возможных битовых шаблонов для 4-битного номера: 0000,0001, ..., 1110,1111.

Теперь просто сгенерируйте случайное число, как вы привыкли, и замените каждый 1111 на границе полубайта с помощью 1 и превратите все остальное в 0.

Корректировка соответственно для других вероятностей.

Ответ 5

Вы получите теоретически оптимальное поведение, то есть сделайте действительно минимальное использование генератора случайных чисел и сможете точно моделировать вероятность p,, если вы подходите к этому, используя арифметическое кодирование.

Арифметическое кодирование представляет собой форму сжатия данных, которая представляет сообщение как промежуток интервала диапазона. Он обеспечивает теоретически оптимальное кодирование и может использовать дробное число бит для каждого входного символа.

Идея такова: представьте, что у вас есть последовательность случайных бит, которые равны 1 с вероятностью p. Для удобства я вместо этого использую q для вероятности того, что бит равен нулю. (q = 1-p). Арифметическое кодирование присваивается каждой битовой части диапазона чисел. Для первого бита назначьте интервал [0, q), если входной сигнал равен 0, а интервал [q, 1), если входной сигнал равен 1. Последующие биты назначают пропорциональные вспомогательные интервалы текущего диапазона. Например, предположим, что q = 1/3. Вход 1 0 0 будет закодирован следующим образом:

Initially       [0, 1),             range = 1
After 1         [0.333, 1),         range = 0.6666        
After 0         [0.333, 0.5555),    range = 0.2222   
After 0         [0.333, 0.407407),  range = 0.074074

Первая цифра, 1, выбирает верхние две трети (1-q) диапазона; вторая цифра, 0, выбирает нижнюю треть этого и т.д. После первого и второго шага интервал перемещается по средней точке; но после третьего шага он полностью находится ниже середины, поэтому может быть выведена первая сжатая цифра: 0. Процесс продолжается, и в качестве терминатора добавляется специальный символ EOF.

Чем это связано с вашей проблемой? Сжатый вывод будет иметь случайные нули и единицы с равной вероятностью. Итак, чтобы получить биты с вероятностью p, просто притворитесь, что вывод вашего RNG является результатом арифметического кодирования, как указано выше, и применить к нему процесс декодирования. То есть, читать биты, как если бы они подразделялись интервал линии на меньшие и меньшие части. Например, после того, как мы прочитаем 01 из RNG, мы будем находиться в диапазоне [0,25, 0,5]. Храните считывающие биты до тех пор, пока достаточный выход не будет "декодирован". Поскольку вы подражаете распаковке, вы получите больше случайных бит, чем вы вставляете.. Поскольку арифметическое кодирование является теоретически оптимальным, нет никакого способа превратить вывод RNG в более предвзятые биты, не жертвуя случайностью: вы получаете максимальный максимум.

Уловка заключается в том, что вы не можете сделать это в нескольких строках кода, и я не знаю библиотеки, на которую я могу указать (хотя, возможно, некоторые из них вы можете использовать). Тем не менее, это довольно просто. выше статьи содержит код для кодера и декодера общего назначения, в C. Это довольно просто и поддерживает несколько входных символов с произвольными вероятностями; в вашем случае возможно гораздо более простая реализация (как ответ Марка Дикинсона), так как вероятностная модель тривиальна. Для расширенного использования потребуется немного больше работы для создания надежной реализации, которая не выполняет много вычислений с плавающей запятой для каждого бита.

Wikipedia также имеет интересное обсуждение арифметического кодирования, рассматриваемого как изменение radix, что является другим способом просмотра вашей задачи.

Ответ 6

Uh, генераторы псевдослучайных чисел, как правило, довольно быстрые. Я не уверен, что это за язык (возможно, Python), но "result.append" (который почти наверняка содержит выделение памяти), скорее всего, медленнее, чем "random_uniform" (что делает небольшую математику).

Если вы хотите оптимизировать производительность этого кода:

• Убедитесь, что это проблема. Оптимизации - это небольшая работа и сложность работы с кодом. Не делайте их, если это необходимо.
• Профиль. Запустите некоторые тесты, чтобы определить, какие части кода на самом деле самые медленные. Это те части, которые вам нужно ускорить.
• Внесите свои изменения и убедитесь, что они на самом деле быстрее. Компиляторы довольно умны; часто чистый код будет скомпилирован в лучший код, что-то сложное, чем может показаться быстрее.

Если вы работаете на компилированном языке (даже с JIT-компилятором), вы получаете удар производительности для каждой передачи управления (если, во время вызова функции и т.д.). Устраните, что вы можете. Распределение памяти также (обычно) довольно дорогое.

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

Я могу только догадываться, где ваши проблемы с производительностью:

• Распределение памяти. Предварительно выделите массив по своему размеру и заполните записи позже. Это гарантирует, что память не будет перераспределена во время добавления записей.
• Отрасли. Возможно, вы сможете избежать "если", произведя результат или что-то подобное. Это будет сильно зависеть от компилятора. Проверьте сборку (или профиль), чтобы убедиться, что она делает то, что вы хотите.
• Числовые типы. Узнайте, какой тип генератор случайных чисел использует изначально, и выполните свою арифметику в этом типе. Например, если генератор естественно возвращает 32-разрядные целые числа без знака, сначала масштабируйте "p" до этого диапазона, а затем используйте его для сравнения.

Кстати, если вы действительно хотите использовать наименьшие бит случайности, используйте "арифметическое кодирование" для декодирования вашего случайного потока. Это не будет быстро.

Ответ 7

Один из способов, который дал бы точный результат, - сначала случайным образом генерировать для k-битового блока число 1 бита, следующего за биномиальным распределением, а затем сгенерировать k-битовое слово с таким количеством бит, используя один из методов здесь. Например, метод mic006 требует только о log k k-битных случайных числах, а my - только один.

Ответ 8

Если p близко к 0, вы можете вычислить вероятность того, что n-й бит является первым битом, равным 1; то вы вычисляете случайное число между 0 и 1 и выбираете n соответственно. Например, если p = 0,005 (0,5%), а случайное число - 0,638128, вы можете вычислить (я угадываю здесь) n = 321, поэтому вы заполняете 321 0 бит и один бит.

Если p близок к 1, используйте 1-p вместо p и установите 1 бит плюс один 0 бит.

Если p не близко к 1 или 0, сделайте таблицу из всех 256 последовательностей из 8 бит, вычислите их кумулятивные вероятности, затем получите случайное число, выполните двоичный поиск в массиве кумулятивных вероятностей, и вы можете установите 8 бит.

Ответ 9

Предполагая, что у вас есть доступ к генератору случайных битов, вы можете сгенерировать значение для сравнения с p по биту и прервать, как только вы сможете доказать, что сгенерированное значение меньше или больше или -equal-to p.

Выполните следующие действия, чтобы создать один элемент в потоке с заданной вероятностью p:

• Начните с 0. в двоичном формате
• Добавить случайный бит; предполагая, что a 1 был вычерчен, вы получите 0.1
• Если результат (в двоичной нотации) доказуемо меньше, чем p, выведите a 1
• Если результат доказуемо больше или равен p, выведите a 0
• В противном случае (если ни одно из них не может быть исключено), перейдите к шагу 2.

Предположим, что p в двоичной нотации 0.1001101...; если этот процесс генерирует любой из 0.0, 0.1000, 0.10010,..., значение больше не может быть больше или равно p; если генерируется какой-либо из 0.11, 0.101, 0.100111,..., значение не может стать меньше p.

Мне кажется, что этот метод использует около двух случайных бит в ожидании. Арифметическое кодирование (как показано в ответе Марка Дикинсона) потребляет не более одного случайного бита на смещенный бит (в среднем) для фиксированного p; стоимость модификации p неясна.

Ответ 10

Что он делает

Эта реализация делает вызов single на случайный модуль ядра устройства через интерфейс специального символьного файла "/dev/urandom", чтобы получить количество случайных данных, необходимых для представления всех значений в заданном разрешении. Максимально возможное разрешение 1/256 ^ 2, поэтому 0,005 может быть представлено:

328/256 ^ 2,

то есть:

разрешение: 256 * 256

x: 328

с ошибкой 0.000004883.

Как это делает

Реализация вычисляет количество бит bits_per_byte, которое представляет собой число равномерно распределенных битов, необходимых для обработки заданного разрешения, то есть представляют все значения @resolution. Затем он делает один вызов устройству рандомизации ( "/dev/urandom", если URANDOM_DEVICE определен, в противном случае он будет использовать дополнительный шум от драйверов устройств по вызову "/dev/random", который может блокироваться, если энтропии недостаточно в битах), чтобы получить необходимое количество равномерно распределенных байтов и заполнить массив rnd_bytes случайных байтов. Наконец, он считывает количество необходимых бит на каждый образец Бернулли из каждого байта bytes_per_byte массива rnd_bytes и сравнивает целочисленное значение этих битов с вероятностью успеха в одном результате Бернулли, заданное x/resolution. Если значение попадает, т.е. Попадает в сегмент длины x/resolution, который мы произвольно выбираем как сегмент [0, x/resolution], тогда мы отмечаем успех и вставляем 1 в результирующий массив.

Чтение с произвольного устройства:

/* if defined use /dev/urandom (will not block),
 * if not defined use /dev/random (may block)*/
#define URANDOM_DEVICE 1

/*
 * @brief   Read @outlen bytes from random device
 *          to array @out.
 */
int
get_random_samples(char *out, size_t outlen)
{
    ssize_t res;
#ifdef URANDOM_DEVICE
    int fd = open("/dev/urandom", O_RDONLY);
    if (fd == -1) return -1;
    res = read(fd, out, outlen);
    if (res < 0) {
        close(fd);
        return -2;
    }
#else
    size_t read_n;
    int fd = open("/dev/random", O_RDONLY);
    if (fd == -1) return -1;
    read_n = 0;
    while (read_n < outlen) {
       res = read(fd, out + read_n, outlen - read_n);
       if (res < 0) {
           close(fd);
           return -3;
       }
       read_n += res;
    }
#endif /* URANDOM_DEVICE */
    close(fd);
    return 0;
}

Заполните вектор образцов Бернулли:

/*
 * @brief   Draw vector of Bernoulli samples.
 * @details @x and @resolution determines probability
 *          of success in Bernoulli distribution
 *          and accuracy of results: p = x/resolution.
 * @param   resolution: number of segments per sample of output array 
 *          as power of 2: max resolution supported is 2^24=16777216
 * @param   x: determines used probability, x = [0, resolution - 1]
 * @param   n: number of samples in result vector
 */
int
get_bernoulli_samples(char *out, uint32_t n, uint32_t resolution, uint32_t x)
{
    int res;
    size_t i, j;
    uint32_t bytes_per_byte, word;
    unsigned char *rnd_bytes;
    uint32_t uniform_byte;
    uint8_t bits_per_byte;

    if (out == NULL || n == 0 || resolution == 0 || x > (resolution - 1))
        return -1;

    bits_per_byte = log_int(resolution);
    bytes_per_byte = bits_per_byte / BITS_PER_BYTE + 
                        (bits_per_byte % BITS_PER_BYTE ? 1 : 0);
    rnd_bytes = malloc(n * bytes_per_byte);
    if (rnd_bytes == NULL)
        return -2;
    res = get_random_samples(rnd_bytes, n * bytes_per_byte);
    if (res < 0)
    {
        free(rnd_bytes);
        return -3;
    }

    i = 0;
    while (i < n)
    {
        /* get Bernoulli sample */
        /* read byte */
        j = 0;
        word = 0;
        while (j < bytes_per_byte)
        {
            word |= (rnd_bytes[i * bytes_per_byte + j] << (BITS_PER_BYTE * j));
            ++j;
        }
        uniform_byte = word & ((1u << bits_per_byte) - 1);
        /* decision */
        if (uniform_byte < x)
            out[i] = 1;
        else
            out[i] = 0;
        ++i;
    }

    free(rnd_bytes);    
    return 0;
}

Использование:

int
main(void)
{
    int res;
    char c[256];

    res = get_bernoulli_samples(c, sizeof(c), 256*256, 328); /* 328/(256^2) = 0.0050 */
    if (res < 0) return -1;

    return 0;
}

Полный код, результаты.