Jak należy podejść do problemu Euler Projektu 213 („Flea Circus”)?

11

Chciałbym rozwiązać Project Euler 213, ale nie wiem od czego zacząć, ponieważ jestem laikiem w dziedzinie statystyki, zauważ, że wymagana jest dokładna odpowiedź, aby metoda Monte Carlo nie zadziałała. Czy mógłbyś polecić mi kilka tematów statystycznych do przeczytania? Proszę nie zamieszczać rozwiązania tutaj.

Flea Circus

Siatka kwadratów 30 × 30 zawiera 900 pcheł, początkowo jedną pchłę na kwadrat. Kiedy dzwonek jest dzwoniony, każda pchła losowo skacze na sąsiedni kwadrat (zwykle 4 możliwości, z wyjątkiem pcheł na krawędzi siatki lub w rogach).

Jaka jest oczekiwana liczba niezajętych kwadratów po 50 dzwonkach dzwonu? Podaj odpowiedź w zaokrągleniu do sześciu miejsc po przecinku.

self-study monte-carlo markov-process

— grokus
źródło

7

Metody Monte Carlo mogą dać bardzo dokładne odpowiedzi, pod warunkiem, że wykonasz wystarczającą liczbę symulacji.

— Rob Hyndman,

3

Jeśli szukasz rozwiązania programistycznego, jedynym rozwiązaniem jest Monte Carlo. Nie widzę żadnego powodu, dla którego nie uzyskałbyś dokładnych odpowiedzi przy użyciu Monte Carlo. Rozwiązanie matematyczne / analityczne może nie być łatwe.

Widziałem dyskusję na temat Monte Carlo i ludzie mówili, że jeśli chcesz osiągnąć 6 miejsc po przecinku, potrwa to zbyt długo lub być może mylę się z innymi podobnymi problemami. Ponieważ dość łatwo jest zakodować podejście Monte Carlo, myślę, że warto najpierw spróbować.

— grokus

4

Nie kwestionuję żadnej z trzech poprzednich odpowiedzi, ale (prosta) analiza w odpowiedzi, którą zaoferowałem, przedstawia te uwagi w perspektywie: jeśli chcesz dokładności z dokładnością do sześciu miejsc po przecinku dla oszacowania liczby, która będzie w setkach, symulacja Monte Carlo potrwa co najmniej rok na maszynie z 10 000 równolegle działających procesorów.

— whuber

Czy wszystkie pchły są uwięzione (tzn. Problem tak naprawdę dotyczy kwadratów z więcej niż jedną pchłą), czy też chodzi o pchły na krawędziach wyskakujące i znikające?

— MissMonicaE

10

Masz rację; Monte Carlo jest niewykonalne. (W naiwnej symulacji - tzn. Takiej, która dokładnie odtwarza sytuację problemową bez żadnych uproszczeń - każda iteracja wymagałaby 900 ruchów pcheł. Szacunkowy procent proporcji pustych komórek wynosi , co sugeruje wariancję Monte - Oszacowanie Carlo po takich iteracjach wynosi około Aby określić odpowiedź z do sześciu miejsc po przecinku, należy ją oszacować z dokładnością do 5.E -7 i, aby osiągnąć poziom ufności 95 +% (powiedzmy), musiałbyś w przybliżeniu zmniejszyć o połowę tę precyzję do 2,5E-7. Rozwiązanie daje $1/e$ $N$ $1/N 1/e (1 - 1/e) = 0.2325\ldots /N$ $\sqrt(0.2325/N) \lt 2.5E-7$ $N > 4E12$ około To byłoby około 3,6E15 ruchów pcheł, z których każdy wymagałby kilku kliknięć procesora. Z jednym dostępnym nowoczesnym procesorem będziesz potrzebował całego roku (wysoce wydajnego) przetwarzania. I nieco niepoprawnie i nadmiernie optymistycznie założyłem, że odpowiedź jest podana w postaci proporcji zamiast liczby: jako liczba będzie potrzebowała jeszcze trzech znaczących liczb, co pociągnie za sobą milion-krotny wzrost obliczeń ... Czy możesz długo czekać?)

Jeśli chodzi o rozwiązanie analityczne, dostępne są pewne uproszczenia. (Można ich również użyć do skrócenia obliczeń Monte Carlo). Oczekiwana liczba pustych komórek jest sumą prawdopodobieństw pustki we wszystkich komórkach. Aby to znaleźć, możesz obliczyć rozkład prawdopodobieństwa liczby zajętości każdej komórki. Rozkłady te uzyskuje się poprzez zsumowanie (niezależnego!) Wkładu każdej pchły. Zmniejsza to problem ze znalezieniem liczby ścieżek o długości 50 wzdłuż siatki 30 na 30 między dowolną parą komórek na tej siatce (jedna jest początkiem pcheł, a druga komórką, dla której chcesz obliczyć prawdopodobieństwo obłożenie pcheł).

