Implémentation Python du modèle Markov caché continu

Depuis que j'ai commencé à publier qiita cette fois, j'écrirai sur l'implémentation du modèle Markov caché continu par python, que j'ai eu du mal à implémenter il y a quelques mois. En mémoire, j'écrirai à partir des bases dans l'ordre que j'ai étudié. De plus, je pense que ce résultat est probablement atteint, mais il y a des éléments suspects, donc si vous remarquez quelque chose, je vous serais reconnaissant si vous pouviez me le faire savoir.

Quel est le modèle de Markov?

Le modèle de Markov est un modèle créé sur la base du processus de Markov, et le processus de Markov est un processus stochastique avec des propriétés de Markov. Markov est sorti beaucoup et je ne pouvais pas le comprendre facilement, alors je vais expliquer ma propre compréhension. Tout d'abord, considérons un diagramme comme la théorie des graphes ci-dessous. スクリーンショット 2020-09-15 19.21.15.png

Sur cette figure, on suppose qu'après un certain laps de temps, l'état change de manière probabiliste dans la direction vers laquelle pointe la flèche, et cette probabilité est le nombre attaché à la flèche. S'il est maintenant dans l'état A, il passera à l'état B avec une probabilité de 0,3 et à l'état C avec une probabilité de 0,7. Si l'état A provient de X ou Y, la probabilité de transition de l'état A à l'état B ou à l'état C ne change pas. En d'autres termes, la propriété selon laquelle la probabilité de transition vers l'état suivant (B ou C) dépend uniquement de l'état actuel (A) est appelée propriété de Markov. Le processus stochastique est écrit comme "une variable stochastique qui change avec le temps.-Wikipedia-" Dans ce cas, la probabilité de prendre un état change avec le temps, on peut donc dire que c'est un processus stochastique. ~~ Excusez-moi, je laisserai la définition stricte à d'autres sites d'information ~~ Un modèle créé de cette manière peut être appelé un modèle de Markov.

Modèle de Markov caché

Le modèle de Markov caché est un modèle comme le modèle de Markov mentionné précédemment, mais c'est un modèle utilisé lorsqu'il a une certaine valeur (sortie) en fonction de l'état au lieu de ne pas connaître l'état. Faisons un exemple concret basé sur le traitement du langage naturel. Tout d'abord, considérons l'exemple de phrase (sortie) "Je réserve un film" (je réserve un film). A cette époque, «je» et «film» sont la nomenclature, «livre» est un verbe et «a» est un acronyme. Ces mots de partie sont considérés comme cachés car ils ne sont pas connus de cette phrase. Supposons que vous souhaitiez que la machine détermine la partie de cette phrase. Cependant, bien que «je» et «a» soient presque uniquement déterminés, «livre» et «film» peuvent également être utilisés comme verbes tels que prendre un livre ou un film, qui est une nomenclature. Par conséquent, examinons s'il est possible de classer les pièces en utilisant le modèle de Markov caché.

En anglais, l'ordre des parties est décidé dans une certaine mesure, comme 5 modèles de phrases. Par conséquent, la facilité de changement des mots de pièce est probablement celle indiquée dans le tableau ci-dessous. (Le prix est fixé par votre propre sentiment)

Début de phrase nom verbe article Fin de phrase
Début de phrase 0 0.6 0.2 0.2 0
nom 0 0.3 0.5 0.1 0.1
verbe 0 0.6 0 0.3 0.1
article 0 1 0 0 0
Fin de phrase 0 0 0 0 0

Considérez la gauche comme l'état actuel et le haut comme l'état suivant. Pour donner un exemple, la probabilité qu'un nez vienne à côté d'un nez est de 0,3, la probabilité qu'un verbe vienne est de 0,5, la probabilité qu'un acronyme vienne de 0,1 et la probabilité que la fin d'une phrase vienne est de 0,1. ~~ Utiliser 0 pour la probabilité n'est pas bon pour un problème de fréquence nulle ~~ Ensuite, on suppose que la sortie est prise avec la probabilité suivante dans chaque état. スクリーンショット 2020-09-16 0.30.15.png

En d'autres termes, lorsqu'il s'agit d'une nomenclature, "I" vaut 0,4, "a" vaut 0,0, "livre" vaut 0,3 et "film" vaut 0,2. La méthode de transition d'état à ce moment est la suivante. (L'état où la probabilité de transition est 0 est supprimé à l'avance)

スクリーンショット 2020-09-16 1.50.16.png

Ensuite, le nombre de mots est n, la chaîne de mots (observable) est $ W = w_1, w_2, ..., w_n $, et la chaîne de partie (non observable: état) est $ T = t_1, t_2,. .., t_n $, la probabilité que ce modèle donnera lors du choix d'un itinéraire est P(w)= \prod_{i=1}^{n}P(t_i|t_{i-1})P(w_i|t_i)

