Недавно я смотрел на симуляцию Монте-Карло и использовал ее для аппроксимации констант, таких как (окружность внутри прямоугольника, пропорциональная область).
Однако я не могу придумать соответствующий метод аппроксимации значения [число Эйлера] с использованием интеграции Монте-Карло.
Есть ли у вас какие-либо указания о том, как это можно сделать?
simulation
monte-carlo
algorithms
random-generation
numerical-integration
статистикаnewbie12345
источник
источник
R
команда2 + mean(exp(-lgamma(ceiling(1/runif(1e5))-1)))
. (Если использование гамма-функции журнала вас беспокоит, замените ее2 + mean(1/factorial(ceiling(1/runif(1e5))-2))
, которая использует только сложение, умножение, деление и усечение, и игнорируйте предупреждения о переполнении.) Что может быть более интересным, так это эффективные симуляции: можете ли вы минимизировать количество вычислительные шаги, необходимые для оценкиОтветы:
Простой и элегантный способ оценки по Монте-Карло описан в этой статье . Бумага на самом деле об обучении e . Следовательно, подход кажется идеально подходящим для вашей цели. Идея основана на упражнении из популярного российского учебника по теории вероятностей Гнеденко. Смотри ex.22 на стр.183e e
Случается так, что , где ξ - случайная величина, которая определяется следующим образом. Это минимальное число п такое , что Σ п я = 1 г я > 1 и г я случайные числа от равномерного распределения на [ 0 , 1 ] . Красиво, не правда ли ?!E[ξ]=e ξ n ∑ni=1ri>1 ri [0,1]
Поскольку это упражнение, я не уверен, что для меня было бы здорово опубликовать решение (доказательство) здесь :) Если вы хотите доказать это сами, вот совет: глава называется «Моменты», которая должна указывать Вы в правильном направлении.
Если вы хотите реализовать это самостоятельно, не читайте дальше!
Это простой алгоритм для моделирования Монте-Карло. Нарисуйте случайный случайный случай, затем еще один и так далее, пока сумма не превысит 1. Количество нарисованных рандомов - ваше первое испытание. Допустим, вы получили:
Затем ваше первое испытание было выполнено 3. Продолжайте выполнять эти испытания, и вы заметите, что в среднем вы получаете .e
Код MATLAB, результат моделирования и гистограмма следуют.
Результат и гистограмма:
ОБНОВЛЕНИЕ: я обновил свой код, чтобы избавиться от массива результатов испытаний, чтобы он не занимал ОЗУ. Я также напечатал оценку PMF.
Обновление 2: вот мое решение Excel. Поместите кнопку в Excel и свяжите ее со следующим макросом VBA:
Введите количество испытаний, например 1000, в ячейку D1 и нажмите кнопку. Вот как должен выглядеть экран после первого запуска:
ОБНОВЛЕНИЕ 3: Серебряная рыбка вдохновила меня на другой путь, не такой элегантный, как первый, но все же крутой. Он рассчитал объемы n-симплексов с использованием последовательностей Соболя .
По совпадению он написал первую книгу о методе Монте-Карло, которую я прочитал еще в старшей школе. На мой взгляд, это лучшее введение в метод.
ОБНОВЛЕНИЕ 4:
Silverfish в комментариях предложил простую реализацию формулы Excel. Вот такой результат вы получите с его подходом после примерно 1 миллиона случайных чисел и 185K испытаний:
Очевидно, что это намного медленнее, чем реализация Excel VBA. Особенно, если вы измените мой код VBA, чтобы не обновлять значения ячеек внутри цикла, а делать это только после сбора всей статистики.
ОБНОВЛЕНИЕ 5
Решение Сианя №3 тесно связано (или даже в некотором смысле то же, что и комментарий jwg в ветке). Сложно сказать, кто придумал эту идею сначала Форсайт или Гнеденко. Оригинальное издание Гнеденко 1950 года на русском языке не имеет разделов «Проблемы» в главах. Таким образом, я не смог найти эту проблему с первого взгляда, где она есть в более поздних выпусках. Возможно это было добавлено позже или похоронено в тексте.
Как я прокомментировал в ответе Сианя, подход Форсайта связан с другой интересной областью: распределением расстояний между пиками (экстремумами) в случайных (IID) последовательностях. Среднее расстояние оказывается равным 3. Последовательность обратного хода в подходе Форсайта заканчивается дном, поэтому, если вы продолжите выборку, вы получите еще одно дно в какой-то момент, затем еще одно и т. Д. Вы можете отслеживать расстояние между ними и строить распределение.
источник
Mean[Table[ Length[NestWhileList[(Random[]+#) &, Random[], #<1&]], {10^6}]]
R
решения, которое я разместил в ответе Сианя, в двадцать раз быстрее:n=10^6; 1. / Mean[UnitStep[Differences[Sort[RandomReal[{0, n}, n + 1]]] - 1]]
Я предлагаю опровергнуть ответ Аксакала. Он беспристрастен и опирается только на метод генерации единичных отклонений.
поэтому мы также можем написать
То есть наша оценка находится путем оценки вероятности того, что конкретное наблюдение опущено из начальной загрузки во многих таких - то есть доли вхождений объекта в начальной загрузке.p m Bj i
В этом приближении есть два источника ошибок. Конечное всегда будет означать, что результаты приблизительны, то есть оценка смещена. Кроме того, будет колебаться вокруг истинного значения, потому что это симуляция.n p^
Я нахожу этот подход несколько очаровательным, потому что студент или другой человек, у которого достаточно мало дел, могут приблизиться к используя колоду карт, кучу маленьких камней или любые другие предметы в руке, в том же ключе, что и человек, которого можно оценить используя компас, прямой край и несколько песчинок. Я думаю, что это хорошо, когда математика может быть отделена от современных удобств, таких как компьютеры.e π
Полученные результаты
Я провел несколько симуляций для различного числа загрузочных повторов. Стандартные ошибки оцениваются с использованием нормальных интервалов.
Обратите внимание, что выбор числа загружаемых объектов устанавливает абсолютный верхний предел точности результатов, поскольку процедура Монте-Карло оценивает а зависит только от . Установка излишне большим просто приведет к перегружению вашего компьютера, либо потому, что вам нужно только «грубое» приближение к либо потому, что смещение будет затоплено дисперсией из-за Монте-Карло. Эти результаты для и точностью до третьего знака после запятой.n p p n n e n=103 p−1≈e
Этот график показывает, что выбор имеет прямые и глубокие последствия для стабильности в . Синяя пунктирная линия показывает а красная линия показывает . Как и ожидалось, увеличение размера выборки дает все более точные оценки .m p^ p e p^
Я написал смущающе длинный R-скрипт для этого. Предложения по улучшению могут быть представлены на обратной стороне 20-долларовой купюры.
источник
Решение 1:
Для распределения Пуассона Следовательно, если , что означает, что вы можете оценить с помощью моделирования Пуассона. И моделирование Пуассона может быть получено из генератора экспоненциального распределения (если не самым эффективным способом).P(λ)
Решение 2:
Другой способ добиться представления константы в виде интеграла состоит в том, чтобы напомнить, что когда то который также является распределением. Таким образом, Второй подход к аппроксимации от Монте-Карло, таким образом, должен моделировать нормальные пары и контролировать частоту времен . В некотором смысле это противоположность аппроксимации Монте-Карло для связанной с частотой времен ...e
Решение 3:
Мой коллега из Уорикского университета М. Поллок указал на другое приближение Монте-Карло, называемое методом Форсайта : идея состоит в том, чтобы запустить последовательность равномерных поколений до . Ожидание соответствующего правила остановки , которое представляет собой число раз, когда равномерная последовательность пошла вниз, равно тогда как вероятность того, что нечетно, равна ! ( Метод Форсайта на самом деле нацелен на моделирование из любой плотности вида , поэтому он более общий, чем аппроксимация и .)u1,u2,... un+1>un N e N e−1 expG(x) e e−1
Быстрая R-реализация метода Форсайта состоит в том, чтобы отказаться от точного следования последовательности униформ в пользу больших блоков, что позволяет выполнять параллельную обработку:
источник
n <- 1e5; 1/mean(n*diff(sort(runif(n+1))) > 1)
Не решение ... просто быстрый комментарий, который слишком длинный для поля комментария.
Аксакал
Аксакал опубликовал решение, в котором мы рассчитываем ожидаемое количество стандартных чертежей Униформы, которые должны быть взяты так, чтобы их сумма превышала 1. В Mathematica моей первой формулировкой было:
РЕДАКТИРОВАТЬ: только что быстро поиграть с этим, и следующий код (тот же метод - также в Mma - просто другой код) примерно в 10 раз быстрее:
Сиань / Вубер
Whuber предложил быстрый классный код для имитации решения Сианя 1:
Версия R:
n <- 1e5; 1/mean(n*diff(sort(runif(n+1))) > 1)
ММА версия:
n=10^6; 1. / Mean[UnitStep[Differences[Sort[RandomReal[{0, n}, n + 1]]] - 1]]
который он отмечает, это в 20 раз быстрее первого кода (или примерно в два раза быстрее, чем новый код выше).
Просто ради интереса, я подумал, что было бы интересно посмотреть, насколько оба подхода эффективны (в статистическом смысле). Для этого я сгенерировал 2000 оценок e, используя:
... оба в Mathematica . Следующая диаграмма сравнивает непараметрическую оценку плотности ядра результирующих наборов данных dataA и dataB.
Таким образом, хотя код whuber (красная кривая) примерно в два раза быстрее, метод не выглядит «надежным».
источник
running four times as many iterations will make them equally accurate
///// ..... Просто поигрался с этим: увеличил количество точек выборки, использованных в методе 1 Сианя, с до 6 x (то есть в 6 раз больше числа баллы) производит кривую, аналогичную Аксаксал.Метод, требующий нечестного количества образцов
Во-первых, вы должны иметь возможность выбирать из нормального распределения. Предполагая, что вы собираетесь исключить использование функции или просмотреть таблицы, полученные из этой функции, вы можете получить приблизительные выборки из нормального распределения через CLT. Например, если вы можете сделать выборку из равномерного (0,1) распределения, то . Как указал Уубер, чтобы подход окончательной оценки приближался к размеру выборки , необходимо, чтобы количество используемых однородных выборок приближалось к когда размер выборки приближался к бесконечности.f(x)=ex x¯12√n√∼˙N(0,1) e ∞ ∞
Теперь, если вы можете производить выборку из нормального распределения с достаточно большими выборками, вы можете получить непротиворечивые оценки плотности . Это можно сделать с помощью гистограмм или сглаживателей ядра (но будьте осторожны, чтобы не использовать ядро Гаусса для следования правилу no !). Чтобы ваши оценки плотности были согласованными, вам нужно, чтобы df (число бинов на гистограмме, обратное к окну для сглаживания) приближался к бесконечности, но медленнее, чем размер выборки.N(0,1) ex
Так что теперь, с большой вычислительной мощностью, вы можете аппроксимировать плотность , то есть . Поскольку , ваша оценка для .N(0,1) ϕ^(x) ϕ((√2))=(2π)−1/2e−1 e=ϕ^(2–√)2π−−√
Если вы хотите полностью сходить с ума, вы можете даже оценить и используя методы, которые вы обсуждали ранее.2–√ 2π−−√
Метод, требующий очень мало выборок, но вызывающий неоправданное количество числовых ошибок
Совершенно глупый, но очень эффективный ответ, основанный на комментарии, который я сделал:
Пусть . Определить, Определить .X∼uniform(−1,1) Yn=|(x¯)n| e^=(1−Yn)−1/Yn
Это будет сходиться очень быстро, но также приведет к крайней числовой ошибке.
Уубер указал, что здесь используется функция power, которая обычно вызывает функцию exp. Это можно обойти путем дискретизации , так что является целым числом, а мощность можно заменить повторным умножением. Было бы необходимо, чтобы при дискретизация становилась все точнее, а дискретизация должна исключать . При этом теоретически оценщик (т. Е. Мир, в котором не существует числовой ошибки) будет сходиться к и довольно быстро!Yn 1/Yn n→∞ Yn Yn=0 e
источник
Вот еще один способ сделать это, хотя это довольно медленно. Я не претендую на эффективность, но предлагаю эту альтернативу в духе полноты.
В отличие от ответа Сианя , для целей этого вопроса я предполагаю, что вы можете сгенерировать и использовать последовательность из равномерных псевдослучайных переменных и затем вам нужно оценить некоторым методом с использованием основных арифметических операций (то есть вы не можете использовать логарифмические или экспоненциальные функции или любые распределения, использующие эти функции). Данный метод мотивирован простым результатом, включающим однородные случайные величины:n U1,⋯,Un∼IID U(0,1) e †
Оценка с использованием этого результата:e сначала мы упорядочиваем значения выборки в порядке убывания, чтобы получить статистику порядка а затем определяем частичные суммы:u(1)⩾⋯⩾u(n)
Теперь пусть а затем оцените путем интерполяции упорядоченных равномерных переменных. Это дает оценку для :m≡min{k|S(k)⩾1} 1/e e
Этот метод имеет небольшое смещение (из-за линейной интерполяции точки отсечки для ), но он является последовательной оценкой для . Метод может быть реализован довольно легко, но он требует сортировки значений, что требует больших вычислительных затрат, чем детерминистский расчет . Этот метод медленный, так как включает в себя сортировку значений.1/e e e
Реализация в R: метод может быть реализован в
R
использованииrunif
для генерации однородных значений. Код выглядит следующим образом:Реализация этого кода дает сходимость к истинному значению , но она очень медленная по сравнению с детерминированными методами.e
источник