2013-10-03 2 views
8

Предположим, у нас есть генератор двоичных случайных чисел, int r();, который будет возвращать нуль или один с возможностью пролонгирования 0.5.Правильный способ генерации случайного поплавка с учетом двоичного генератора случайных чисел?

Я смотрел на Boost.Random, и они производят, скажем, 32 бита и сделать что-то вроде этого (псевдокод):

x = double(rand_int32()); 
return min + x/(2^32) * (max - min); 

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

Что бы быстрый способ создания равномерно распределенной float или double в полуоткрытый диапазон [min, max), предполагая IEEE754? Акцент здесь лежит на правильности распространения, а не скорости.

Для правильного определения правильности правильное распределение будет равным тому, которое мы получили бы, если бы мы взяли бесконечно точный равномерно распределенный генератор случайных чисел и для каждого числа, которое мы бы округлили до ближайшего представления IEEE754, если это представление все равно будет в пределах [min, max), в противном случае число не будет учитываться для распределения.

P.S .: Меня также интересовали бы правильные решения для открытых диапазонов.

+0

Я предполагаю, что «возьмите бесконечно точное равномерное генератор случайных чисел на диапазоне' [min, max) '", потому что без ограничений диапазона не существует равномерных генераторов случайных чисел. :) – Yakk

+0

@Yakk Это создало бы неправильное распределение на самых концах из-за округления, '[min - n * ULP, max + n * ULP)' было бы хорошо, хотя для некоторого значения 'n', что я слишком ленив, чтобы подумать, но, вероятно, это 1. – orlp

+0

@nightcracker Чтобы убедиться, что я правильно понимаю: вы хотите использовать dsitrbution, где некоторые битовые шаблоны с плавающей запятой гораздо более распространены, чем другие (из-за того, что они находятся дальше от 0), НЕ где каждая битовая диаграмма с плавающей запятой равновероятна? –

ответ

1

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

Таким образом, мы должны определить «equi-distribution». Тем не менее, считая, что каждый номер IEEE754 является просто представителем вероятности лежания в интервале, определяемом округлением IEEE754, процедура первого генерирования равнораспределенных «чисел» и округления до IEEE754 будет генерировать (по определению) equi-distribution "номеров IEEE754.

Следовательно, я считаю, что приведенная выше формула будет произвольной вблизи такого распределения, если мы просто выберем достаточно высокую точность. Если мы ограничим проблему поиском числа в [0,1], это означает ограничение набора деноминированных номеров IEEE 754, которые являются взаимно однозначными и 53-битными целыми числами. Таким образом, это должно быть быстрым и правильным, чтобы генерировать только мантис по 53-битовому генератору случайных чисел.

IEEE 754 арифметика всегда «арифметика с бесконечной точностью с последующим округлением», то есть число IEEE754, представляющий б это один будучи ближе к Ь (иначе говоря, вы можете думать о * Ь, вычисленной на бесконечное точность, затем округляется до номера IEEE754). Поэтому я считаю, что min + (max-min) * x, где x - деноминированное число, является допустимым подходом.