Peut être représenté par. En d'autres termes, la probabilité de passer de l'état précédent à l'état actuelP(t_i|t_{i-1})La probabilité de transition qui est et la probabilité que le mot soit sorti dans l'état d'une certaine partieP(w_i|t_i)Il peut être obtenu en additionnant le produit des probabilités d'apparition. Dans ce modèle, l'état initial (état initial) était le début de la phrase. Il est évident que le premier état est cette fois le début de la phrase, mais lors de l'adaptation à d'autres données, de quel état le premier état est susceptible de commencer (probabilité initiale) est également fortement lié à la probabilité. Par conséquent, la probabilité résultante de ce modèle est $ \ pi_i $ pour la probabilité initiale, $ a_ {i, j} $ pour la probabilité de transition et $ b_j (w_t) $ pour la fonction de sortie (cependant, i, j représentent un certain état. t indique l'ordre de sortie.) Et la probabilité lorsque l'ordre des états est $ I = i_1, i_2, ..., i_n $ est P(w)=\pi_{i_0}\prod_{t=1}^{n}a_{i_{t-1},i_t}b_{i_t}(w_t)

(Cependant, dans l'exemple de la phrase utilisée cette fois, $ w_0 $ est vide, $ w_1 $ est'I ', ..., $ w_n $ est vide.) Ensuite, la partie de cette phrase peut être prédite en trouvant la méthode de transition avec la probabilité la plus élevée. En passant, dans le cas de cet exemple, on peut dire qu'il a été biaisé de manière facile à comprendre, mais on pense que la probabilité de transition dans l'ordre du début de la phrase → nez → verbe → acronyme → nez → fin de phrase est la plus élevée, et les mots des parties sont correctement classés. ..

De côté La probabilité initiale continuera à être $ \ pi $ , La probabilité de transition est $ a_ {i, j} $, et la probabilité d'apparition est $ b_j (w) $.

Modèle de Markov caché continu

Puisque le modèle de Markov caché mentionné précédemment pensait à quatre mots, la sortie peut être considérée comme une sortie à quatre valeurs. La sortie était donc une valeur discrète. Par conséquent, on peut dire qu'il s'agit d'un modèle de Markov caché discret. Par conséquent, le modèle de Markov caché lors de l'utilisation de valeurs continues telles que la température est appelé le modèle de Markov caché continu. Dans le modèle de Markov discret mentionné précédemment, la probabilité d'apparition de chaque état est calculée individuellement, mais comme il n'est pas réaliste d'attribuer toutes les valeurs à des valeurs continues, nous définirons la probabilité d'apparition en utilisant une distribution gaussienne mixte. .. Avant d'entrer dans l'implémentation du modèle de Markov caché continu, écrivons d'abord sur la distribution gaussienne mixte.

Distribution gaussienne / distribution gaussienne mixte

La distribution gaussienne mixte est une représentation de la distribution de probabilité en ajoutant un certain nombre de distributions gaussiennes. En premier lieu, la distribution gaussienne est également appelée distribution normale. N(x| \mu,\sigma)=\frac{1}{\sqrt{2 \pi\sigma^2}}\exp(-\frac{(x-\mu)^2}{2\sigma^2})

Il est représenté par. 300px-Standard_deviation_diagram.svg.png

-Image empruntée à wikipedia- Comme vous pouvez le voir sur ce graphique, plus ce x est proche de la moyenne $ \ mu $ (0 sur cette figure), plus la valeur est élevée. Et si la dispersion est grande, elle sera plate. En d'autres termes, si cette valeur est considérée comme une probabilité, plus la valeur de x est proche de la moyenne $ \ mu $, plus la valeur sera élevée, et si la variance est grande, il sera possible de prendre cette valeur même si elle est loin de la moyenne. Je peux dire.

De plus, la distribution gaussienne mixte revient à ajouter plusieurs distributions normales ensemble, et si K distributions gaussiennes sont combinées, la formule de la distribution gaussienne mixte est P(x|\pi,\mu,\sigma)=\sum_{k=1}^{K}\pi_kN(x|\mu_k,\sigma_k)

Cependant, il existe une condition de $ \ sum_ {k = 1} ^ K \ pi_k = 1 $.

Il est représenté par. Si vous essayez d'illustrer cela スクリーンショット 2020-09-17 23.27.20.png

Ça ressemble à ça.

Clustering avec distribution gaussienne mixte

À l'avenir, lors de la mise en œuvre du modèle de Markov caché continu, il est nécessaire d'identifier la distribution gaussienne à partir des points de données. Par conséquent, j'aimerais voir le flux en implémentant le clustering doux avec une distribution gaussienne mixte.

Commençons par tracer les points en fonction des trois distributions gaussiennes. スクリーンショット 2020-09-18 14.57.21.png

Voici l'intrigue. Ces points sont dans un état où l'on sait à partir de quelle distribution gaussienne ils ont été générés. (Données complètes) Cependant, lors du clustering, je ne sais pas à partir de quelle distribution il a été généré. (Données incomplètes) Par conséquent, à partir de maintenant, nous ajusterons les paramètres de la distribution gaussienne mixte constituée des trois distributions gaussiennes afin qu'elles puissent être au mieux divisées en trois selon la distribution gaussienne.

Comment pouvons-nous bien le diviser en trois? Cela signifie que la probabilité simultanée de tous les points selon la distribution gaussienne mixte souhaitée doit être la plus élevée. En d'autres termes ln(p(X|\pi,\mu,\sigma))= ln(\prod_{i=1}^{n}\sum_{k=1}^K \pi_kN(x_i|\mu_k,\sigma_k))=\sum_{i=1}^nln(\sum_{k=1}^K\pi_kN(x_i|\mu_k,\sigma_k))

$ X = [x_1, x_2, ..., x_n] $ ・ ・ ・ Chaque point de données $ K $ ・ ・ ・ Nombre de modèles (nombre de clusters)

Ajustez les paramètres pour que ce soit le plus élevé. Ceci est également appelé la fonction de vraisemblance logarithmique. Je voudrais changer les paramètres afin de maximiser cette fonction, mais comme la somme de la somme apparaît et qu'il est difficile de la résoudre analytiquement, je vais la résoudre en utilisant l'algorithme EM.

En tant que méthode d'algorithme EM Étape E: calcul de la valeur attendue Étape M: Ajustez les paramètres pour maximiser la valeur attendue C'est une méthode pour trouver la solution optimale en répétant.

Le flux cette fois est

  1. Déterminez la valeur initiale de $ \ pi, \ mu, \ sigma $
  2. Calculez le taux de charge (étape E)
  3. Ajustez les valeurs des paramètres $ \ pi, \ mu, \ sigma $.
  4. S'il n'y a pas de différence de probabilité, le processus se termine. Si cela change toujours, revenez à 2.

Je vais le faire dans le flux.

valeur initiale##

La valeur initiale initiale est décidée de manière appropriée. (Bien sûr, si vous connaissez cette valeur d'une manière ou d'une autre, il semble qu'elle convergera plus rapidement si vous l'ajustez.)

E étape

Vient ensuite l'étape E. Le terme taux de charge vient ici, qui est la probabilité qu'un point x se produise à partir du modèle k. Compte tenu du graphique ci-dessus, je ne savais pas de quel modèle de distribution gaussienne venait. Considérez la variable latente z indiquant si cela est né d'un modèle pour cela. Si $ z_ {i, k} = 1 $, alors le point $ x_i $ est un point né du modèle k.

Quand cela devient une expression p(z_{i,k}=1|x_i)=\frac{p(z_{i,x}=1,x_i)}{p(x_i)}

Aussi, p(x_i|x_{i,k}=1)=\frac{p(z_{i,k}=1,x_i)}{p(z_{i,k}=1)}

Parce que p(z_{i,k}=1|x_i)=\frac{p(x_i|z_{i,k}=1)p(z_{i,z}=1)}{p(x_i)}

Est. (Théorème de Bayes) Aussi, puisque $ p (x_i) $ est une valeur marginalisée par rapport à z, elle est exprimée par la formule de distribution gaussienne mixte. p(z_{i,k}=1|x_i)=\frac{p(x_i|z_{i,k}=1)p(z_{i,k}=1)}{\sum_{k=1}^{K}\pi_kN(x|\mu_k,\sigma_k)}=\gamma(z_{i,k})

Ce sera. ici,p(x_i|z_{i,k}=1)Est du modèle kx_iDepuis la probabilité qui est née, la distribution normale du modèle kN(x_i|\mu_k,\sigma_k)est. Et parce que $ p (z_ {i, k} = 1) = \ pi_k $ \gamma(z_{i,k})=\frac{\pi_kN(x|\mu_k,\sigma_k)}{\sum_{k=1}^{K}\pi_kN(x|\mu_k,\sigma_k)} Ce sera.

Je n'ai pas compris pourquoi $ p (z_ {i, k} = 1) = \ pi_k $, mais c'est une valeur qui n'a rien à voir avec x, et à quel point elle appartient au modèle k. C'est un sentiment que cela peut certainement être le cas s'il s'agit d'un montant à comparer, que ce soit facile ou non.)

La raison pour laquelle on appelle cela le taux de charge est qu'il montre à quel point le modèle k occupe la proportion de la distribution gaussienne mixte. Si vous regardez le chiffre d'il y a deux ans, vous pouvez l'imaginer d'une manière ou d'une autre.

M étape

Vient ensuite la mise à jour des paramètres. Comme pour chaque politique, pour $ \ mu et \ sigma $, la fonction de vraisemblance la plus élevée est suffisante, nous allons donc procéder à la politique de différenciation et de recherche du point où la valeur devient 0. De plus, puisque $ \ pi $ a la condition $ \ sum_ {k = 1} ^ K \ pi_k = 1 $, il sera résolu par la méthode du multiplicateur indéterminé de Lagrange.

La dérivation de cette différenciation est Explication approfondie de l'algorithme EM, donc d'autres personnes l'expliquent d'une manière très facile à comprendre, donc je ne résumerai que les résultats. ~~ Je l'ai personnellement dérivé, mais c'était très difficile ~~ résultat \mu_k*=\frac{1}{N_k}\sum_{i=1}^N\gamma(z_{i,k})x_i \sigma_k*=\frac{1}{N_k}\sum_{i=1}^N\gamma(z_{i,k})(x_i-\mu_k)(x_i-\mu_k)^T \pi_k*=\frac{N_k}{N}=\frac{\sum_{n=1}^N\gamma(z_{i,k})}{N}

Le résultat est.

Implémentation de la distribution gaussienne mixte

Depuis que la dérivation a été résumée, j'ai en fait mis en œuvre le clustering de la distribution gaussienne mixte.

ex_gmm.py


import numpy as np
import matplotlib.pyplot as plt



#difine
data_num = 100 
model_num = 3
dim = 2
color = ['#ff0000','#00ff00','#0000ff']
epsilon = 1.0e-3
check = False # if not check == False

#gaussian calculator
def cal_gauss(x,mu,sigma):
	return np.exp(-np.dot((x-mu),np.dot(np.linalg.inv(sigma),(x-mu).T))/2)/np.sqrt(2*np.pi*abs(np.linalg.det(sigma)))



#init create param
mu =[]
sigma = []
for i in range(model_num):
	mu.append([np.random.rand()*10 ,np.random.rand()*10])
	temp = np.random.rand()-np.random.rand()
	sigma.append([[np.random.rand()+1,temp],[temp,np.random.rand()+1]])
	


#create data
values = []

for i in range(model_num):
	values.append(np.random.multivariate_normal(mu[i],sigma[i],data_num))
for i in range(model_num):
	plt.scatter(values[i][:,0],values[i][:,1],c=color[i],marker='+')

plt.show()



#create training data
train_values =values[0]
for i in range(model_num-1):
	train_values = np.concatenate([train_values,values[i+1]])
plt.scatter(train_values[:,0],train_values[:,1],marker='+')




#init model param
model_mu =[]
model_sigma =[]
 
for i in range(model_num):
	model_mu.append([np.random.rand()*10 ,np.random.rand()*10])
	temp = np.random.rand()-np.random.rand()
	model_sigma.append([[np.random.rand()+1,temp],[temp,np.random.rand()+1]])
	
model_mu = np.array(model_mu)
model_sigma = np.array(model_sigma)
for i in range(model_num):
	plt.scatter(model_mu[i,0],model_mu[i,1],c=color[i],marker='o') #plot pre train mu
temp = []
for i in range(model_num):
	temp.append(np.random.rand())
model_pi = np.zeros_like(temp)
for i in range(model_num):
	model_pi[i] = temp[i]/np.sum(temp) 

plt.show()



#training
old_L = 1
new_L = 0
gamma = np.zeros((data_num*model_num,model_num))
while(epsilon<abs(new_L-old_L)):

	#gamma calculation
	for i in range(data_num*model_num):
		result_gauss = []
		for j in range(model_num):
			result_gauss.append(cal_gauss(train_values[i,:],model_mu[j],model_sigma[j]))	
		result_gauss = np.array(result_gauss)
		old_L = np.sum(model_pi*result_gauss) #old_L calculation	
		for j in range(model_num):
			gamma[i,j] = (model_pi[j]*result_gauss[j])/np.sum(model_pi*result_gauss,axis=0)

	N_k =np.sum(gamma,axis=0)
	
	#calculation sigma
	for i in range(model_num):
		temp = np.zeros((dim,dim))
		for j in range(data_num*model_num):	
			temp+=gamma[j,i]*np.dot((train_values[j,:]-model_mu[i]).reshape(-1,1),(train_values[j,:]-model_mu[i]).reshape(1,-1))
	
		model_sigma[i] = temp/N_k[i]


	#calculation mu
	for i in range(model_num):
		model_mu[i] = np.sum(gamma[:,i].reshape(-1,1)*train_values,axis=0)/N_k[i]


	#calculation pi
	for i in range(model_num):
		model_pi[i]=N_k[i]/(data_num*model_num)

	#check plot
	if(check == True ):
		gamma_crasta = np.argmax(gamma,axis=1)	
		for i in range(data_num*model_num):
			color_temp = color[gamma_crasta[i]]
			plt.scatter(train_values[i,0],train_values[i,1],c=color_temp,marker='+')
		
		for i in range(model_num):
			plt.scatter(model_mu[i,0],model_mu[i,1],c=color[i],marker='o')

		plt.show()
	##new_Lcalculation
	for i in range(data_num*model_num):
		result_gauss = []
		for j in range(model_num):
			result_gauss.append(cal_gauss(train_values[i,:],model_mu[j],sigma[j]))
		result_gauss = np.array(result_gauss)
		new_L = np.sum(model_pi*result_gauss)
	
gamma_crasta = np.argmax(gamma,axis=1)	
for i in range(data_num*model_num):
	color_temp = color[gamma_crasta[i]]
	plt.scatter(train_values[i,0],train_values[i,1],c=color_temp,marker='+')

for i in range(model_num):
	plt.scatter(model_mu[i,0],model_mu[i,1],c=color[i],marker='o')

plt.show()


Cela devrait fonctionner comme ça. Dans ce programme, tous les paramètres sont sortis de manière aléatoire, il peut donc ne pas être possible de bien regrouper. ~~ En fait, lorsque de telles données arrivent, elles ne peuvent pas être regroupées ... ~~

Implémentation du modèle de Markov caché continu

Cela a été un détour, mais j'aimerais voir désormais l'implémentation du modèle Markov continu. Lorsque l'information d'observation y = {$ y_1, y_2, ..., y_t $}, qui est un modèle de Markov continu, est obtenue, la probabilité $ p (y | M) $ que ce modèle rappelle cette valeur est déterminée par l'algorithme en treillis. peut faire. L'algorithme de treillis $ p (y | M) = \ sum_ {i} \ alpha (i, T) $ (Cependant, ce i est l'état final)

Est requis à. Et ce $ \ alpha (i, t) $ est

\begin{eqnarray}
\alpha_{i,t}=\left\{ \begin{array}{ll}
\ \pi_i & (t=0) \\
\ \sum_j \alpha(j,t-1)a_{j,i}・ B_{i}(y_t) & (|l-k|=1) \\
\end{array} \right.
\end{eqnarray}

est. Et $ a_ {j, i} $ et $ b_ {j, i} (y_t) $ qui apparaissent ici sont il y a peu de temps, mais ils montrent la probabilité de transition et la probabilité d'apparition qui apparaissaient également dans le modèle de Markov caché. Je suis. (Cette histoire n'a rien à voir avec l'apprentissage ...) Ce $ \ alpha (i, t) $ trace la probabilité depuis l'avant (fonction avant), mais définit une fonction qui trace cela depuis la fin (fonction arrière).

\begin{eqnarray}
\beta_{i,T}=\left\{ \begin{array}{ll}
\ 1 & (État final) \\
\ 0 & (Pas dans l'état final) \\
\end{array} \right.
\end{eqnarray}
\beta(i,t)=\sum_{j}a_{i,j}・ B_{i}(y_t)・\beta(j,t+1) \qquad (t=T,T-1,...,1)

Les paramètres du modèle du modèle de Markov caché peuvent être appris à l'aide de l'algorithme de Baum Welch, qui est un type d'algorithme EM. L'algorithme de Baumwelch est utilisé pour calculer en E étapes. $ \ Gamma (i, j, t) = \ frac {\ alpha (i, t-1) ・ a_ {i, j} ・ b_ {j} ・ \ beta (j, t)} {p (y | M )} $

Cette valeur. Cette valeur est considérée comme la probabilité de passer de l'état i à j au temps t et de sortir $ y_t $ divisée par la probabilité que y se produise dans ce modèle.

Et la formule d'estimation (mise à jour) des paramètres est

\pi*=\frac{\sum_j\gamma(i,j,1)}{\sum_i\sum_j\gamma(i,j,1)}

a_{i,j}*=\frac{\sum_t\gamma(i,j,t)}{\sum_t\sum_j\gamma(i,j,t)}

Suivant, Parce que b a une distribution gaussienne

\mu_{i,j}*=\frac{\sum_{t=1}^T\gamma(i,j,t)y_t}{\sum_{t=1}^T\gamma(i,j,t)}

\sigma_{i,j}=\frac{\sum_{t=1}^T \gamma(i,j,t)(y_t-\mu_j)(y_t-\mu_j)^t}{\sum_{t=1}^T\gamma(i,j,t)}

Est requis à.

Ensuite, c'est un code. Un fichier d'implémentation de la distribution gaussienne car elle a été ajustée séquentiellement. Il est divisé en trois parties: un fichier modèle Markov caché et un fichier exécutable. (Au début, je prévoyais de le faire avec GMM-HMM, il y a donc beaucoup de pièces qui sont inutiles à implémenter. Veuillez patienter un moment car je le mettrai à jour à nouveau dès que GMM-HMM fonctionnera.)

class_b.py


import numpy as np
import matplotlib.pyplot as plt
import math
import sys

#construct
#c	=	mixing coefficient	
#mu	=	mean
#sigma	=	dispersion
#m	=	model num
#size	=	datasize
#x	=	data

##method
#gaussian(x)	:	calculation gaussian:x=data
#updata(x)	:	update gmm :x=data
#output(x)	:	gmm output:x=data (size=1)
class function_b:
	upsilon = 1.0e-50
	def __init__(self,c,mu,sigma,m,size,x):

		self.c = np.array(c)
		self.c = self.c.reshape((-1,1))

		self.mu = np.array(mu)
		self.mu = self.mu.reshape((-1,1))

		self.sigma = np.array(sigma)
		self.sigma = self.sigma.reshape((-1,1))

		self.size = size
		self.model_num = m
		self.x = np.array(x).reshape(1,-1)
		self.ganma = np.zeros((self.model_num,self.size))
		self.sum_gaus = 0	
		self.N_k=np.zeros(self.model_num)
		

	def gaussian(self,x):
		return np.exp(-((x-self.mu)**2)/(2*self.sigma))/np.sqrt(2*math.pi*self.sigma)
		
		
	def update(self,pre_ganma):
		#x=[x1~xt]
		self.ganma = np.zeros((self.model_num,self.size))
	
		#calculate ganma
		for m in range(self.model_num):
			self.ganma[m] = pre_ganma*((self.c[m]*self.gaussian(self.x)[m]+function_b.upsilon)/np.sum(self.c*self.gaussian(self.x)+function_b.upsilon,axis=0))+function_b.upsilon
		
		for m in range(self.model_num):
			self.N_k[m] = np.sum(self.ganma[m])
		#calculate c
		for m in range(self.model_num):
			self.c[m]=self.N_k[m]/np.sum(self.N_k)
		#calculate sigma
		for m in range(self.model_num):
			self.sigma[m]=np.sum(self.ganma[m]*((self.x-self.mu[m])**2))/self.N_k[m]

		#calculate mu
		for m in range(self.model_num):
			self.mu[m]=np.sum(self.ganma[m]*self.x)/self.N_k[m]
		
	def output(self,x):
		return np.sum(self.c.reshape((1,-1))*self.gaussian(x).reshape((1,-1)))	


class_hmm_gmm.py


import numpy as np
import sys
import matplotlib.pyplot as plt
from class_b import function_b 
import copy
#construct
#s	=	state num
#x	=	data
#pi	=	init prob
#a 	=	transition prob
#c	=	gmm mixing coefficient
#mu	=	gmm mean
#sigma	=	gmm dispersion
#m	=	gmm model num

#method
#cal_alpha()	:	calculation alpha
#cal_beta()	:	calculation beta
#cal_exp()	:	calculation expected value(?)
#update()	:	update prams
#makestatelist	:	make statelist 
#viterbi	:	viterbi algorithm


class Hmm_Gmm():
	upsilon=1.0e-300
	def __init__(self,s,x,pi,a,c,mu,sigma,m):
		self.s = s
		self.x = x
		self.a = np.array(a)
		self.size=len(x)
		self.b=[]	
	
		for i in range(s):	
			self.b.append(function_b(c[i],mu[i],sigma[i],m,self.size,x))
		self.pi = np.array(pi)	
	
		self.alpha = np.zeros((self.s,self.size))
		self.beta = np.zeros((self.s,self.size))		
		self.prob = None
		self.kusai = np.zeros((self.s,self.s,self.size-1))
		self.ganma = np.zeros((self.s,self.size))
		self.preganma = None
		self.delta = None
		self.phi = None
		self.F = None
		self.P = None
		self.max_statelist=None

	def cal_alpha(self):
		
	
		for t in range(self.size):
			for j in range(self.s):
				if t==0:
					self.alpha[j,t]=self.pi[j]*self.b[j].output(self.x[t])
				else:
					self.alpha[j,t]=0
					for i in range(self.s): 
						self.alpha[j,t]+=self.alpha[i,t-1]*self.a[i,j]*self.b[j].output(self.x[t])
		self.prob =0
		for i in range(self.s):
			self.prob+=self.alpha[i,self.size-1]+Hmm_Gmm.upsilon
		
 
	def cal_beta(self):
		for t in reversed(range(self.size)):
			for i in range(self.s):
				
				if t ==self.size-1:
					self.beta[i,t]=1	
				else:
					self.beta[i,t]=0	
					for j in range(self.s):
						self.beta[i,t]+=self.beta[j,t+1]*self.a[i,j]*self.b[j].output(self.x[t+1])
	
	def cal_exp(self):
		for t in range(self.size-1):
			for i in range(self.s):
				for j in range(self.s):
								
					self.kusai[i,j,t] = ((self.alpha[i,t]*self.a[i,j]*self.b[j].output(self.x[t+1])*self.beta[j,t+1])+(Hmm_Gmm.upsilon/(self.s)))/self.prob
			
		for t in range(self.size):
			for i in range(self.s):
				#self.ganma[i,t] = (self.alpha[i,t]*self.beta[i,t])/self.prob
				#self.ganma[i,t] = (self.alpha[i,t]*self.beta[i,t])/self.prob + Hmm_Gmm.uplilon
				self.ganma[i,t] = ((self.alpha[i,t]*self.beta[i,t])+(Hmm_Gmm.upsilon))/self.prob
	def update(self):
		#cal exp

		self.cal_alpha()
		self.cal_beta()
		self.cal_exp() 
		#cal pi
		self.pi = self.ganma[:,0]
			
		#cal a
		for i in range(self.s):
			for j in range(self.s):
				
				#self.a[i,j]=np.sum(self.ganma[i,:])

				self.a[i,j]=np.sum(self.kusai[i,j,:])/np.sum(self.ganma[i,:])
				#self.a[i,j]=np.sum(self.kusai[i,j,:])

	
		#cal b
		#cal preganma
		self.preganma = np.zeros((self.s,self.size))	
		for t in range(self.size):	
			temp = self.alpha[:,t]*self.beta[:,t]+self.upsilon
			self.preganma[:,t] = temp/np.sum(temp)
		for i in range(self.s):
			self.b[i].update(np.array(self.preganma[i]))	
	
	def makestatelist(self,f,size):
		statelist = np.zeros(size)
		for t in reversed(range(size)):
			if t == size-1:
				statelist[t] = f
			else:	
				statelist[t] = self.phi[statelist[t+1],t+1]
		return statelist

	def viterbi(self,x):
		data=np.array(x)
		self.delta = np.zeros((self.s,data.shape[0]))
		self.phi = np.zeros((self.s,data.shape[0]))
		
		for t in range(data.shape[0]):
			for i in range(self.s):	
				if t == 0:
					self.phi[i,t] = np.nan
					self.delta[i,t] = np.log(self.pi[i])	
				else: 
					self.phi[i,t] = np.argmax(self.delta[:,t-1]+np.log(self.a[:,i]))					
					self.delta[i,t] = np.max(self.delta[:,t-1]+np.log(self.a[:,i]*self.b[i].output(data[t])+self.upsilon))	
		self.F = np.argmax(self.delta[:,data.shape[0]-1])
		self.P = np.max(self.delta[:,data.shape[0]-1])	
		#print(self.makestatelist(self.F,data.shape[0]))
		self.max_statelist=self.makestatelist(self.F,data.shape[0])
		

ex_hmm_demo.py



import numpy as np
import matplotlib.pyplot as plt
import sys 
from class_hmm_gmm import Hmm_Gmm 

model_num = 1
state_num = 4
color=['#ff0000','#00ff00','#0000ff','#888888']
value_size = int(10*model_num*state_num)


create_mu = []
create_sigma= []
create_pi =[]
for i in range(state_num):
	temp = []
	temp1 = []
	temp2=[]
	for j in range(model_num):
		temp.extend([(np.random.rand()*30-np.random.rand()*30)])
		temp1.extend([np.random.rand()+1])
		temp2.extend([np.random.rand()])
		
	create_mu.append(temp)
	create_sigma.append(temp1)	
	create_pi.append((temp2/np.sum(temp2)).tolist())

#create data
values = []
for i in range(state_num):
	temp = []
	for it in range(int(value_size/state_num/model_num)):
		for j in range(model_num):
			temp.extend(np.random.normal(create_mu[i][j],create_sigma[i][j],1).tolist())
	values.append(temp)

x_plot = np.linspace(0,value_size,value_size)
for i in range(state_num):
	plt.scatter(x_plot[i*len(values[0]):i*len(values[0])+len(values[0])],values[i],c=color[i],marker='+')

for i in range(state_num):
	for j in range(model_num):
		plt.plot(x_plot[i*len(values[0]):i*len(values[0])+len(values[0])],np.full(len(values[0]),create_mu[i][j]),c=color[i])
	

plt.show()


train_values=[]
for i in range(state_num):
	train_values.extend(values[i])

train_mu = []
train_sigma= []
train_pi =[]
for i in range(state_num):
	temp = []
	temp1 = []
	temp2=[]
	for j in range(model_num):
		temp.extend([np.random.rand()*30-np.random.rand()*30])
		temp1.extend([np.random.rand()+1])
		temp2.extend([np.random.rand()])
		
	train_mu.append(temp)
	train_sigma.append(temp1)	
	train_pi.append((temp2/np.sum(temp2)).tolist())


##backup train
b_train_mu = train_mu
b_train_sigma = train_sigma

plt.scatter(x_plot,train_values,c="#000000",marker='+')

for i in range(state_num):
	for j in range(model_num):
		plt.plot(x_plot[i*len(values[0]):i*len(values[0])+len(values[0])],np.full(len(values[0]),train_mu[i][j]),c=color[i])
	
plt.show()
#init a
a = []

for j in range(state_num):
	temp = []
	for i in range(state_num):	
		temp.append(np.random.rand())	
		
	a.append(temp)

for i in range(state_num):
	for j in range(state_num):	
		#a[j][i] = a[j][i]/np.sum(np.array(a)[:,i])
		a = (np.array(a)/np.sum(np.array(a),axis=0)).tolist()

#init init_prob
init_prob = []

for j in range(state_num):
	temp2.extend([np.random.rand()])
#temp2[0]+=1
init_prob.extend((temp2/np.sum(temp2)).tolist())
#init Hmm_Gmm model
hmm_gmm = Hmm_Gmm(state_num,train_values,init_prob,a,train_pi,train_mu,train_sigma,model_num)

#train roop
print(create_mu)
print(a)
count = 0
while(count != 100):	
	hmm_gmm.update()	
	count +=1
	
#viterbi
hmm_gmm.viterbi(train_values)
print(hmm_gmm.max_statelist)
#*************************************
plt.scatter(x_plot,train_values,c="#000000",marker='+')

##plot of init train mu 	

for i in range(state_num):
	for j in range(model_num):
		plt.plot(x_plot[i*len(values[0]):i*len(values[0])+len(values[0])],np.full(len(values[0]),b_train_mu[i][j]),c="#dddddd")

for i in range(state_num):
	for j in range(model_num):
		plt.plot(x_plot[i*len(values[0]):i*len(values[0])+len(values[0])],np.full(len(values[0]),hmm_gmm.b[i].mu[j]),c=color[i])
	
plt.show()


résultat スクリーンショット 2020-09-19 2.45.03.png Les points en sortie de la valeur initiale $ \ mu $ et de la distribution gaussienne. Cette fois, quatre états sont supposés et des points sont générés à partir de chacune des quatre distributions gaussiennes. Depuis le début, chaque état est appelé état 1 (rouge), état 2 (vert), état 3 (bleu) et état 4 (gris).

スクリーンショット 2020-09-19 2.45.12.png

La valeur moyenne $ \ mu $ a été réinitialisée avec un nombre aléatoire. Nous continuerons d'en apprendre davantage à ce sujet. C'est le résultat d'une mise à jour 100 fois. スクリーンショット 2020-09-19 2.45.29.png

C'est un résultat d'apprentissage. Les états sont dans l'ordre de l'état 2, état 1, état 4 et état 3 depuis le début, mais en conséquence, je pense que cela fonctionne bien.

Finalement##

Cette fois, j'ai essayé de résumer le modèle Markov caché à ma manière, mais le montant est plus important que ce à quoi je m'attendais, et il y a de nombreux endroits où j'oublie comment penser même si je l'ai implémenté il y a quelques mois, il n'est donc pas possible de sortir tôt. J'ai vraiment senti que c'était important. En partie à cause de cela, je n'ai pas été en mesure d'expliquer l'algorithme bitterbi ou la fonction de la liste d'états. Je le résumerai à une date ultérieure. ~~ (Il y a aussi une raison pour laquelle je n'ai pas pu entrer à cause de mon manque de capacité d'écriture ...) ~~ De plus, la distribution gaussienne peut ne pas s'intégrer correctement dans ce programme, je voudrais donc examiner un peu plus si elle est due à mon propre programme ou si elle est tombée dans une solution localement optimale. J'ai également essayé d'implémenter GMM-HMM (changer model_num le rend GMM), mais cela n'a pas bien fonctionné, et où cela s'est-il encore mal passé (ou est-il à nouveau tombé dans la solution optimale locale)? Je voudrais le revoir. Je suis désolé que cela ne soit pas correctement terminé De plus, j'ai trouvé un article qui semble décrire une méthode pour résoudre un problème qui s'inscrit dans une solution localement optimale, je vais donc l'étudier également.

[référence]##

https://mieruca-ai.com/ai/markov_model_hmm/ [Explication technique] Modèle de Markov et modèle de Markov caché https://qiita.com/kenmatsu4/items/59ea3e5dfa3d4c161efb Explication approfondie de l'algorithme EM

Recommended Posts

Implémentation Python du modèle Markov caché continu
PRML Chapter 13 Estimation la plus probable Implémentation Python du modèle de Markov caché
Implémentation du modèle de sujet d'espace continu
Implémentation d'estimation la plus probable du modèle de sujet en python
Implémentation Python du filtre à particules
Implémentation d'estimation bayésienne de variante du modèle de sujet en python
Implémentation du tri rapide en Python
Implémentation du jeu de vie en Python
Implémentation des notifications de bureau à l'aide de Python
Implémentation Python de l'arborescence de segments non récursive
Implémentation de Light CNN (Python Keras)
Implémentation du tri original en Python
Implémentation de la méthode Dyxtra par python
Explication du modèle d'optimisation de la production par Python
PRML Chapitre 14 Implémentation Python de modèle mixte conditionnel
Pourquoi l'implémentation Python d'ISUCON 5 a utilisé Bottle
Les bases de Python ①
PRML Chapitre 11 Implémentation Python Monte Carlo Chaîne de Markov
Bases de python ①
Copie de python
[Entretien de codage] Implémentation de la machine cryptographique Enigma (Python)
Implémentation du modèle Deep Learning pour la reconnaissance d'images
Explication de la distance d'édition et de l'implémentation en Python
Introduction de Python
Exemple d'implémentation d'un système de traitement LISP simple (version Python)
Comparaison de l'implémentation de plusieurs moyennes mobiles exponentielles (DEMA, TEMA) en Python
[# 2] Créez Minecraft avec Python. ~ Dessin du modèle et implémentation du lecteur ~
Implémentation python de la classe de régression linéaire bayésienne
Un mémorandum sur la mise en œuvre des recommandations en Python
[Python] Opération d'énumération
Liste des modules python
Implémentation RNN en python
Implémentation ValueObject en Python
Unification de l'environnement Python
[python] comportement d'argmax
Utilisation des locaux Python ()
Installation de Python 3.3 rc1
# 4 [python] Bases des fonctions
Connaissance de base de Python
Résumé des arguments Python
Bases de python: sortie
Installation de matplotlib (Python 3.3.2)
Application de Python 3 vars
Implémentation SVM en python
Divers traitements de Python
Implémentation Python du mode de fusion CSS3 et discussion sur l'espace colorimétrique
Une implémentation Python simple de la méthode k-voisinage (k-NN)
Comment écrire un exemple d'implémentation E14 Python en temps réel hors ligne
[Avec une explication simple] Implémentation Scratch d'une machine Boltsman profonde avec Python ②
[Avec une explication simple] Implémentation Scratch d'une machine Boltzmann profonde avec Python ①