Вот если сместить учет пар на одно деление вправо, получим:
[‘A’, ‘GG’, ‘TG’, ‘AT’, ‘GG’, ‘A’]
Это то, о чем я говорил во втором вопросе. По идее уникальных делений на пары только два, так что может все таки можно как-то найти решение через этот способ?
Я тогда тоже поделюсь своими пока незаконченными рассуждениями. Можно разделить все нежелательные исходы на три группы, которые не пересекаются:
Для первого и третьего случая посчитать нетрудно. В первом случае это просто 3^{100}, а в третьем — 3^{98}. Сейчас я пытаюсь понять, как можно найти число комбинаций для второго случая. Пробую посчитать это число для коротких последовательностей и потом вывести для общего случая, при длине в n букв.
А это прям красиво получилось. Я попытался продолжить твой ход, но пришел к еще одной проблеме. Мы не учитываем порядок пар. Он тоже может меняться. И казалось бы, что можно просто умножить на число перестановок этих 49 пар и 1 нуклеотида, но тогда вопрос — сколько таких перестановок?
У нас какие-то пары будут встречаться по несколько раз, такие варианты надо будет отсекать, но их количество зависит от того, какая пара повторяется и сколько раз она повторяется. Может получиться, что все 49 пар являются ‘AA’, а может получиться, что ‘AA’ повторяется 40 раз и ‘TT’ повторяется 9 раз. В каждом случае считать по-разному.
Вчера у меня стоял похожий вопрос и тогда я подумал, а что если через какое-то распределение выразить то, сколько будет какая-то пара повторяться? Тогда к каждому такому случаю можно было бы написать свое произведение, но это звучит как-то нереально, не могу даже представить, что там выйдет.
ED: я сейчас понял, что нам порядок и не важен, тогда все нормально
если мы считаем все возможности пар, то это уже включает повторяющиеся пары и их перестановки. Так что домножать не нужно)
Как насчет такого конечного ответа:
49 пар без ‘GG’. Потом 49-я пара либо имеет G на конце (тогда последним может быть любой из остальных трех), либо не имеет. И возвести в квадрат, потому что и при первом делении, и при втором, это все должно выполняться. А сами деления симметричны, поэтому можно взять одинаковые вероятности для них.
Во второй скобке поставил 15, потому что слева степень 49-я, значит уже все пары без ‘GG’, то есть выбирать осталось только среди 15 вариантов. 16 должно стоять если слева будет 48-я степень, мне кажется.
Если они одинаковые, то надо помножить на два, а не делить.
Не совсем понял откуда ты взял 12/15. Если продолжать мою мысль, то финальный ответ:
Где k=1 или k=2 в зависимости от того перемножаем мы на два за счет другого деления или нет. На 4 множим т.к. первый нуклеотид (в последовательности из 100) может быть любым из четырех.
Ответ @tardigrade 0.999\,999\,245
Если посчитать мое с k=2: 0.730\, 180\, 873
Давай теперь подумаем насколько это физически имеет смысл?
Откуда деление? Я возвел в квадрат, потому что это можно представить так: мы делим все на пары либо первым способом, либо вторым. Теперь они, вроде как, независимы. Нам надо, чтобы выполнялось “не” для обоих, поэтому надо их перемножить. Поскольку случаи симметричны, получилось возведение в квадрат.
Вероятность того, что последняя пара не заканчивается на G. Тогда последний нуклеотид может быть любым, то есть умножение на 1.
В 4 раза увеличится количество событий, как всех, так и нужных, поэтому я думаю, надо без умножения. И второй аргумент — две формулировки задачи эквивалентны.
Вот это самое интересное. Со вчерашнего дня у меня крутится мысль о том, чтобы написать код для симуляции. Можно написать функцию, которая рандомно выбирает букву из листа [‘A’, ‘T’, ‘G’, ‘C’], и ставит их подряд пока не получится 100 букв. Потом смотрит есть ли в этой строке ‘NGG’ и увеличивает счетчик нужных событий. Только есть опасения, что симуляция окажется слишком долгой.
Ой, да, квадрат.
А, ну да, логично
Я в том же пути рассуждал, но получил 51/64 суммарно.
О, давай Монте Карло запустим! Современные компьютеры достаточно быстрые для этого, даже на ноуте быстро будет!
В случае, когда последняя пара заканчивается на ‘G’, ты умножил на \frac{1}{4}, но надо было на \frac{3}{4}. И еще у меня в знаменателях стоит 15, а не 16, потому что в предыдущей скобке мы уже учли вероятность того, что все пары, включая последнюю, не являются ‘GG’. Тогда уже всего вариантов для последней пары не 16, а 15.
Как правильно отметил @Anton
нам достаточно решить задачу “какова вероятность того, что в последовательности из 99 нуклеотидов не будет подряд идущих двух G”.
Давайте решим более общую задачу – какова вероятность того, что в последовательности из n нуклеотидов не будет подряд идущих двух G. Для этого найдем количество последовательностей из n нуклеотидов, таких, что в них нет GG, потом поделим на 4^n – на количество всех возможных последовательностей из n нуклеотидов. Так и получим вероятность.
Чтобы это найти, введем две вспомогательные функции: пусть f(n) – количество последовательностей длины n в которых нет GG и в которой последний нуклеотид G; g(n) – количество последовательностей длины n в которых нет GG и последний нуклеотид не G. Тогда ответом к нашей задаче будет число
Попробуем найти рекуррентную зависимость этих двух функций.
С одной стороны, мы можем получить последовательность длины n которая не заканчивается на G, если к последовательности длины n-1 заканчивающееся на G припишем нуклеотид A,C,T (3 способа) или к последовательности длины n-1 не заканчивающееся на G припишем нуклеотид A,C,T (3 способа). Следовательно
С другой стороны, мы можем получить последовательность длины n которая заканчивается на G, если к последовательности длины n-1 не заканчивающееся на G припишем нуклеотид G (1 способ). Тогда имеем
Соединим два равенства и распишем базовые случаи для n=1,n=2. Получаем
Выпишем явную форму данного рекуррентного уравнения (можно сделать через характеристические уравнения, а можно через вольфрам). Получается вот такое несложное выражение 
Далее выражаем f(n) = g(n-1) .
Выходит, что ответ равен
prev = 3
cur = 12
for i in range(3, 100):
temp = cur
cur = 3 * temp + 3 * prev
prev = temp
ans = (cur + prev) / 4**99
print(1 - ans)
посчитал. получил что вероятность нахождения NGG равна \approx 0.9948085, что в целом похоже на правду
Параллельно, я сделал функцию, которая считает по твоему точному результату:
from math import sqrt
def recurrent(n):
denominator = 7 * (3 + sqrt(21)) * 4**n
term1 = 2*sqrt(21)*(3/2-sqrt(21)/2)**n
term2 = 2*sqrt(21)*(3/2-sqrt(21)/2)**(n-1)
term3 = (21+5*sqrt(21))*(3/2+sqrt(21)/2)**n
term4 = (21+5*sqrt(21))*(3/2+sqrt(21)/2)**(n-1)
numerator = term1 + term2 + term3 + term4
return numerator/denominator
def solve(n):
return 1 - recurrent(n)
# arman = sum(res:=[solve(100) for i in range(10)])/len(res)
# print(arman)
Получаю 0.9950793852734019 (официальное решение / то, что указала Балжан дает 0.999999245)
Дальше, я запустил Монте-Карло на задачу.
Cначала написал функцию проверки наличия фрагмента в последовательности
def isPresent(sequence, fragment):
'''
Inputs:
- sequence (str): a DNA sequence (or just a general string)
- fragment (str): two chars for which we perform the check
Output:
bool - whether el is present in sequence or not
'''
# TO DO: generalize function to general fragment???
last = None
for nucleotide in sequence:
if not last:
last = nucleotide
continue
if last + nucleotide == fragment:
return True
else:
last = nucleotide
return False
# SMAL DEBUG
# t1 = 'CATG'
# t2 = 'CGGT'
# t3 = 'CAGATGC'
# t4 = 'CAGATGG'
# for t in [t1, t2, t3, t4]:
# print(isPresent(t, 'GG'))
Дальше строим рандомную цепочку
import random
NUCLEOTIDES = ['A', 'G', 'C', 'T']
def build_sequence(elts, length):
'''
builds a string composed randomly from elts of given length
'''
seq = ''
for _ in range(length):
seq += random.choice(elts)
return seq
И наконец сама Монте Карло. Строим 1000 последовательностей десять раз:
def estimate_probability(length, elts, fragment, iterations):
'''
In N iterations, estimates probability that fragment is present in a random sequence of length made from elts
'''
fragmentPresent = 0
fragmentNotPresent = 0
for i in range(iterations):
print(f"building seq {i} / {iterations}")
randomSequence = build_sequence(elts, length)
if isPresent(randomSequence, fragment):
fragmentPresent += 1
else:
fragmentNotPresent += 1
return (fragmentPresent)/(fragmentPresent + fragmentNotPresent)
def average_probability(length, elts, fragment, iterations, runs):
'''
Perform N runs of estimate_probability
'''
return sum(probs:=[estimate_probability(length, elts, fragment, iterations) for _ in range(runs)])/len(probs)
length = 100
monteCarlo = average_probability(length, NUCLEOTIDES, 'GG', 1000, 10)
theory = theory.solve(length)
print(f"\n{monteCarlo} (MonteCarlo)\n{theory} (Theory)\n{theory-monteCarlo} (diff)\n")
В итоге, получаем
0.9953999999999998 (MonteCarlo)
0.9950793852734019 (Theory)
-0.0003206147265979453 (diff)
Я запустил Монте Карло еще пару раз. Среднее значение вероятности получалось проводя 10 экспериментов, в которых создавалось N последовательностей ДНК, длиной в 100 нуклеотидов:
| N | Время работы | Вероятность | Истинное значение (теор.) | Отклонение |
|---|---|---|---|---|
| 100 | 0.074 с | 0.991 | 0.9950793852734019 | 0.004079385273401903 |
| 1 000 | 0.633 с | 0.9949 | 0.9950793852734019 | 0.0001793852734018886 |
| 10 000 | 5.924 с | 0.99506 | 0.9950793852734019 | 1.938527340195062e-05 |
| 100 000 | 59.021 с | 0.9950130000000001 | 0.9950793852734019 | 6.638527340174782e-05 |
| 1 000 000 | 5849.2 с | 0.9950823200000001 | 0.9950793852734019 | -2.9347265981805037e-06 |
В общем, я думаю мы нашли правильное решение))
P.S. 5849.2 секунды это 97.5 минут
Думаю будет полезно пояснить откуда возникли эти страшные корни. На самом деле все просто и в процессе выведения явных решений линейных рекурсий нет ничего сложного.
Пусть нам дано рекуррентное соотношение
и известны значения a_1,a_2,p,q. Пусть x_1, x_2 – корни уравнения x^2 - px - q=0. Данное уравнение называется характеристическим уравнением рекуррентной последовательности. Тогда должны существовать такие A и B, что
Для нахождения A, B просто подставим n=1 и n=2 и решим систему из двух уравнений, откуда получим явную формулу
Конкретно в нашем случае мы имеем
Отсель получаем характеристическое уравнение – x^2 - 3x - 3 = 0 \implies x_{1,2} = \frac{3}{2} \pm \frac{\sqrt{21}}{2}. Следовательно
Остается найти A,B. Подставляем n=1, n=2
Решаем простенькую систему и получаем, что
Далее подставляем в (*) и вуаля!
Подробнее про рекурсию можно почитать тутъ
Это все точно биология???!!!
ага, просто в более жестокой форме
сколько лет пользовался аском и не натыкался на это прекрасное обсуждение… спасибо что подняли