Как следует подходить к проблеме проекта Эйлера 213 («Блошиный цирк»)?

11

Я хотел бы решить Project Euler 213, но не знаю, с чего начать, потому что я непрофессионал в области статистики, обратите внимание, что требуется точный ответ, чтобы метод Монте-Карло не работал. Не могли бы вы порекомендовать мне некоторые темы статистики? Пожалуйста, не размещайте решение здесь.

Блошиный цирк

Сетка квадратов 30 × 30 содержит 900 блох, первоначально по одной блохе на квадрат. Когда звонит колокол, каждая блоха прыгает на соседний квадрат случайным образом (обычно 4 варианта, кроме блох на краю сетки или в углах).

Каково ожидаемое количество незанятых квадратов после 50 звонков колокола? Дайте ответ, округленный до шести знаков после запятой.

grokus
источник
7
Методы Монте-Карло могут дать очень точные ответы при условии, что вы выполняете достаточно моделирования.
Роб Хиндман
3
Если вам нужно программное решение, Монте-Карло - единственный подход. Я не вижу никаких причин, почему вы не получите точные ответы, используя Монте-Карло. Математическое / аналитическое решение может быть нелегким.
Я видел дискуссию о Монте-Карло, и люди говорили, что если вы хотите получить 6 знаков после запятой, это займет слишком много времени, или, возможно, меня смущают другие подобные проблемы. Поскольку довольно легко закодировать подход Монте-Карло, я думаю, стоит сначала попробовать.
Грокус
4
Я не оспариваю ни один из трех предыдущих ответов, но (простой) анализ в ответе, который я предложил, ставит эти замечания в перспективе: если вам нужна точность в шесть знаков после запятой для оценки числа, которое будет исчисляться сотнями, Симуляция Монте-Карло займет не менее года на машине с 10 000 процессоров, работающих параллельно.
whuber
Все ли блохи пойманы в ловушку (то есть проблема действительно в квадратах с более чем одной блохой на них) или в том, что блохи по краям выпрыгивают и исчезают?
MissMonicaE

Ответы:

10

Вы правы; Монте-Карло невыполнимо. (В наивном моделировании, то есть в том, которое точно воспроизводит проблемную ситуацию без каких-либо упрощений, каждая итерация будет включать в себя 900 ходов блох. Грубая оценка доли пустых ячеек равна , что подразумевает дисперсию Монте -Карловая оценка после N таких итераций составляет примерно 1 / N 1 / e ( 1 - 1 / e ) = 0,2325 / N1/еN1/N1/е(1-1/е)знак равно0,2325.../N, Чтобы закрепить ответ до шести десятичных разрядов, вам нужно будет оценить его с точностью до 5.E-7, и, чтобы получить уверенность в 95 +% (скажем), вам придется примерно вдвое уменьшить эту точность до 2.5E-7 , Решая даетN>4E12, приблизительно. Это будет около 3,6E15 блох, каждый из которых занимает несколько тактов процессора. При наличии одного современного процессора вам понадобится целый год (высокоэффективных) вычислений. И я несколько неправильно и чрезмерно оптимистично предположил, что ответ дается в виде пропорции, а не подсчета: для подсчета ему понадобятся еще три значащие цифры, что повлечет за собой увеличение вычислений в миллион раз ... Вы можете ждать долго?)(0,2325/N)<2.5Е-7N>4Е12

Что касается аналитического решения, существуют некоторые упрощения. (Они также могут использоваться для сокращения вычислений по методу Монте-Карло.) Ожидаемое количество пустых ячеек - это сумма вероятностей пустоты по всем ячейкам. Чтобы найти это, вы можете вычислить распределение вероятностей чисел занятости каждой ячейки. Эти распределения получены путем суммирования (независимых!) Вкладов от каждой блохи. Это сводит вашу проблему к нахождению количества путей длиной 50 по сетке 30 на 30 между любой данной парой ячеек в этой сетке (одна - источник блох, а другая - ячейка, для которой вы хотите вычислить вероятность размещение блох).

Whuber
источник
2
Просто для забавы, я сделал расчет грубой силы в Mathematica. Его ответ - это отношение целого числа 21 574 к целому числу 21 571; как десятичное число это удобно близко к 900 / e, как и ожидалось (но, поскольку нас просят не публиковать решение, я не буду давать больше подробностей).
whuber
6

