Probabilidades - quão alto você pode ir?

10

Eu já fiz uma pergunta sobre como calcular uma probabilidade com rapidez e precisão. No entanto, evidentemente, foi muito fácil, pois foi fornecida uma solução de formulário fechado! Aqui está uma versão mais difícil.

Esta tarefa é sobre escrever código para calcular uma probabilidade exata e rapidamente . A saída deve ser uma probabilidade precisa escrita como uma fração em sua forma mais reduzida. Ou seja, nunca deve produzir, 4/8mas sim 1/2.

Para um número inteiro positivo n, considere uma sequência uniformemente aleatória de 1s e -1s de comprimento ne chame-a de A. Agora concatene Auma cópia de si mesma. Ou seja, A[1] = A[n+1]se estiver indexando a partir de 1, A[2] = A[n+2]e assim por diante. Aagora tem comprimento 2n. Agora considere também uma segunda sequência aleatória de comprimento ncujos primeiros nvalores sejam -1, 0 ou 1 com probabilidade 1 / 4,1 / 2, 1/4 cada e chame-a de B.

Agora considere o produto interno de Bcom A[1+j,...,n+j]para diferente j =0,1,2,....

Por exemplo, considere n=3. Valores possíveis para Ae Bpoderiam ser A = [-1,1,1,-1,...]e B=[0,1,-1]. Nesse caso, os dois primeiros produtos internos são 0e 2.

Tarefa

Para cada um j, começando com j=1, seu código deve gerar a probabilidade de que todos os primeiros j+1produtos internos sejam zero para todos n=j,...,50.

Copiando a tabela produzida por Martin Büttner para j=1os seguintes resultados de amostra.

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

Ponto

Sua pontuação é a maior que jseu código é concluído em 1 minuto no meu computador. Para esclarecer um pouco, cada jum recebe um minuto. Observe que o código de programação dinâmica na pergunta vinculada anterior fará isso facilmente j=1.

Desempate

Se duas entradas obtiverem a mesma jpontuação, a entrada vencedora será a que atingir a maior pontuação nem um minuto na minha máquina j. Se as duas melhores entradas também forem iguais nesse critério, o vencedor será a resposta enviada primeiro.

Línguas e bibliotecas

Você pode usar qualquer idioma e bibliotecas disponíveis gratuitamente que desejar. Devo ser capaz de executar seu código, portanto, inclua uma explicação completa de como executar / compilar seu código no Linux, se possível.

Minha máquina Os horários serão executados na minha máquina. Esta é uma instalação padrão do ubuntu em um processador AMD FX-8350 de oito núcleos. Isso também significa que eu preciso poder executar seu código.

Entradas vencedoras

  • j=2em Python por Mitch Schwartz.
  • j=2em Python por feersum. Atualmente, a entrada mais rápida.
Comunidade
fonte
Se a pergunta não estiver clara, informe-me para que eu possa corrigi-la rapidamente.
2
Você é de longe a minha pergunta favorita. Então, novamente, eu tenho uma coisa para calcular valores exatamente e rapidamente .
06
@primo Obrigado! Isso significa que podemos esperar uma resposta no RPython? :)
Você poderia colocar a diferença entre essa pergunta e a outra?
kirbyfan64sos
@ kirbyfan64sos A outra é essencialmente a mesma pergunta apenas para `j = 1`.

Respostas:

3

Python 2, j = 2

Tentei encontrar uma espécie de 'formulário fechado' para j = 2. Talvez eu pudesse criar uma imagem do MathJax, embora fosse realmente feio com todo o manuseio do índice. Eu escrevi esse código não otimizado apenas para testar a fórmula. Demora cerca de 1 segundo para concluir. Os resultados correspondem ao código de Mitch Schwartz.

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)

Considere uma sequência em que o i-ésimo membro é ese A [i] == A [i + 1] ou nse A [i]! = A [i + 1]. ino programa é o número de ns. Se ifor par, a sequência deve começar e terminar com e. Se ifor ímpar, a sequência deve começar e terminar com n. As sequências são ainda classificadas pelo número de execuções de es ou s consecutivos n. Existem j+ 1 de um e jdo outro.

Quando a ideia passeio aleatório é estendido a 3 dimensões, não é um problema lamentável que existem 4 possíveis direcções para andarem (um para cada uma das ee, en, ne, or nn) o que significa que não são linearmente dependente. Portanto, o kíndice soma o número de etapas executadas em uma das direções (1, 1, 1). Depois, haverá um número exato de etapas que devem ser seguidas nas outras três direções para cancelá-lo.

P (n, p) fornece o número de partições ordenadas de n objetos em p partes. W (l, d) fornece o número de maneiras para uma caminhada aleatória de letapas percorrer uma distância exata d. Como antes, há 1 chance de se mover para a esquerda, 1 chance de se mover para a direita e 2 para permanecer parado.

feersum
fonte
Obrigado! Uma imagem da fórmula seria realmente ótima.
Obrigado pela explicação. Você faz parecer simples! Acabei de ver seu comentário para o qual você poderia fazer uma solução j=3. Isso seria incrível!
3

Python, j = 2

A abordagem de programação dinâmica, j = 1na minha resposta à pergunta anterior, não precisa de muitas modificações para funcionar mais alto j, mas fica lenta rapidamente. Tabela para referência:

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

E o código:

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()

Aqui nós estamos mantendo o controle dos dois primeiros elementos A, os dois últimos elementos B(onde b2é o último elemento), e os produtos internos de (A[:n], B), (A[1:n], B[:-1])e (A[2:n], B[:-2]).

Mitch Schwartz
fonte
.... atingiu N_MAX com 21,20 segundos restantes