— Whuber
źródło

2

Dla zabawy wykonałem obliczenia z użyciem siły brutalnej w Mathematica. Jego odpowiedzią jest stosunek liczby całkowitej 21 574 cyfr do liczby całkowitej 21 571 cyfr; w systemie dziesiętnym jest wygodnie zbliżony do 900 / e zgodnie z oczekiwaniami (ale ponieważ jesteśmy proszeni o nie publikowanie rozwiązania, nie podam więcej szczegółów).

— whuber

6

Czy nie możesz iterować po prawdopodobieństwie zajęcia komórek dla każdej pcheł. Oznacza to, że pchła k jest początkowo w komórce (i (k), j (k)) z prawdopodobieństwem 1. Po 1 iteracji ma prawdopodobieństwo 1/4 w każdej z 4 sąsiednich komórek (zakładając, że nie jest na krawędzi ani w róg). Następnie, w następnej iteracji, każda z tych ćwiartek zostaje „rozmazana” z kolei. Po 50 iteracjach masz macierz prawdopodobieństw okupacji dla pcheł k. Powtórz ponad 900 pcheł (jeśli skorzystasz z symetrii, zmniejszy to prawie o współczynnik 8) i dodaj prawdopodobieństwa (nie musisz przechowywać ich wszystkich naraz, tylko macierz bieżącej pchły (hmm, chyba że jesteś bardzo sprytne, możesz potrzebować dodatkowej działającej macierzy) i bieżącej sumy macierzy). Wydaje mi się, że istnieje wiele sposobów na przyspieszenie tego tu i tam.

Nie wymaga to żadnej symulacji. Wymaga to jednak sporo obliczeń; nie powinno być bardzo trudno obliczyć rozmiar symulacji wymagany do uzyskania odpowiedzi z nieco lepszą dokładnością niż 6 dp z dużym prawdopodobieństwem i dowiedzieć się, które podejście będzie szybsze. Oczekuję, że takie podejście przewyższy symulację o pewien margines.

— Glen_b - Przywróć Monikę
źródło

2

Twoja odpowiedź na nieco inne pytanie niż pytanie. Pytanie brzmi: spodziewana liczba komórek, które byłyby puste po 50 skokach. Popraw mnie, jeśli się mylę, ale nie widzę bezpośredniej ścieżki z prawdopodobieństwa, że pchła trafi na określony kwadrat po 50 skokach do odpowiedzi, ile komórek miałoby być pustych.

— Andy W

1

@Andy W - świetny komentarz; jeszcze Monte Carlo może być użyte do zrobienia tego ostatniego kroku ;-)

4

@Andy W: Tak naprawdę trudność polegała na uzyskaniu wszystkich tych prawdopodobieństw. Zamiast dodawać je do każdej komórki, pomnóż ich uzupełnienia: to prawdopodobieństwo, że komórka będzie pusta. Suma tych wartości we wszystkich komórkach daje odpowiedź. Podejście Glen_b bije symulację o siedem lub osiem rzędów wielkości ;-).

— whuber

@ whuber, dziękuję za wyjaśnienie. Rzeczywiście uzyskanie tych prawdopodobieństw w niecałą minutę byłoby trudne. To zabawna łamigłówka i dziękuję za Twój wkład.

— Andy W

5

Chociaż nie sprzeciwiam się praktycznej niemożności (lub niepraktyczności) rozwiązania tego problemu z Monte Carlo z dokładnością do 6 miejsc po przecinku wskazanej przez whubera , sądzę , że można uzyskać rozdzielczość z sześciocyfrową dokładnością.

Po pierwsze, zgodnie z Glen_b , cząstki są wymienialne w trybie stacjonarnym, a zatem wystarczające (jak w wystarczającym stopniu ) jest monitorowanie zajętości różnych komórek, ponieważ stanowi to również proces Markowa. Rozkład zajętości w następnym kroku jest zakończony, określony przez zajętości w bieżącym czasie . Napisanie macierzy przejścia jest zdecydowanie niepraktyczne, ale symulacja przejścia jest prosta. $t+1$ $t$ $K$

Po drugie, jak zauważył shabbychef można śledzić proces użytkowanie na 450 nieparzyste (parzyste) lub kwadratów, która pozostaje na nieparzystych kwadratów kiedy tylko rozważa nawet razy, czyli kwadratu Markowa macierzy . $K^2$

