Prawdopodobieństwa - jak wysoko możesz się udać?

10

Wcześniej zadałem pytanie, jak szybko i dokładnie obliczyć prawdopodobieństwo. Jednak najwyraźniej było to zbyt łatwe, ponieważ podano rozwiązanie w formie zamkniętej! Oto trudniejsza wersja.

To zadanie dotyczy pisania kodu w celu dokładnego i szybkiego obliczenia prawdopodobieństwa . Wynik powinien być precyzyjnym prawdopodobieństwem zapisanym jako ułamek w najbardziej zredukowanej formie. Oznacza to, że nigdy nie powinien generować, 4/8ale raczej 1/2.

Dla pewnej dodatniej liczby całkowitej nrozważ jednolicie losowy ciąg 1s i -1s długości ni nazwij ją A. Teraz połącz się Az kopią samego siebie. To znaczy, A[1] = A[n+1]jeśli indeksowanie od 1 A[2] = A[n+2]i tak dalej. Ateraz ma długość 2n. Teraz rozważ także drugi losowy ciąg długości, nktórego pierwsze nwartości to -1, 0 lub 1 z prawdopodobieństwem 1 / 4,1 / 2, 1/4 każdy i nazwij go B.

Teraz rozważ wewnętrzny produkt Bz A[1+j,...,n+j]dla różnych j =0,1,2,....

Rozważmy na przykład n=3. Możliwe wartości Ai Bmogą być A = [-1,1,1,-1,...]i B=[0,1,-1]. W tym przypadku pierwsze dwa produkty wewnętrzne to 0i 2.

Zadanie

Dla każdego j, zaczynając od j=1, twój kod powinien wyświetlać prawdopodobieństwo, że wszystkie pierwsze j+1wewnętrzne produkty są zerowe dla każdego n=j,...,50.

Kopiując tabelę wyprodukowaną przez Martina Büttnera j=1, mamy następujące przykładowe wyniki.

n   P(n)
1   1/2
2   3/8
3   7/32
4   89/512
5   269/2048
6   903/8192
7   3035/32768
8   169801/2097152

Wynik

Twój wynik jest najwyższy, jaki jTwój kod wypełnia w ciągu 1 minuty na moim komputerze. Aby wyjaśnić trochę, każda jdostaje jedną minutę. Zauważ, że kod programowania dynamicznego w poprzednim połączonym pytaniu zrobi to łatwo j=1.

Łamacz krawatów

Jeśli dwa zgłoszenia otrzymają ten sam jwynik, zwycięskim wejściem będzie ten, który osiągnie najwyższą wartość nw ciągu jednej minuty na mojej maszynie j. Jeśli dwa najlepsze zgłoszenia są równe również w tym kryterium, zwycięzcą będzie odpowiedź przesłana jako pierwsza.

Języki i biblioteki

Możesz korzystać z dowolnego darmowego języka i bibliotek, które ci się podobają. Muszę być w stanie uruchomić Twój kod, więc proszę podać pełne wyjaśnienie, jak uruchomić / skompilować kod w systemie Linux, jeśli to w ogóle możliwe.

Moja maszyna Czasy zostaną uruchomione na moim komputerze. Jest to standardowa instalacja ubuntu na ośmiordzeniowym procesorze AMD FX-8350. Oznacza to również, że muszę być w stanie uruchomić Twój kod.

Zwycięskie zgłoszenia

  • j=2w Pythonie autor: Mitch Schwartz.
  • j=2w Python przez feersum. Obecnie najszybszy wpis.
Społeczność
źródło
Jeśli pytanie jest niejasne w jakikolwiek sposób, daj mi znać, żebym mógł je szybko rozwiązać.
2
Jesteś zdecydowanie moim ulubionym pytającym. Z drugiej strony, mam coś do precyzyjnego i szybkiego obliczania wartości .
primo
@primo Dziękujemy! Czy to oznacza, że ​​możemy spodziewać się odpowiedzi w RPython? :)
Czy potrafisz odróżnić to pytanie od drugiego?
kirbyfan64sos
@ kirbyfan64sos Drugi to zasadniczo to samo pytanie tylko dla `j = 1`.

Odpowiedzi:

3

Python 2, j = 2

Próbowałem znaleźć rodzaj „zamkniętej formy” dla j = 2. Być może mógłbym zrobić z niego obraz MathJax, chociaż byłoby to naprawdę brzydkie z tym całym indeksowaniem. Napisałem ten niezoptymalizowany kod tylko w celu przetestowania formuły. Wykonanie zajmuje około 1 sekundy. Wyniki są zgodne z kodem Mitcha Schwartza.

ch = lambda n, k: n>=k>=0 and fac[n]/fac[k]/fac[n-k]
W = lambda l, d: ch(2*l, l+d)
P = lambda n, p: n==0 or ch(n-1, p-1)
ir = lambda a,b: xrange(a,b+1)

N = 50
fac = [1]
for i in ir(1,4*N): fac += [i * fac[-1]]

