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.
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.
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)
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
Peut être représenté par. En d'autres termes, la probabilité de passer de l'état précédent à l'état actuel
(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.
Il est représenté par.
-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
Cependant, il existe une condition de $ \ sum_ {k = 1} ^ K \ pi_k = 1 $.
Il est représenté par. Si vous essayez d'illustrer cela
Ç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.
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
$ 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
- Déterminez la valeur initiale de $ \ pi, \ mu, \ sigma $
- Calculez le taux de charge (étape E)
- Ajustez les valeurs des paramètres $ \ pi, \ mu, \ sigma $.
- 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
Aussi,
Parce que
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.
Ce sera.
ici,
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
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
Suivant, Parce que b a une distribution gaussienne
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
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).
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.
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