Не могли бы вы перебрать вероятности занятия клеток для каждой блохи. То есть блоха k изначально находится в ячейке (i (k), j (k)) с вероятностью 1. После 1 итерации он имеет вероятность 1/4 в каждой из 4 соседних ячеек (при условии, что он не на краю или в угол). Затем на следующей итерации каждый из этих кварталов по очереди «размазывается». После 50 итераций у вас есть матрица вероятностей занятия для блох k. Повторите все 900 блох (если вы используете преимущества симметрии, это уменьшает почти в 8 раз) и добавьте вероятности (вам не нужно хранить их все сразу, только матрицу текущей блохи (хмм, если вы не очень умный, вам может понадобиться дополнительная рабочая матрица) и текущая сумма матриц). Мне кажется, что есть много способов ускорить это здесь и там.

Это не требует никакой симуляции вообще. Тем не менее, это довольно много вычислений; не должно быть очень сложно определить размер симуляции, необходимый для того, чтобы дать ответы с несколько большей точностью, чем 6 dp с высокой вероятностью, и выяснить, какой подход будет быстрее. Я ожидаю, что этот подход побьет симуляцию с некоторым запасом.

Glen_b - Восстановить Монику
источник
2
Ваш ответ немного отличается от вопроса, который задает вопрос. Вопрос в том, чтобы узнать ожидаемое количество ячеек, которые будут пустыми после 50 прыжков. Поправьте меня, если я ошибаюсь, но я не вижу прямого пути от вероятности, что блоха окажется в определенном квадрате после 50 прыжков к ответу, сколько ячеек должно быть пустым.
Энди W
1
@Andy W - отличный комментарий; тем не менее, Монте-Карло можно использовать для этого последнего шага ;-)
4
@ Andy W: На самом деле, самая сложная часть - получить все эти вероятности. Вместо добавления их в каждую ячейку, умножьте их дополнения: это вероятность того, что ячейка будет пустой. Сумма этих значений по всем ячейкам дает ответ. Подход Glen_b превосходит симуляцию на семь или восемь порядков ;-).
whuber
@whuber, спасибо за объяснение. Действительно, получить эти вероятности менее чем за минуту было бы сложно. Это веселая головоломка и спасибо за ваш вклад.
Энди W
5

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

Во-первых, после Glen_b частицы могут обмениваться в стационарном режиме, поэтому достаточно (как при достаточности ) контролировать занятость разных ячеек, поскольку это также составляет марковский процесс. Распределение занятых мест на следующем шаге времени завершено, определяется занятостью в текущий момент времени t . Написание матрицы переходов K определенно нецелесообразно, но моделировать переход просто.T+1TК

Во-вторых, как отметил Шаббычеф , можно следить за процессом заполнения на 450 нечетных (или четных) квадратах, который остается на нечетных квадратах при рассмотрении только четных времен, т. Е. Квадратовой матрицы Маркова .К2

В- третьих, исходная задача рассматривает только частоту нулевых , после 50 марковских переходов. Принимая во внимание , что начальная точка имеет очень высокое значение для стационарного распределения вероятностей цепи Маркова ( Х ( т ) ) , и при условии , что фокус на одном среднем по всем клеткам, р 0 = 1п^050(Икс(T))можно считать, что реализация цепочки(X(t))в момент времениt=50является реализацией из стационарного распределения вероятностей. Это значительно снижает стоимость вычислений, поскольку мы можем смоделировать непосредственно из этого стационарного распределенияπ, которое является многочленным распределением с вероятностями, пропорциональными 2, 3 и 4 на четном угле, другими ячейками на краю и внутренними ячейками. соответственно.

п^0знак равно1450Σязнак равно1450я0(Икся(50))
(Икс(T))Tзнак равно50π

Σязнак равно1450(1-πя)450
166.1069
pot=rep(c(rep(c(0,1),15),rep(c(1,0),15)),15)*c(2,
    rep(3,28),2,rep(c(3,rep(4,28),3),28),2,rep(3,28),2)
pot=pot/sum(pot)
sum((1-pot)^450)-450
[1] 166.1069

166,11

Как прокомментировал whuber , оценки должны быть умножены на 2, чтобы правильно ответить на вопрос, следовательно, окончательное значение 332,2137,