for n in ir(2, N):
    s = 0
    for i in ir(0,n+1):
     for j in ir(0,min(i,n+1-i)):
      for k in ir(0,n+i%2-i-j):
       t=(2-(k==0))*W(n+i%2-i-j,k)*W(i-(j+i%2),k)*W(j,k)**2*P(i,j+i%2)*P(n+1-i,j+1-i%2)
       s+=t
    denp = 3 * n - 1
    while denp and not s&1: denp -= 1; s>>=1
    print n, '%d/%d'%(s,1<<denp)

Rozważ sekwencję, w której ósmym elementem jest, ejeśli A [i] == A [i + 1] lub njeśli A [i]! = A [i + 1]. iw programie jest liczba ns. Jeśli ijest parzyste, sekwencja musi zaczynać się i kończyć na e. Jeśli ijest nieparzyste, sekwencja musi zaczynać się i kończyć na n. Sekwencje są dalej klasyfikowane według liczby przebiegów kolejnych es lub ns. Jest j+ 1 jednego i jdrugiego.

Gdy błądzenia losowego pomysł zostanie przedłużony do 3 wymiarach, jest niefortunne problem, że są 4 możliwe kierunki chodzić (po jednym dla każdego ee, en, ne, lub nn), co oznacza, że nie są liniowo zależne. Tak więc kindeks sumuje się przez liczbę kroków wykonanych w jednym z kierunków (1, 1, 1). Następnie będzie dokładna liczba kroków, które należy wykonać w pozostałych 3 kierunkach, aby go anulować.

P (n, p) podaje liczbę uporządkowanych partycji n obiektów na części p. W (l, d) podaje liczbę sposobów losowego przejścia lkroków na pokonanie odległości dokładnie d. Tak jak poprzednio, istnieje 1 szansa na ruch w lewo, 1 szansa na ruch w prawo i 2 na utrzymanie pozycji.

feersum
źródło
Dziękuję Ci! Obraz tej formuły byłby naprawdę świetny.
Dziękuję za wyjaśnienie. Sprawiasz, że wygląda to prosto! Właśnie zobaczyłem twój komentarz, dla którego możesz znaleźć rozwiązanie j=3. To by było wspaniałe!
3

Python, j = 2

Podejście do programowania dynamicznego j = 1w mojej odpowiedzi na poprzednie pytanie nie wymaga wielu modyfikacji, aby pracować na wyższym poziomie j, ale szybko zwalnia. Tabela w celach informacyjnych:

n   p(n)

2   3/8
3   11/64
4   71/512
5   323/4096
6   501/8192
7   2927/65536
8   76519/2097152
9   490655/16777216
10  207313/8388608

I kod:

from time import*
from fractions import*
from collections import*

def main():
    N_MAX=50

    T=time()

    n=2
    Y=defaultdict(lambda:0)
    numer=0

    for a1 in [1]:
        for b1 in (1,0):
            for a2 in (1,-1):
                for b2 in (1,0,0,-1):
                    if not a1*b1+a2*b2 and not a2*b1+a1*b2:
                        numer+=1
                    Y[(a1,a2,b1,b2,a1*b1+a2*b2,a2*b1,0)]+=1

    thresh=N_MAX-1

    while time() <= T+60:
        print('%d %s'%(n,Fraction(numer,8**n/4)))

        if thresh<2:
            print('reached N_MAX with %.2f seconds remaining'%(T+60-time()))
            return

        n+=1
        X=Y
        Y=defaultdict(lambda:0)
        numer=0

        for a1,a2,b1,b2,s,t,u in X:
            if not ( abs(s)<thresh and abs(t)<thresh+1 and abs(u)<thresh+2 ):
                continue

            c=X[(a1,a2,b1,b2,s,t,u)]

            # 1,1

            if not s+1 and not t+b2+a1 and not u+b1+a1*b2+a2: numer+=c
            Y[(a1,a2,b2,1,s+1,t+b2,u+b1)]+=c

            # -1,1

            if not s-1 and not t-b2+a1 and not u-b1+a1*b2+a2: numer+=c
            Y[(a1,a2,b2,1,s-1,t-b2,u-b1)]+=c

            # 1,-1

            if not s-1 and not t+b2-a1 and not u+b1+a1*b2-a2: numer+=c
            Y[(a1,a2,b2,-1,s-1,t+b2,u+b1)]+=c

            # -1,-1

            if not s+1 and not t-b2-a1 and not u-b1+a1*b2-a2: numer+=c
            Y[(a1,a2,b2,-1,s+1,t-b2,u-b1)]+=c

            # 1,0

            c+=c

            if not s and not t+b2 and not u+b1+a1*b2: numer+=c
            Y[(a1,a2,b2,0,s,t+b2,u+b1)]+=c

            # -1,0

            if not s and not t-b2 and not u-b1+a1*b2: numer+=c
            Y[(a1,a2,b2,0,s,t-b2,u-b1)]+=c

        thresh-=1

main()

Tutaj jesteśmy śledzenie dwóch pierwszych elementów A, dwa ostatnie elementy B(gdzie b2jest ostatni element) i wewnętrzne produkty (A[:n], B), (A[1:n], B[:-1])oraz (A[2:n], B[:-2]).

Mitch Schwartz
źródło
.... osiągnął N_MAX z pozostałą 21,20 sekundą