Вероятности - как высоко вы можете пойти?

10

Ранее я задавал вопрос о том, как быстро и точно вычислить вероятность. Тем не менее, очевидно, что это было слишком легко, так как было дано решение в закрытой форме! Вот более сложная версия.

Эта задача о написании кода для точного и быстрого вычисления вероятности . Вывод должен быть точной вероятностью, записанной в виде дроби в наиболее сокращенной форме. То есть это никогда не должно выводиться, 4/8а скорее 1/2.

Для некоторого положительного целого числа nрассмотрим равномерно случайную строку длиной 1 с и 1 с nи назовите ее A. Теперь объедините ее Aс копией. То есть A[1] = A[n+1]если индексировать с 1 A[2] = A[n+2]и тд. Aтеперь имеет длину 2n. Теперь также рассмотрим вторую случайную строку длины n, первые nзначения которой равны -1, 0 или 1 с вероятностью 1 / 4,1 / 2, 1/4 каждая и назовем ее B.

Теперь рассмотрим внутреннее произведение Bс A[1+j,...,n+j]по разному j =0,1,2,....

Например, рассмотрим n=3. Возможные значения Aи Bмогут быть A = [-1,1,1,-1,...]и B=[0,1,-1]. В этом случае первые два внутренних продукта являются 0и 2.

задача

Для каждого j, начиная с j=1, ваш код должен выводить вероятность того, что все первые j+1внутренние продукты равны нулю для каждого n=j,...,50.

Копируя таблицу, произведенную Мартином Бюттнером, j=1мы получаем следующий пример.

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

Гол

Ваша оценка - самая большая сумма, которую jваш код выполняет за 1 минуту на моем компьютере. Чтобы уточнить немного, каждый jполучает одну минуту. Обратите внимание, что динамический программный код в предыдущем связанном вопросе сделает это легко для j=1.

Стяжка

Если две записи получают одинаковый jсчет, то выигравшая запись будет той, которая достигнет максимума nза одну минуту на моей машине j. Если по этому критерию две лучшие записи равны, то победителем будет ответ, представленный первым.

Языки и библиотеки

Вы можете использовать любой свободно доступный язык и библиотеки, которые вам нравятся. Я должен быть в состоянии запустить ваш код, поэтому, пожалуйста, включите полное объяснение того, как запустить / скомпилировать ваш код в Linux, если это вообще возможно.

Моя машина Время будет запущено на моей машине. Это стандартная установка Ubuntu на восьмиъядерный процессор AMD FX-8350. Это также означает, что мне нужно иметь возможность запустить ваш код.

Победные записи

  • j=2в питоне Митч Шварц.
  • j=2в Python с помощью feersum. На данный момент самая быстрая запись.
Сообщество
источник
Если вопрос неясен каким-либо образом, пожалуйста, дайте мне знать, чтобы я мог быстро это исправить.
2
Вы, безусловно, мой любимый вопрос. С другой стороны, у меня есть возможность точно и быстро вычислить значения .
Примо
@primo Спасибо! Значит ли это, что мы можем ожидать ответа в RPython? :)
Не могли бы вы поставить разницу между этим вопросом и другим?
kirbyfan64sos
@ kirbyfan64sos Другой по сути тот же вопрос только для `j = 1`.

Ответы:

3

Python 2, j = 2

Я пытался найти что-то вроде «закрытой формы» для j = 2. Возможно, я мог бы сделать это на MathJax-образе, хотя это было бы очень уродливо, если бы теребили индекс. Я написал этот неоптимизированный код только для проверки формулы. Это займет около 1 секунды. Результаты совпадают с кодом Митча Шварца.

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)

Рассмотрим последовательность, в которой i-й член равен, eесли A [i] == A [i + 1] или nесли A [i]! = A [i + 1]. iв программе есть число nс. Если iчетное, последовательность должна начинаться и заканчиваться на e. Если iнечетно, последовательность должна начинаться и заканчиваться n. Последовательности дополнительно классифицируются по количеству последовательных последовательностей es или ns. Есть j+1 от одного и jдругого.

Когда идея случайной ходьбы увеличиваются до 3 -х измерениях, есть неудачная проблема , что существует 4 возможных направления , чтобы идти в ( по одному для каждого из ee, en, ne, или nn) означает , что они не являются линейно зависимыми. Таким образом, kиндекс суммирует количество шагов, выполненных в одном из направлений (1, 1, 1). Затем будет точное количество шагов, которые необходимо предпринять в трех других направлениях, чтобы отменить его.

P (n, p) дает количество упорядоченных разбиений n объектов на p частей. W (l, d) дает количество способов случайного обхода lшагов, чтобы пройти расстояние точно d. Как и прежде, есть 1 шанс двигаться влево, 1 шанс двигаться вправо и 2 остаться на месте.

feersum
источник
Спасибо! Изображение формулы было бы действительно здорово.
Спасибо за объяснение. Вы делаете это выглядит просто! Я только что видел ваш комментарий, для которого вы могли бы найти решение j=3. Это было бы удивительно!
3

Python, J = 2

Подход динамического программирования j = 1в моем ответе на предыдущий вопрос не требует много модификаций для работы на более высоком уровне j, но быстро становится медленным. Таблица для справки:

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

И код:

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

Здесь мы отслеживанием первых двух элементов A, последние два элемента B(где b2это последний элемент), и внутренние продукты (A[:n], B), (A[1:n], B[:-1])и (A[2:n], B[:-2]).

Митч Шварц
источник
.... достиг N_MAX с оставшимися