Нужно мне реализовать метод Монте-Карло для подсчёта интегралов. Для этого нужны случайные числа с равномерным распределением.
Я так понял, стандартная функция rand() языка C не подходит, поскольку её алгоритм не обеспечивает идеально равномерного распределения. Хотелось бы, чтобы распределением было именно равномерным.
Ну и вообще, на будущее (может, пригодится, хотелось бы поэкспериментировать), не помешали бы и нормальное распределение, и другие распределения...
Как-то писать самому код, моделировать %) не хочется на это время тратить. Есть куча других дел. Ведь есть же наверное библиотеки, в которых всё уже есть, достаточно подключить и использовать?
Язык C, ну или C++ на крайний случай. Python не хочу... мне он не нравится. Другие языки вроде java тоже лучше не приплетать... (нет времени изучать что-то новое)
Нашёл библиотеку boost. Как вариант вообще. Но она для C++ только. А есть что для чистого C?
Но ещё смутила строка "Note: The current implementation is buggy, because it may not fill all of the mantissa with random bits. "
http://www.boost.org/doc/libs/1_45_0/doc/html/boost/uniform_real.html
Это, простите меня, неприемлемо. Проверю, конечно, может, интегралы считаться и будут правильно... но как так текущая реализация не работает правильно??? А смысл ещё вообще выкладывать тогда в общий доступ? %)(% туплю. Или это просто документация устаревшая может быть?
LinuxForum.Fans
Генераторы случайных чисел
Автор: Battle_Coder
7 сообщений
3750 просмотров
Ну и вообще, на будущее (может, пригодится, хотелось бы поэкспериментировать), не помешали бы и нормальное распределение, и другие распределения...
есть /dev/random для криптостойких СЧ, и /dev/urandom для не криптостойких.
вам нужен классический детерминированный ГПСЧ. читать здесь:
Дональд Э. Кнут. Глава 3. Случайные числа // Искусство программирования = The Art of Computer Programming. — 3-е изд. — М.: Вильямс, 2000. — Т. 2. Получисленные алгоритмы. — 832 с. — 7000 экз. — ISBN 5-8459-0081-6 (рус.) ISBN 0-201-89684-2 (англ.)
я использую
text
Xi+1 = (29 * Xi + 0x7A4537ED)MOD2^32умножение выполняется по модулю 2^64
Это, простите меня, неприемлемо.
почему?
А смысл ещё вообще выкладывать тогда в общий доступ?
для монте-карло - то, что доктор прописал. Если у вас куча свободного времени и сверх-быстрый CPU, вы можете рассчитать md4 от вышеприведённой ПСП.
Или это просто документация устаревшая может быть?
дык, насколько я знаю, с 70х годов прошлого века ничего на этих фронтах не поменялось.
Да, ещё есть генератор цепного кода - сдвигаете влево на 1 бит регистр, причём в младший бит вдвигаете сумму по модулю 2 от старшего, и ещё какого-то бита. (какого? зависит от разрядности N). Если N не степень двойки (точнее не x^2^2 или как-то так, не помню точно), то имеется 2 цикла - один длинной 2^N-1 и второй длинной 1, из числа 0000000...00000. Очень эффективно на аппаратном уровне. К сожалению, если Xi значительно меньше 2^N, то Xi+1 ~= 2*Xi. Но для монте-карло и прочих 99% приложений это не важно.
почему?
да как-то странно %) заявлять о том, что функция неработоспособна... а почему нельзя пофиксить?
Если самое то, говорите - то здорово. Я, кстати, уже проверил - работает... Считает правильно (примерно). Подходит. Описание настораживает. =)
Ну да. Криптостойкость мне не нужна... Наверное, вы правы. =)
да как-то странно %) заявлять о том, что функция неработоспособна... а почему нельзя пофиксить?
а в том-то и дело: пофиксить можно, но будет работать раз так в 10 медленнее. при этом, точность вычисления интегралов не изменится. Спрашивается - зачем?
Если самое то, говорите - то здорово. Я, кстати, уже проверил - работает... Считает правильно (примерно). Подходит. Описание настораживает.
ну такое у них чувство юмора. :lol:
а "точно" - это вам в другие алгоритмы ;) монте-карло вправе выдать вам неверный результат, просто вероятность этого невелика. Но таки есть.
ЗЫЖ вы спрашивали про другие распределения. специальных ГПСЧ для них нет, другие распределения получаются из равномерного с помощью преобразований. формулы не помню, они в справочниках или легко гуглятся.
а в том-то и дело: пофиксить можно, но будет работать раз так в 10 медленнее. при этом, точность вычисления интегралов не изменится. Спрашивается - зачем?
Ах, вот оно что... Вопрос тогда снят.
а "точно" - это вам в другие алгоритмы wink монте-карло вправе выдать вам неверный результат, просто вероятность этого невелика. Но таки есть.
ЗЫЖ вы спрашивали про другие распределения. специальных ГПСЧ для них нет, другие распределения получаются из равномерного с помощью преобразований. формулы не помню, они в справочниках или легко гуглятся.
Да нет... в boost они вроде бы есть. Но получить их действительно можно, особой проблемы не представляет. Да и пока мне не нужно...
А монте-карло мне нужен потому что вроде как его трудоёмкость меньше зависит от размерности, чем в других алгоритмах. А мне это критично... нужно считать очень многомерные интегралы (мне так научник сказал).
Пока посчитал только простой интеграл (одномерный, по отрезку). Вот думаю... а кратные считать как? %) приводить к повторным? Или наверное всё-таки не так, а как-то по-другому? Ведь если приводить к повторным интеграл N-й кратности, то получится, что нужно считать N интегралов... тут и точность падает и трудоёмкость вычислений увеличивается.
Или наверное надо как-то по-другому? %) чтобы СРАЗУ посчитать, а не N раз считать?
Да нет... в boost они вроде бы есть.
ну тем лучше. smile
А монте-карло мне нужен потому что вроде как его трудоёмкость меньше зависит от размерности, чем в других алгоритмах. А мне это критично... нужно считать очень многомерные интегралы (мне так научник сказал).
правильно сказал. ЕМНИП сложность вычисления по методу монте-карло растёт намного медленнее с увеличением числа измерений.
Пока посчитал только простой интеграл (одномерный, по отрезку). Вот думаю... а кратные считать как? %)
в смысле "кратные"?
считают многомерные так-же как одномерные. Разница лишь в том, что при вычислении одномерного интеграла вы забиваете случайными точками отрезок (вычисляете функции f(x) для случайных x), а при вычислении многомерных вам надо забить точками многомерный гиперкуб, т.е. вычислять функции f(x,y,z) для трёхмерного интеграла, и для трёх случайных координат(x,y,z) точки внутри куба. Т.е. для получения одной точки внутри гиперкуба с N измерениями вам достаточно вычислить N случайных координат, что намного проще, чем возится с аппроксимированием гипер-поверхностей (например для одномерного интеграла вы можете легко аппроксимировать функцию кусочками парабол, но для N мерного это на порядки сложнее).
А вот как вычисляются именно кратные - я не помню.
Войдите или зарегистрируйтесь, чтобы ответить.