Сиань
источник
1
+1 Очень проницательно. Я считаю, что вам нужно удвоить свой окончательный ответ, потому что вопрос задает обо всех 900 ячейках.
whuber
1
Я полагаю, что вы, возможно, начинаете дальше от стационарного распределения, чем вы думаете. Расчеты методом грубой силы, которые я первоначально делал, вычисляли 50-ую степень матрицы перехода, используя точную (рациональную) арифметику. Из него я получил значение 330.4725035083710 .... Возможно, я допустил ошибку .... У меня действительно была ошибка, и теперь я получаю 330.7211540144080 .... Обширная проверка показывает, что матрица перехода верна.
whuber
@whuber: Спасибо, это действительно возможно. Я попытался найти аргумент связи, чтобы определить скорость стационарности, но не смог. Симуляция Монте-Карло с оригинальным процессом дала мне 333,96 на 10⁶ реплик и 57 часов вычислений. Без дополнительной гарантии на точность.
Сиань
1
Вот мои рассуждения. Матрица перехода для 50 шагов является 50-й степенью матрицы перехода, откуда ее собственные значения являются 50-й степенью собственных значений. Только собственные векторы, соответствующие значениям, чьи 50-ые степени имеют любой заметный размер, будут отображаться в качестве компонентов в конце ваших 50 шагов. Более того, эти 50-е полномочия информируют нас об относительной ошибке, допущенной при остановке на 50-м шаге, а не по-настоящему достижении устойчивого состояния.
whuber
1
900×900
4

Аналитический подход может быть утомительным, и я не продумал все тонкости, но вот подход, который вы можете рассмотреть. Так как вас интересует ожидаемое количество ячеек, которые будут пустыми после 50 колец, вам нужно определить цепочку Маркова над «Нет блох в клетке», а не положение блох (см. Ответ Glen_b, который моделирует положение блоха как цепочка Маркова. Как отметил Энди в комментариях к этому ответу, такой подход может не дать того, что вы хотите.)

В частности, пусть:

nij(t)яJ

Затем цепочка Маркова начинается со следующего состояния:

niJ(0)знак равно1яJ

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

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

Сообщество
источник
1
NяJ
@whuber Нет, вам не нужно поддерживать блох в качестве цепочки Маркова. Думайте о том, что я предлагаю как случайную прогулку по клетке. Первоначально ячейка находится в положении «1», откуда она может перейти в 0, 1, 2, 3, 4 или 5. Вероятность перехода состояния зависит от состояний соседних ячеек. Таким образом, предлагаемая цепь находится в переопределенном пространстве состояний (количество клеток для каждой клетки), а не в самой позиции блох. Имеет ли это смысл?
1
Это имеет смысл, но это похоже на шаг назад, потому что теперь число штатов не намного больше? В одной модели 900 состояний - положение одной блохи - и не более четырех переходов из каждого. Вычисление должно быть сделано только для одной блохи, потому что они все двигаются независимо. По-вашему, состояние описывается занятостью ячейки и занятостью ее до четырех соседей. Это будет чрезвычайно большое количество штатов, а также очень большое количество переходов между штатами. Должно быть, я неправильно понимаю, каково ваше новое пространство состояний.
whuber
{NяJ}
2

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

shabbychef
источник
1
Это хорошее наблюдение. Тем не менее, я нашел, что это беспокоит больше, чем стоит явно использовать это. Большая часть программирования сводится к настройке матрицы перехода. Как только вы это сделаете, просто возведите в квадрат и работайте с этим. При использовании разреженных матриц удаление половины нулей не экономит время в любом случае.
whuber
@whuber: Я подозреваю, что смысл этих проблем в том, чтобы изучить методы решения проблем, а не потреблять много вычислительных циклов. Симметрия, четность и т. Д. - это классические приемы из книги Ларсона по решению проблем.
Шаббычеф
1
Неплохо подмечено. В конечном счете, необходимо некоторое суждение. Проект Euler, кажется, подчеркивает компромиссы между математическим пониманием и вычислительной эффективностью. Glen_b упомянул симметрии, которые стоит использовать в первую очередь, потому что от них можно извлечь еще больше. Более того, используя разреженную матричную арифметику, вы автоматически получите двукратное усиление (знаете ли вы о четности или нет!).
whuber
1

Я подозреваю, что некоторые знания о цепях Маркова с дискретным временем могут оказаться полезными.

Саймон Бирн
источник
3
Это должен был быть комментарий, но я думаю, что мы можем дедушку в этом месте.
gung - Восстановить Монику
Это автоматически помечается как низкое качество, возможно потому, что оно очень короткое. Вы можете расширить это?
gung - Восстановить Монику
Я не понимаю, почему: вопрос требует темы, которая может быть полезной, и эта тема, на мой взгляд, наиболее актуальна.
Саймон Бирн
1
Это было отмечено как низкое качество . Я проголосовал, что все в порядке. Если вы посмотрите на другие ответы в этой теме, они все значительно длиннее. Стандарты со временем развивались, но сегодня это будет рассматриваться как комментарий, даже если упоминается «тема, которая может быть полезной». Как я уже сказал, я думал, что это может быть сделано как есть. Пытаетесь ли вы расширить его, зависит от вас. Я просто дал вам знать.
gung - Восстановить Монику