Po trzecie, oryginalny problem uważa tylko częstotliwości zerowej po przejściach Markowa. Uwzględniając fakt, że punkt początkowy ma bardzo duże znaczenie dla stacjonarnego rozkładu prawdopodobieństwa łańcucha Markowa oraz fakt, że koncentrują się na pojedynczej średnią dla wszystkich $\hat{p}_0$ $50$ $(\mathbf{X}^{(t)})$ możemy uznać, że realizacja łańcuchaw czasiejest realizacją ze stacjonarnego rozkładu prawdopodobieństwa. Zapewnia to znaczną redukcję kosztów obliczeniowych, ponieważ możemy symulować bezpośrednio z tego rozkładu stacjonarnego, który jest rozkładem wielomianowym z prawdopodobieństwami proporcjonalnymi do 2, 3 i 4 w parzystym narożniku, innymi komórkami na krawędzi i komórkami wewnętrznymi odpowiednio.

{\hat{p}}_{0} = \frac{1}{450} \sum_{i = 1}^{450} I_{0} (X_{i}^{(50)})

$\hat{p}_0=\frac{1}{450}\sum_{i=1}^{450}\mathbb{I}_0(X_i^ {(50)})$

(X^{(t)})

$(\mathbf{X}^{(t)})$

t = 50

$t=50$

π

$\mathbf{\pi}$

\sum_{i = 1}^{450} (1 - π_{i})^{450}

$\sum_{i=1}^{450} (1-\pi_i)^{450}$

166.1069

$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$

Jak skomentował Whuber , szacunki należy pomnożyć przez 2, aby poprawnie odpowiedzieć na pytanie, stąd ostateczna wartość 332,2137,

— Xi'an
źródło

1

+1 Bardzo wnikliwe. Uważam, że musisz podwoić ostateczną odpowiedź, ponieważ pytanie dotyczy wszystkich 900 komórek.

— whuber

1

Myślę, że możesz zacząć od dystrybucji stacjonarnej, niż myślisz. Obliczenia siły brutalnej, które pierwotnie wykonałem, obliczyły 50. moc macierzy przejściowej przy użyciu dokładnej (racjonalnej) arytmetyki. Z niego uzyskałem wartość 330.4725035083710 ... Być może popełniłem błąd .... Miałem błąd, a teraz otrzymuję 330.7211540144080 .... Szczegółowe testy sugerują, że macierz przejścia jest poprawna.

— whuber

@whuber: Dzięki, to rzeczywiście możliwość. Próbowałem znaleźć argument łączący w celu ustalenia prędkości do stacjonarności, ale nie mogłem. Symulacja Monte Carlo z oryginalnym procesem dała mi 333,96 ponad 10⁶ replik i 57 godzin obliczeń. Bez dalszej gwarancji precyzji.

— Xi'an

1

Oto moje rozumowanie. Macierz przejściowa dla 50 kroków jest 50. mocą macierzy przejściowej, skąd jej wartości własne są 50. mocami wartości własnych. Tylko wektory własne odpowiadające wartościom, których 50 moce mają znaczny rozmiar, pojawią się jako składniki na końcu twoich 50 kroków. Co więcej, te 50 moce informują nas o względnym błędzie popełnionym przez zatrzymanie się na 50 stopniu zamiast prawdziwego osiągnięcia ustalonego stanu.

— whuber

1

900 \times 900

$900\times 900$

4

Podejście analityczne może być nużące i nie zastanawiałem się nad zawiłościami, ale oto podejście, które warto rozważyć. Ponieważ interesuje Cię oczekiwana liczba komórek pustych po 50 pierścieniach, musisz zdefiniować łańcuch markowa nad „liczbą pcheł w komórce”, a nie pozycję pcheł (patrz odpowiedź Glen_b, która modeluje pozycję pchła jako łańcuch markowa. Jak zauważył Andy w komentarzach do tej odpowiedzi, takie podejście może nie osiągnąć tego, czego chcesz.)

W szczególności pozwól:

$n_{ij}(t)$ $i$ $j$

Następnie łańcuch markowa zaczyna się od następującego stanu:

$n_{ij}(0) =1$ $i$ $j$

Ponieważ pchły przenoszą się do jednej z czterech sąsiednich komórek, stan komórki zmienia się w zależności od liczby pcheł znajdujących się w komórce docelowej i liczby pcheł w czterech sąsiednich komórkach oraz prawdopodobieństwa, że zostaną przeniesione do tej komórki. Za pomocą tej obserwacji możesz zapisać prawdopodobieństwo przejścia stanu dla każdej komórki w zależności od stanu tej komórki i stanu sąsiednich komórek.

