Probabilités - jusqu'où pouvez-vous aller?

10

J'ai déjà posé une question sur la façon de calculer une probabilité rapidement et avec précision. Cependant, c'était évidemment trop facile car une solution sous forme fermée a été donnée! Voici une version plus difficile.

Cette tâche consiste à écrire du code pour calculer une probabilité exactement et rapidement . La sortie doit être une probabilité précise écrite sous forme de fraction dans sa forme la plus réduite. C'est-à-dire qu'il ne devrait jamais sortir 4/8mais plutôt 1/2.

Pour un entier positif n, considérez une chaîne uniformément aléatoire de 1 et -1 de longueur net appelez-la A. Maintenant, concaténez à Aune copie d'elle-même. C'est A[1] = A[n+1]si l'indexation à partir de 1, A[2] = A[n+2]et ainsi de suite. Aa maintenant une longueur 2n. Considérons maintenant également une deuxième chaîne aléatoire de longueur ndont les premières nvaleurs sont -1, 0 ou 1 avec probabilité 1 / 4,1 / 2, 1/4 chacune et appelons-la B.

Considérez maintenant le produit intérieur de Bavec A[1+j,...,n+j]pour différent j =0,1,2,....

Par exemple, réfléchissez n=3. Les valeurs possibles pour Aet Bpourraient être A = [-1,1,1,-1,...]et B=[0,1,-1]. Dans ce cas, les deux premiers produits intérieurs sont 0et 2.

Tâche

Pour chacun j, en commençant par j=1, votre code devrait afficher la probabilité que tous les premiers j+1produits internes soient nuls pour tous n=j,...,50.

En copiant le tableau produit par Martin Büttner pour j=1nous avons les exemples de résultats suivants.

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

But

Votre score est le plus élevé que jvotre code complète en 1 minute sur mon ordinateur. Pour clarifier un peu, chacun jobtient une minute. Notez que le code de programmation dynamique de la question liée précédente le fera facilement j=1.

Tie break

Si deux entrées obtiennent le même jscore, l'entrée gagnante sera celle qui atteint le plus haut nen une minute sur ma machine pour cela j. Si les deux meilleures entrées sont également égales sur ce critère, le gagnant sera la réponse soumise en premier.

Langues et bibliothèques

Vous pouvez utiliser n'importe quelle langue et bibliothèque librement disponibles. Je dois être en mesure d'exécuter votre code, veuillez donc inclure une explication complète sur la façon d'exécuter / compiler votre code sous Linux si possible.

Ma machine Les horaires seront exécutés sur ma machine. Il s'agit d'une installation ubuntu standard sur un processeur AMD FX-8350 à huit cœurs. Cela signifie également que je dois pouvoir exécuter votre code.

Entrées gagnantes

  • j=2en Python par Mitch Schwartz.
  • j=2en Python par feersum. Actuellement l'entrée la plus rapide.
Communauté
la source
Si la question n'est pas claire de quelque façon que ce soit, faites-le-moi savoir afin que je puisse la résoudre rapidement.
2
Vous êtes de loin mon questionneur préféré. Là encore, j'ai une chose pour calculer les valeurs exactement et rapidement .
primo
@primo Merci! Est-ce à dire que nous pouvons nous attendre à une réponse dans RPython? :)
Pourriez-vous mettre la différence entre cette question et l'autre?
kirbyfan64sos
@ kirbyfan64sos L'autre est essentiellement la même question pour «j = 1» seulement.

Réponses:

3

Python 2, j = 2

J'ai essayé de trouver une sorte de «forme fermée» pour j = 2. Je pourrais peut-être en faire une image MathJax, même si ce serait vraiment moche avec tous les tripotages d'index. J'ai écrit ce code non optimisé uniquement pour tester la formule. Il faut environ 1 seconde pour terminer. Les résultats correspondent au code 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)

Considérons une séquence où le ième membre est esi A [i] == A [i + 1] ou nsi A [i]! = A [i + 1]. idans le programme est le nombre de ns. Si iest pair, la séquence doit commencer et se terminer par e. Si iest impair, la séquence doit commencer et se terminer par n. Les séquences sont en outre classées selon le nombre d'exécutions de es ou ns consécutifs . Il y a j+1 de l'un et jde l'autre.

Lorsque l'idée de marche aléatoire est étendue à 3 dimensions, il y a un problème regrettable qu'il ya 4 directions possibles pour marcher dans (un pour chacun ee, en, neou nn) ils ne sont pas linéairement dépendants qui moyens. Ainsi, l' kindice totalise le nombre de pas effectués dans l'une des directions (1, 1, 1). Ensuite, il y aura un nombre exact d'étapes qui doivent être prises dans les 3 autres directions pour l'annuler.

P (n, p) donne le nombre de partitions ordonnées de n objets en p parties. W (l, d) donne le nombre de façons dont une marche aléatoire de lpas parcourt une distance exacte d. Comme auparavant, il y a 1 chance de se déplacer à gauche, 1 chance de se déplacer à droite et 2 de rester en place.

feersum
la source
Je vous remercie! Une image de la formule serait vraiment géniale.
Merci pour l'explication. Vous le faites paraître simple! Je viens de voir votre commentaire pour lequel vous pourriez trouver une solution j=3. Ce serait génial!
3

Python, j = 2

L'approche de programmation dynamique pour j = 1dans ma réponse à la question précédente n'a pas besoin de beaucoup de modifications pour fonctionner plus haut j, mais devient lente rapidement. Tableau de référence:

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

Et le code:

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

Nous gardons ici la trace des deux premiers éléments de A, les deux derniers éléments de B(où b2est le dernier élément), et les produits internes de (A[:n], B), (A[1:n], B[:-1])et (A[2:n], B[:-2]).

Mitch Schwartz
la source
.... atteint N_MAX avec 21,20 secondes restantes