(Примечание. Как явствует из моего комментария, я сначала не знал, что вы указываете на случай, когда min и max отличаются от 0,1. Денормализованные числа обладают свойством, что они равномерно распределены. распределение equi путем сопоставления 53 бит с мантиссой. Затем вы можете использовать арифметику с плавающей запятой, из-за того, что она правильна до отказа машины. Если вы используете обратное сопоставление, вы восстановите равнораспределение.

См. Этот вопрос для другого аспекта этой проблемы: Scaling Int uniform random range into Double one

+0

Эта проблема именно поэтому я задал этот вопрос. Есть еще много номеров IEEE754, близких к 0, например, таких, которые близки к 1. – orlp

+0

Да. Но я хотел сказать следующее: номер IEEE 754 является представителем его интервала округления (в смысле отношений эквивалентности), а случайное число представляет собой событие «с реальным числом, падающим в этот интервал». Наделенное этим понятием, вы отбираете равноценное распределение. Просто так, что ваши интервалы выборки не равны друг другу, но это не наносит вреда тому факту, что вы выбрали равноценное распределение. –

3

Вот правильный подход без каких-либо усилий.

Начнем с класса бигума, а затем - рациональной обертки указанных бонусов.

Мы производим ряд «достаточно большой, чем» наш [min, max) диапазон, так что округление нашего smaller_min и bigger_max производит значения с плавающей точкой за пределами этого диапазона, в нашем рациональном построене на bignum.

Теперь мы подразделим диапазон на две части отлично по середине (что мы можем сделать, поскольку у нас есть рациональная система bignum). Мы выбираем одну из двух частей случайным образом.

Если после округления верх и низ выбранного диапазона будут (A) за пределами [min, max) (с той же стороны, заметьте!), Вы отклоняете и перезапускаете с самого начала.

Если (B) верхняя и нижняя части ваших рядов округлены до одинаковых double (или float, если вы возвращаете поплавок), все готово, и вы возвращаете это значение.

В противном случае (C) вы рекурсируете на этот новый, меньший диапазон (подразделяйте, выбирайте случайным образом, проверяйте).

Нет гарантий, что эта процедура останавливается, потому что вы можете постоянно переходить к «краю» между двумя закруглениями double s, или вы можете постоянно выбирать значения за пределами диапазона [min, max). Вероятность этого события (никогда не останавливается), однако, равна нулю (предполагая хороший генератор случайных чисел и [min, max) ненулевого размера).

Это также работает для (min, max) или даже подбирает номер в закругленном достаточно толстом наборе Кантора. Пока мера допустимого диапазона реалов, округленная до правильных значений с плавающей запятой, не равна нулю, а диапазон имеет компактную поддержку, эта процедура может быть выполнена и имеет вероятность 100% завершения, но не жесткая верхняя связанное с временем, которое может потребоваться.

+0

Звучит солидно. Есть ли алгоритм округления, который вы бы рекомендовали для рационального преобразования> IEEE754? – orlp

+0

@ nightcracker Округление будет работать (дайте правильное распределение): это действительно зависит от того, какой «регион» вы хотите, чтобы каждый IEEE754 представлял. Думаю, мы могли бы взглянуть на спецификацию IEEE754, чтобы выяснить, как она проходит, когда вы делаете что-то вроде умножения? Это дало бы нам понять, в каком регионе они «должны» представлять. – Yakk

+0

Арифметика IEEE754 всегда «арифметика с бесконечной точностью, за которой следует округление», т. Е. Число IEEE754, представляющее a * b, является самым близким к a * b (иначе говоря, вы можете думать о * b, вычисленном с бесконечной точностью, затем округленным к закрытию номера IEEE754). Поэтому я считаю, что min + (max-min) * x, где x - деноминированное число, является допустимым подходом. –

1

std::uniform_real_distribution.

У S.T.L. есть really good talk. начиная с конференции Going Native в этом году, которая объясняет, почему вы должны использовать стандартные дистрибутивы, когда это возможно. Короче говоря, ручной код имеет тенденцию быть смехотворно низкого качества (думаю, std::rand() % 100) или иметь более тонкие недостатки однородности, например, в (std::rand() * 1.0/RAND_MAX) * 99, что является примером, приведенным в разговоре, и является частным случаем кода, опубликованного в вопрос.

EDIT: Я посмотрел на libstdC++ реализации 's из std::uniform_real_distribution, и это то, что я нашел:

Реализация производит ряд в диапазоне [dist_min, dist_max), используя простое линейное преобразование из некоторого числа, полученного в диапазоне [0, 1) , Он генерирует этот номер источника, используя std::generate_canonical, the implementation of which my be found here (в конце файла). std::generate_canonical определяет количество раз (обозначенное как k) диапазон распределения, выраженный как целое число и обозначаемый здесь как r *, будет вписываться в мантиссу целевого типа. То, что он тогда делает, состоит в том, чтобы генерировать одно число в [0, r) для каждого сегмента r-сегмента мантиссы и, используя арифметику, заполнять каждый сегмент соответственно.Формула для результирующего значения может быть выражена как

Σ(i=0, k-1, X/(r^i)) 

, где X является случайной величиной в [0, r). Каждое деление на диапазон эквивалентно сдвигу на количество бит, используемых для его представления (то есть, log2(r)), и таким образом заполняет соответствующий сегмент мантиссы. Таким образом, используется вся точность целевого типа, и поскольку диапазон результатов равен [0, 1), показатель остается 0 ** (по модулю смещения), и вы не получаете проблем с однородностью, которые у вас возникают, когда вы начинаете беспорядок с показателем.

Я не буду верить, что этот метод криптографически безопасен (и у меня есть подозрения о возможных ошибках «один за другим» при вычислении размера r), но я полагаю, что он значительно более надежный с точки зрения единообразия чем вы делали Boost, и определенно лучше, чем возиться с std::rand.

Может быть, стоит отметить, что код подталкивания фактически является вырожденный случай этого алгоритма, где k = 1, а это означает, что она эквивалентна , если диапазон входного требует, по меньшей мере, 23 бита, чтобы представить его размер (ПЭО 754 одно- точность) или не менее 52 бит (с двойной точностью). Это означает минимальный диапазон ~ 8,4 миллиона или ~ 4,5e15, соответственно. В свете этой информации, я не думаю, что если вы используете двоичный генератор, то реализация Boost - это , а собирается разрезать ее.

После краткого обзора libc++’s implementation, похоже, что они используют один и тот же алгоритм, реализованный несколько иначе.

(*) r - это фактически диапазон ввода плюс один. Это позволяет использовать значение ёртга max в качестве допустимого ввода.

(**) Строго говоря, кодированный показатель не является 0, так как IEEE 754 кодирует неявное начало 1 перед основанием знака. Концептуально, однако, это не имеет отношения к этому алгоритму.

+3

-1 При всем моем уважении это совсем не полезно. В моем вопросе явно указано, что я не доверяю реализации Boost.Random, а 'std :: uniform_real_distribution' скопирован в значительной степени из Boost.Random. И 'std :: uniform_real_distribution' - это просто имя. Я ищу фактический алгоритм. – orlp

+0

@ nightcracker Интерфейс стандартной случайной библиотеки определенно сильно зависит от Boost, но это не обязательно означает, что реализация. В случае с libstdC++, похоже, используется 'std :: generate_canonical' для отображения вывода генератора на [0, 1) и линейное преобразование от него до [min, max). Реализация 'std :: generate_canonical' довольно интересна, и я увижу, могу ли я обновить объяснение. –

1

AFAIK, правильный (и, вероятно, самый быстрый) способ состоит в том, чтобы сначала создать 64-битное беззнаковое целое число, где 52 бит бит - это случайные биты, а показатель - 1023, который, если тип записывается в (IEEE 754) double будет равномерно распределенным случайным значением в диапазоне [1.0, 2.0]. Таким образом, последним шагом является вычитание 1.0 из этого, что приводит к равномерно распределенному случайному двойному значению в диапазоне [0.0, 1.0].

В псевдокоде:

rndDouble = bitCastUInt64ToDouble (1023 < < 52 | rndUInt6 0xfffffffffffff) - 1,0

Этот метод упоминается здесь: http://xoroshiro.di.unimi.it (см «Создание единых двойников в единичном интервале ")

EDIT: Рекомендованный метод сменил до: (x >> 11) * (1./(UINT64_C (1) < < 53))

См. Выше ссылку для получения более подробной информации.

+0

Просто захотелось добавить, что если вы хотите сгенерировать одиночный 32-битный float (что означает «float» означает в большинстве языков), формула будет выглядеть следующим образом: rndFloat = bitCastUInt32ToDouble (127 << 23 | rndUInt32 >> 9) - 1.0 (здесь, rndUInt32 >> 9 - это то же, что и rndUInt32 & 0x1ffffff) –

Смежные вопросы