import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import statistics
import random
import math
import gdown
from itertools import *
random.seed(10)
#Генерируем два набора чисел размером n1 n2 из нормальных распределений с заданными параметрами
n1=1000
n2=500
set1 = np.random.normal(0, 1, n1)
set2 = np.random.normal(1, 1, n2)
#Смешиваем их в один набор, сразу сортируем его
data = np.sort(np.concatenate((set1, set2)))
sns.histplot(set1, color='orange', fill=True, label="N(0,1)")
sns.histplot(set2, color='violet', fill=True, label="N(1,1)")
plt.legend()
plt.show()
#Задаем функции для Е- и М-шага
def e_step(data, border):
mu1 = np.mean(data[:border])
sigma1 = np.std(data[:border])
mu2 = np.mean(data[border:])
sigma2 = np.std(data[border:])
return mu1, sigma1, mu2, sigma2
def m_step(data, mu1, mu2, sigma1, sigma2):
max_likelihood = -1000
sum1 = 0
sum2 = 0
for m in range(len(data)):
sum1 = 0
sum2 = 0
current_likelihood = 0
for i in range(m):
sum1 += (data[i]-mu1)**2/(2*sigma1)**2
for j in range(m+1, len(data)):
sum2 += (data[j]-mu2)**2/(2*sigma2)**2
current_likelihood = -sum1-sum2
if current_likelihood > max_likelihood:
max_likelihood = current_likelihood
border = m
return border
def density(mu, sigma, x):
dens = math.e**(-(x-mu)**2/(2*sigma**2)) / (sigma*(2*math.pi)**0.5)
return dens
#соберем шаги в ЕМ-алгоритм, в явном виде вычислим вероятности отнести число к 1 или 2 распределению
def EM_algorithm(data):
list_mu1 = []
list_mu2 = []
size_set1 = []
size_set2 = []
border = len(data)//2
for k in range(1,21):
mu1, sigma1, mu2, sigma2 = e_step(data, border)
border = m_step(data, mu1, mu2, sigma1, sigma2)
probabilities=[]
for number in data:
prob_number=[]
prob1 = (border*density(mu1, sigma1, number))/(border*density(mu1, sigma1, number) + (len(data)-border)*density(mu2, sigma2, number))
#prob 2 можно было бы вычислить 1-prob1, но тот факт, что эти два числа сложатся в единицу будет гарантией, что я не ошиблась в формуле
prob2 = ((len(data)-border)*density(mu2, sigma2, number))/(border*density(mu1, sigma1, number) + (len(data)-border)*density(mu2, sigma2, number))
prob_number.append(prob1)
prob_number.append(prob2)
probabilities.append(prob_number)
return (border, mu1, mu2, probabilities)
#добавим визуализацию каждые пять шагов, чтобы "следить за прогрессом"
def EM_algorithm_visual(data):
list_mu1 = []
list_mu2 = []
size_set1 = []
size_set2 = []
border = len(data)//2
for k in range(1,21):
mu1, sigma1, mu2, sigma2 = e_step(data, border)
border = m_step(data, mu1, mu2, sigma1, sigma2)
list_mu1.append(mu1)
list_mu2.append(mu2)
size_set1.append(border)
if k==1 or k%5==0:
print("ITERATION", k)
print( "mu1", mu1)
print( "mu2", mu2)
print( "sigma1", sigma1)
print( "sigma2", sigma2)
print( "size of set1", border)
print( "size of set2", len(data)-border)
supposed_set1 = data[:border]
supposed_set2 = data[border:]
sns.histplot(set1, color='violet', fill=True, label='true set 1')
sns.histplot(supposed_set1, color='blue', fill=True, label='predicted set 1')
plt.legend()
plt.show()
sns.histplot(set2, color='orange', fill=True, label='true set 2')
sns.histplot(supposed_set2, color='red', fill=True, label='predicted set 2')
plt.legend()
plt.show()
plt.plot(list_mu1, label='mu1')
plt.xlabel("Iterations")
plt.legend()
plt.show()
plt.plot(list_mu2, label='mu2')
plt.xlabel("Iterations")
plt.legend()
plt.show()
plt.plot(size_set1, label='size of set1')
plt.xlabel("Iterations")
plt.legend()
plt.show()
return (border, mu1, mu2, sigma1, sigma2)
EM_algorithm_visual(data)
ITERATION 1 mu1 -0.5446096914284532 mu2 1.1614921889368086 sigma1 0.6134187648088477 sigma2 0.6841361482367456 size of set1 745 size of set2 755
ITERATION 5 mu1 -0.5776942822258133 mu2 1.126412508120323 sigma1 0.6038146278666493 sigma2 0.6934039853078701 size of set1 708 size of set2 792
ITERATION 10 mu1 -0.6234521947825068 mu2 1.0812309336382575 sigma1 0.5901912605873494 sigma2 0.7051116019114307 size of set1 679 size of set2 821
ITERATION 15 mu1 -0.6326491388001118 mu2 1.0722247516967893 sigma1 0.5876076525623156 sigma2 0.7075963925550735 size of set1 671 size of set2 829
ITERATION 20 mu1 -0.6349570522545914 mu2 1.0699796363154732 sigma1 0.586961402949591 sigma2 0.7082181907851283 size of set1 670 size of set2 830
Для одномерного случая ЕМ-алгоритм может только подобрать оптимальную границу, числа до которой он назовет набором 1, а после которой - набором 2, это хорошо видно на гистограммах. В случае с нормальными распределениями с заданными параметрами ( N(0,1) и N(1,1) ), мне кажется, перекрывание хвостов слишком велико, чтобы алгоритм работал корректно.
def generate_data(n1, n2):
N = n1 + n2
set1 = np.random.normal(0, 1, n1)
set2 = np.random.normal(1, 1, n2)
data = np.sort(np.concatenate((set1, set2)))
return data
data1=generate_data(10,5)
size1, mu1, mu2, probabilities = EM_algorithm(data1)
print("ЕМ алгоритм предсказал:")
print("Размер выборки 1 - ", size1)
print("Размер выборки 2 - ", len(data1)-size1)
print("Мат ожидание mu1 - ", mu1)
print("Мат ожидание mu2 - ", mu2)
print("Вероятности каждого числа прийти из распределения 1 и распределения 2 -", probabilities)
print("--------")
print("Правильные ответы:")
print("Размер выборки 1 - 10")
print("Размер выборки 2 - 5")
print("Мат ожидание mu1 - 0")
print("Мат ожидание mu2 - 1")
ЕМ алгоритм предсказал: Размер выборки 1 - 6 Размер выборки 2 - 9 Мат ожидание mu1 - -0.6829857214402241 Мат ожидание mu2 - 1.0014445989374807 Вероятности каждого числа прийти из распределения 1 и распределения 2 - [[0.9999921428571301, 7.857142869908976e-06], [0.9999836512015062, 1.6348798493701047e-05], [0.9956000788845433, 0.004399921115456775], [0.7418664074921703, 0.2581335925078298], [0.7342243629976319, 0.26577563700236806], [0.4468239305866354, 0.5531760694133646], [0.3267454258266011, 0.6732545741733988], [0.1407233476274981, 0.859276652372502], [0.13716351333451324, 0.8628364866654867], [0.07951492538534938, 0.9204850746146506], [0.061929154718052035, 0.9380708452819481], [0.037366132966073086, 0.962633867033927], [0.013961615429088208, 0.9860383845709118], [0.0067091176907595395, 0.9932908823092403], [0.005846620332958259, 0.9941533796670418]] -------- Правильные ответы: Размер выборки 1 - 10 Размер выборки 2 - 5 Мат ожидание mu1 - 0 Мат ожидание mu2 - 1
data2=generate_data(100,50)
size1, mu1, mu2, probabilities = EM_algorithm(data2)
print("ЕМ алгоритм предсказал:")
print("Размер выборки 1 - ", size1)
print("Размер выборки 2 - ", len(data2)-size1)
print("Мат ожидание mu1 - ", mu1)
print("Мат ожидание mu2 - ", mu2)
#на html страничке эти вероятности занимают очень много места - выведу только первые десять для нагляности
print("Вероятности каждого числа прийти из распределения 1 и распределения 2 -", probabilities[:10])
print("--------")
print("Правильные ответы:")
print("Размер выборки 1 - 100")
print("Размер выборки 2 - 50")
print("Мат ожидание mu1 - 0")
print("Мат ожидание mu2 - 1")
ЕМ алгоритм предсказал: Размер выборки 1 - 59 Размер выборки 2 - 91 Мат ожидание mu1 - -0.8762737659851891 Мат ожидание mu2 - 1.100233376643748 Вероятности каждого числа прийти из распределения 1 и распределения 2 - [[0.9986954239297332, 0.0013045760702667792], [0.9982283807271941, 0.0017716192728059288], [0.9977429574948607, 0.002257042505139206], [0.9973009269498366, 0.0026990730501634487], [0.9972936448126304, 0.0027063551873695974], [0.9971186994635394, 0.0028813005364605637], [0.9961998714200562, 0.003800128579943771], [0.9958954153431863, 0.004104584656813744], [0.9950755783272539, 0.00492442167274607], [0.9939922918117275, 0.0060077081882724155]] -------- Правильные ответы: Размер выборки 1 - 100 Размер выборки 2 - 50 Мат ожидание mu1 - 0 Мат ожидание mu2 - 1
data3=generate_data(1000,500)
size1, mu1, mu2, probabilities = EM_algorithm(data3)
print("ЕМ алгоритм предсказал:")
print("Размер выборки 1 - ", size1)
print("Размер выборки 2 - ", len(data3)-size1)
print("Мат ожидание mu1 - ", mu1)
print("Мат ожидание mu2 - ", mu2)
#на html страничке эти вероятности занимают очень много места - выведу только первые десять для нагляности
print("Вероятности каждого числа прийти из распределения 1 и распределения 2 -", probabilities[:10])
print("--------")
print("Правильные ответы:")
print("Размер выборки 1 - 1000")
print("Размер выборки 2 - 500")
print("Мат ожидание mu1 - 0")
print("Мат ожидание mu2 - 1")
ЕМ алгоритм предсказал: Размер выборки 1 - 734 Размер выборки 2 - 766 Мат ожидание mu1 - -0.5710396451858859 Мат ожидание mu2 - 1.234044588284195 Вероятности каждого числа прийти из распределения 1 и распределения 2 - [[0.9999945232521766, 5.47674782342883e-06], [0.9999851435687469, 1.4856431253075374e-05], [0.999982411689097, 1.758831090291378e-05], [0.9999792692835755, 2.0730716424456953e-05], [0.9999576531124774, 4.234688752252021e-05], [0.999940462615171, 5.95373848289209e-05], [0.9999345308555373, 6.546914446272247e-05], [0.9999226197704194, 7.738022958061423e-05], [0.9999153004441598, 8.469955584018715e-05], [0.999892850823654, 0.00010714917634596944]] -------- Правильные ответы: Размер выборки 1 - 1000 Размер выборки 2 - 500 Мат ожидание mu1 - 0 Мат ожидание mu2 - 1
Три упражнения продемонстрировали, что проблема была не в малом размере выборки (с увеличением количества элементов качество не улучшилось), а все же в перекрывании данных - средние лежат на расстоянии в стандартное отклонение, из-за чего нельзя четко разграничить их одной прямой
(Пример, показывающий, что алгоритм все-таки работает)
На наборах данных, у которых средние отличаются больше, ЕМ-алгоритм работает корректнее.
n1 = 1000
n2 = 500
N = n1 + n2
random.seed(10)
set1 = np.random.normal(0, 1, n1)
set2 = np.random.normal(5, 1, n2)
data_test = np.sort(np.concatenate((set1, set2)))
sns.histplot(set1, color='orange', fill=True, label="N(0,1)")
sns.histplot(set2, color='violet', fill=True, label="N(5,1)")
plt.legend()
plt.show()
EM_algorithm_visual(data_test)
ITERATION 1 mu1 -0.3819378123653121 mu2 3.773742028443601 sigma1 0.7301003842571164 sigma2 1.973739541634068 size of set1 779 size of set2 721
ITERATION 5 mu1 -0.20717026444257006 mu2 4.338384765370443 sigma1 0.8053136889456561 sigma2 1.6399068708665174 size of set1 912 size of set2 588
ITERATION 10 mu1 -0.03139007787213524 mu2 4.894348094878663 sigma1 0.9236995761806746 sigma2 1.1366982074305512 size of set1 980 size of set2 520
ITERATION 15 mu1 -0.004637704382946664 mu2 4.958027117860511 sigma1 0.949204162362557 sigma2 1.0697594037542155 size of set1 986 size of set2 514
ITERATION 20 mu1 -0.004637704382946664 mu2 4.958027117860511 sigma1 0.949204162362557 sigma2 1.0697594037542155 size of set1 986 size of set2 514
(986, -0.004637704382946664, 4.958027117860511, 0.949204162362557, 1.0697594037542155)
#Задаем исходные параметры модели
alphabet = 'ab'
real_init_prob = {'a' : 0.5,
'b' : 0.5}
real_prob = {"aa" : 0.6,
"ab" : 0.4,
"ba" : 0.3,
"bb" : 0.7}
# Функция для генерации последовательности
def generate_sequence(length):
sequence=""
letter = 'a' if random.random() < real_init_prob['a'] else 'b'
for _ in range(length):
sequence+=letter
prev_letter=letter
if prev_letter == 'a':
letter = 'a' if random.random() < real_prob['aa'] else 'b'
else:
letter = 'a' if random.random() < real_prob['ba'] else 'b'
return sequence
#сразу зададим все функции, которые нужны будут, чтобы оценить параметры марковской модели любого порядка, а не только первого
#функция подсчета количества к-меров до определенной позиции
def count_kmers(k, alphabet, sequences, end):
#создаем словарь со всеми возможными к-мерами заданной длиным из букв нашего алфавита
kmers=[]
for i in product(alphabet, repeat=k):
kmer = "".join(i)
kmers.append(kmer)
kmers_dict = dict.fromkeys(kmers, 0)
#считаем количество каждого к-мера
for seq in sequences:
for kmer in kmers:
counter=0
for i in range(end-k+1):
if seq[i:i+k]==kmer:
counter+=1
kmers_dict[kmer] += counter
return kmers_dict
#функция подсчета частоты к-меров до определенной позиции
def kmers_freq(k, sequences, end):
kmers_dict = count_kmers(k, alphabet, sequences, end)
prefix_dict = count_kmers(k-1, alphabet, sequences, end-1)
kmer_freq_dict = dict.fromkeys(kmers_dict.keys(), 0)
for kmer in kmers_dict.keys():
prefix = kmer[:-1]
kmer_freq_dict[kmer] = "{:.3f}".format(kmers_dict[kmer] / prefix_dict[prefix])
return kmer_freq_dict
#функция оценки параметров бернуллиевской модели
def estimate_bernoulli_model(sequences, alphabet, end):
freq_dict = dict.fromkeys(list(alphabet), 0)
for seq in sequences:
for letter in freq_dict.keys():
freq_dict[letter] += seq.count(letter, 0, end)
for letter in freq_dict.keys():
freq_dict[letter] = freq_dict[letter]/(len(sequences)*end)
return freq_dict
#функция оценки параметров марковской модели любого порядка
def estimate_markov_model(order, sequences):
#считаем распределения вероятностей для первых order-1 символов
probs = [estimate_bernoulli_model(sequences, alphabet, 1)]
for i in range(2, order+1):
probs.append( kmers_freq(i, sequences, i+1) )
kmers_freq_dict1 = estimate_bernoulli_model(sequences, alphabet, 1)
#распределения вероятностей для остальных символов
probs.append( kmers_freq(order+1, sequences, len(sequences[0])) )
return probs
# Генерируем 1000 последовательностей длиной 3000 нуклеотидов
sequences = [generate_sequence(3000) for _ in range(1000)]
estimate_markov_model(1, sequences)
[{'a': 0.497, 'b': 0.503},
{'aa': '0.599', 'ab': '0.401', 'ba': '0.300', 'bb': '0.700'}]
Параметры хорошо соотносятся с теми, что я задала выше:
{'a' : 0.5, 'b' : 0.5}
{"aa" : 0.6, "ab" : 0.4, "ba" : 0.3, "bb" : 0.7}
Задаем необходимые функции
def ln_likelihood_Markov(probabilities, order, sequences):
L=0
for seq in sequences:
#print(seq)
L+=math.log(len(seq), math.e)
for i in range(order):
#print(seq[0:i+1])
L+=math.log(float(probabilities[i][seq[0:i+1]]), math.e)
#print("-------")
for i in range(order, len(seq)):
L+=math.log(float(probabilities[order][seq[i-order:i]+seq[i]]), math.e)
return L
def ln_likelihood_bernoulli(probabilities, sequences):
L = 0
for seq in sequences:
L+=math.log(len(seq), math.e)
for letter in seq:
L+=math.log(probabilities[letter], math.e)
return L
def indep_params(alphabet, order):
return (len(alphabet)**(order+1))-1
def BIC(n, k, lnL):
#k - число параметров модели
#n - размер обучающей выборки
#lnL - ln(likelihood)
return k * math.log(n, math.e) - 2*lnL
def AIC(k, lnL):
#n - размер обучающей выборки
#lnL - ln(likelihood)
return 2*k - 2*lnL
#ln_likelihood_Markov(kmers_freq_dict, 3, test_seq)
Оцениваем параметры для бернуллиевской модели и марковский моделей 1-3 порядков. Из построения мы знаем, что все последовательности строились одинаково, поэтому перемешивать их для разделения на тестовую и тренировочную выборки необязательно.
train_set = sequences[:800]
test_set = sequences[800:]
alphabet = 'ab'
params = estimate_bernoulli_model(train_set, alphabet, len(train_set[0]))
k = indep_params(alphabet, 0)
lnL = ln_likelihood_bernoulli(params, test_set)
bic = BIC(len(train_set), k, lnL)
aic = AIC(k, lnL)
print("AIC", aic)
print("BIC", bic)
AIC 816418.6689926077 BIC 816423.3536043355
params = estimate_markov_model(order=1, sequences=train_set)
k = indep_params(alphabet, 1)
lnL = ln_likelihood_Markov(params, order=1, sequences=test_set)
bic = BIC(len(train_set), k, lnL)
aic = AIC(k, lnL)
print("AIC", aic)
print("BIC", bic)
AIC 761385.1288494183 BIC 761399.1826846013
params = estimate_markov_model(order=2, sequences=train_set)
k = indep_params(alphabet, order=2)
lnL = ln_likelihood_Markov(params, order=2, sequences=test_set)
bic = BIC(len(train_set), k, lnL)
aic = AIC(k, lnL)
print("AIC", aic)
print("BIC", bic)
AIC 761392.584674897 BIC 761425.3769569907
params = estimate_markov_model(order=3, sequences=train_set)
k = indep_params(alphabet, order=3)
lnL = ln_likelihood_Markov(params, order=3, sequences=test_set)
bic = BIC(len(train_set), k, lnL)
aic = AIC(k, lnL)
print("AIC", aic)
print("BIC", bic)
AIC 761410.1336455714 BIC 761480.4028214865
AIC
Бернуллиевская модель: 816418.6689926077
Марковская модель 1-го порядка: 761385.1288494183 (Победитель!)
Марковская модель 2-го порядка: 761392.584674897
Марковская модель 3-го порядка: 761410.1336455714
BIC
Бернуллиевская модель: 816423.3536043355
Марковская модель 1-го порядка: 761399.1826846013 (Победитель!)
Марковская модель 2-го порядка: 761425.3769569907
Марковская модель 3-го порядка: 761480.4028214865
Действительно, баланса между правдоподобием и количеством параметров достигла марковская модель 1-го порядка, которая и была использована для генерации последовательностей :)
my_input = open('genome_bacillus.fna')
lines = my_input.readlines()
genome_seq=""
for line in lines:
if line[0]!=">":
line=line.replace("\n", "")
genome_seq+=line
print(len(genome_seq))
genome_pieces=[]
for i in range(len(genome_seq)//1000):
genome_pieces.append(genome_seq[i*1000:(i+1)*1000])
print(len(genome_pieces))
4215606 4215
В случае с геномом уже стоит перемешать полученные последовательности
random.shuffle(genome_pieces)
train_pieces = genome_pieces[:3000] #будем обучать на 4000 последовательностях, а тестировать на остальных
test_pieces = genome_pieces[3000:]
alphabet = "ATGC"
params = estimate_bernoulli_model(train_pieces, alphabet, len(train_pieces[0]))
k = indep_params(alphabet, 0)
lnL = ln_likelihood_bernoulli(params, test_pieces)
bic = BIC(len(train_pieces), k, lnL)
aic = AIC(k, lnL)
print("AIC", aic)
print("BIC", bic)
{'A': 0.281867, 'T': 0.28293033333333334, 'G': 0.21687633333333334, 'C': 0.21832633333333334}
AIC 3331344.947942344
BIC 3331362.967045047
for i in range(1, 6):
print("Марковская модель", i, '-го порядка')
params = estimate_markov_model(order=i, sequences=train_pieces)
k = indep_params(alphabet, order=i)
lnL = ln_likelihood_Markov(params, order=i, sequences=test_pieces)
bic = BIC(len(train_pieces), k, lnL)
aic = AIC(k, lnL)
print("AIC", aic)
print("BIC", bic)
Марковская модель 1 -го порядка AIC 3290410.9096070216 BIC 3290501.0051205363 Марковская модель 2 -го порядка AIC 3267764.051151953 BIC 3268142.4523087153 Марковская модель 3 -го порядка AIC 3257974.5465277163 BIC 3259506.1702574673 Марковская модель 4 -го порядка AIC 3255188.2092244066 BIC 3261332.723246113 Марковская модель 5 -го порядка
--------------------------------------------------------------------------- ValueError Traceback (most recent call last) <ipython-input-158-6c266ffd3b81> in <cell line: 1>() 3 params = estimate_markov_model(order=i, sequences=train_pieces) 4 k = indep_params(alphabet, order=i) ----> 5 lnL = ln_likelihood_Markov(params, order=i, sequences=test_pieces) 6 7 bic = BIC(len(train_pieces), k, lnL) <ipython-input-133-bbbb0db358eb> in ln_likelihood_Markov(probabilities, order, sequences) 6 for i in range(order): 7 #print(seq[0:i+1]) ----> 8 L+=math.log(float(probabilities[i][seq[0:i+1]]), math.e) 9 #print("-------") 10 for i in range(order, len(seq)): ValueError: math domain error
С марковской моделью пятого порядка возникла ошибка, так как в тренировочной выборке не встретился один из 5-меров, а его частота была приравнена 0. Но в тестовой выборке этот 5-мер возник, а логарифм от него взять не вышло. Можно было бы добавить "псевдочастоты", равные не 0, а, например, 10^(-10). Это бы несильно повлияло на остальную статистику, но позволило бы избежать ошибки. Но так как по BIC уже видно, что модель четвертого порядка уступает третьей, я решила остановится.
AIC
Бернуллиевская модель: 3331344.947942344
Марковская модель 1-го порядка: 3290410.9096070216
Марковская модель 2-го порядка: 3267764.051151953
Марковская модель 3-го порядка: 3257974.5465277163
Марковская модель 4-го порядка: 3255188.2092244066 (Победитель!)
BIC
Бернуллиевская модель: 3331362.967045047
Марковская модель 1-го порядка: 3290501.0051205363
Марковская модель 2-го порядка: 3268142.4523087153
Марковская модель 3-го порядка: 3259506.1702574673 (Победитель!)
Марковская модель 4-го порядка: 3261332.723246113
По AIC оценке геном лучше всего описывается марковской моделью 4-го порядка, по BIC - марковской моделью 3-го порядка.