Jeśli chcesz, mogę rozszerzyć odpowiedź, ale to wraz z podstawowym wprowadzeniem do łańcuchów Markowa powinno zacząć.

— Społeczność
źródło

1

n_{i j}

$n_{ij}$

@whuber Nie, nie musisz utrzymywać pozycji pcheł jako łańcucha markowa. Pomyśl o tym, co proponuję jako przypadkowy spacer do komórki. Komórka początkowo znajduje się w pozycji „1”, skąd może przejść do 0, 1, 2, 3, 4 lub 5. Prawdopodobieństwo przejścia stanu zależy od stanów sąsiednich komórek. Zatem proponowany łańcuch znajduje się w ponownie zdefiniowanej przestrzeni stanów (liczba komórek dla każdej komórki), a nie w samej pozycji pcheł. Czy to ma sens?

1

Ma to sens, ale wydaje się krokiem wstecz, ponieważ czy liczba stanów nie jest teraz znacznie większa? W jednym modelu jest 900 stanów - pozycja pojedynczej pchły - i nie więcej niż cztery przejścia z każdego z nich. Obliczenia należy wykonać tylko dla jednej pchły, ponieważ wszystkie poruszają się niezależnie. W twoim wydaje się, że stan jest opisany przez zajętość komórki wraz z zajmowaniem jej przez maksymalnie czterech sąsiadów. Byłaby to bardzo duża liczba stanów, a także bardzo duża liczba przejść między nimi. Muszę nie rozumieć, jaka jest twoja nowa przestrzeń państwowa.

— whuber

{n_{i j}}

$\{n_{ij}\}$

2

jeśli masz zamiar iść drogą numeryczną, prosta obserwacja: problem wydaje się podlegać czerwono-czarnej parzystości (pchła na czerwonym kwadracie zawsze przesuwa się na czarny kwadrat i odwrotnie). Może to pomóc zmniejszyć rozmiar problemu o połowę (rozważ tylko dwa ruchy na raz i, na przykład, patrz tylko na pchły na czerwonych kwadratach).

— shabbychef
źródło

1

To miła obserwacja. Uważam jednak, że bardziej przeszkadza, niż warto to wyraźnie wykorzystywać. Większość programowania polega na ustawieniu macierzy przejścia. Gdy to zrobisz, po prostu to wyprostuj i pracuj z tym. Dzięki zastosowaniu rzadkich macierzy usunięcie połowy zer i tak nie oszczędza czasu.

— whuber

@whuber: Podejrzewam, że sednem tych problemów jest nauczenie się technik rozwiązywania problemów, zamiast zużywania dużej liczby cykli obliczeniowych. Symetria, parzystość itp. To klasyczne techniki z książki Larsona na temat rozwiązywania problemów.

— shabbychef

1

Trafne spostrzeżenie. Ostatecznie potrzebny jest osąd. Projekt Euler wydaje się podkreślać kompromisy między wglądem matematycznym a wydajnością obliczeniową. Glen_b wspomniał o symetriach, które warto najpierw wykorzystać, ponieważ można z nich uzyskać więcej. Co więcej, stosując rzadką arytmetykę macierzy, automatycznie osiągniesz podwójny zysk (niezależnie od tego, czy znasz parzystość, czy nie!).

— whuber

1

Podejrzewam, że pewna znajomość łańcuchów Markowa w dyskretnym czasie może okazać się przydatna.

— Simon Byrne
źródło

3

To powinien był być komentarz, ale myślę, że w tym momencie możemy to zrobić.

— gung - Przywróć Monikę

Jest to automatycznie oznaczane jako niskiej jakości, prawdopodobnie dlatego, że jest tak krótkie. Czy możesz to rozwinąć?

— gung - Przywróć Monikę

Nie rozumiem, dlaczego: pytanie dotyczy tematów, które mogą być przydatne, i jest to temat, który moim zdaniem jest najbardziej odpowiedni.

— Simon Byrne

1

Oznaczono to jako niską jakość . Głosowałem, że było OK. Jeśli spojrzysz na inne odpowiedzi na ten wątek, wszystkie są znacznie dłuższe. Standardy ewoluowały w miarę upływu czasu, ale dziś byłoby to uważane za komentarz, nawet jeśli wspomina o „temacie, który może być przydatny”. Tak jak powiedziałem, myślałem, że można to zrobić tak jak jest. To, czy spróbujesz ją rozwinąć, zależy od ciebie. Po prostu dałem ci znać.

— gung - Przywróć Monikę