Python-Implementierung eines kontinuierlichen Hidden-Markov-Modells
Da ich dieses Mal angefangen habe, Qiita zu veröffentlichen, werde ich über die Implementierung des kontinuierlichen Hidden-Markov-Modells durch Python schreiben, das ich vor einigen Monaten nur schwer implementieren konnte. Als Memorandum werde ich aus den Grundlagen in der Reihenfolge schreiben, in der ich studiert habe. Ich denke auch, dass dieses Ergebnis wahrscheinlich erreicht wird, aber es gibt einige verdächtige Teile. Wenn Sie also etwas bemerken, würde ich es begrüßen, wenn Sie es mich wissen lassen könnten.
Was ist das Markov-Modell?
Das Markov-Modell ist ein Modell, das auf der Grundlage des Markov-Prozesses erstellt wurde, und der Markov-Prozess ist ein stochastischer Prozess mit Markov-Eigenschaften.
Markov kam viel heraus und ich konnte es nicht leicht verstehen, also werde ich mein eigenes Verständnis erklären.
Betrachten Sie zunächst ein Diagramm wie die folgende Graphentheorie.
In dieser Figur wird angenommen, dass sich der Zustand nach einer bestimmten Zeit wahrscheinlich in der Richtung ändert, in die der Pfeil zeigt, und dass die Wahrscheinlichkeit die dem Pfeil zugeordnete Zahl ist. Wenn es jetzt Zustand A ist, geht es mit einer Wahrscheinlichkeit von 0,3 in Zustand B und mit einer Wahrscheinlichkeit von 0,7 in Zustand C über. Wenn der A-Zustand von X oder Y kommt, ändert sich die Übergangswahrscheinlichkeit von Zustand A zu Zustand B oder Zustand C nicht. Mit anderen Worten, die Eigenschaft, dass die Übergangswahrscheinlichkeit zum nächsten Zustand (B oder C) nur vom aktuellen Zustand (A) abhängt, wird als Markov-Eigenschaft bezeichnet. Der stochastische Prozess wird als "eine stochastische Variable, die sich mit der Zeit ändert.-Wikipedia-" geschrieben. In diesem Fall ändert sich die Wahrscheinlichkeit, einen Zustand anzunehmen, mit der Zeit, so dass von einem stochastischen Prozess gesprochen werden kann. ~~ Entschuldigung, ich werde die strenge Definition anderen Informationsseiten überlassen ~~ Ein auf diese Weise erstelltes Modell kann als Markov-Modell bezeichnet werden.
Verstecktes Markov Modell
Das versteckte Markov-Modell ist ein Modell wie das zuvor erwähnte Markov-Modell, es wird jedoch verwendet, wenn es je nach Status einen bestimmten Wert (Ausgabe) hat, anstatt den Status nicht zu kennen. Lassen Sie uns ein konkretes Beispiel machen, das auf der Verarbeitung natürlicher Sprache basiert. Betrachten Sie zunächst den Beispielsatz (Ausgabe) "Ich buche einen Film" (ich reserviere einen Film). Zu diesem Zeitpunkt sind "Ich" und "Film" Nomenklatur, "Buch" ist ein Verb und "a" ist ein Akronym. Solche Teilwörter gelten als verborgen, da sie aus diesem Satz nicht bekannt sind. Angenommen, Sie möchten, dass die Maschine den Teil dieses Satzes entscheidet. Obwohl "Ich" und "A" fast eindeutig bestimmt sind, können "Buch" und "Film" auch als Verben verwendet werden, z. B. als Buch oder Film, was eine Nomenklatur ist. Überlegen wir uns daher, ob es möglich ist, Teile mithilfe des Hidden-Markov-Modells zu klassifizieren.
Im Englischen wird die Reihenfolge der Teile bis zu einem gewissen Grad festgelegt, z. B. 5 Satzmuster. Daher ist es wahrscheinlich, dass das Ändern von Teilwörtern einfach ist, wie in der folgenden Tabelle gezeigt. (Der Preis wird nach Ihrem eigenen Gefühl festgelegt.)
| Satzanfang | Substantiv | Verb | Artikel | Ende des Satzes | |
|---|---|---|---|---|---|
| Satzanfang | 0 | 0.6 | 0.2 | 0.2 | 0 |
| Substantiv | 0 | 0.3 | 0.5 | 0.1 | 0.1 |
| Verb | 0 | 0.6 | 0 | 0.3 | 0.1 |
| Artikel | 0 | 1 | 0 | 0 | 0 |
| Ende des Satzes | 0 | 0 | 0 | 0 | 0 |
Stellen Sie sich die Linke als den aktuellen Zustand und die Spitze als den nächsten Zustand vor. Zum Beispiel beträgt die Wahrscheinlichkeit, dass eine Nase neben eine Nase kommt, 0,3, die Wahrscheinlichkeit, dass ein Verb kommt, 0,5, die Wahrscheinlichkeit, dass ein Akronym kommt, 0,1 und die Wahrscheinlichkeit, dass das Ende eines Satzes kommt, ist 0,1. ~~ Die Verwendung von 0 für die Wahrscheinlichkeit ist nicht gut für das Nullfrequenzproblem ~~
Dann wird angenommen, dass die Ausgabe in jedem Zustand mit der folgenden Wahrscheinlichkeit genommen wird.
Mit anderen Worten, wenn es sich um eine Nomenklatur handelt, ist "I" 0,4, "a" ist 0,0, "Buch" ist 0,3 und "Film" ist 0,2. Das Verfahren des Zustandsübergangs zu diesem Zeitpunkt ist wie folgt. (Der Zustand, in dem die Übergangswahrscheinlichkeit 0 ist, wird im Voraus gelöscht.)
Dann ist die Anzahl der Wörter n, die Wortkette (beobachtbar) ist $ W = w_1, w_2, ..., w_n $ und die Teilzeichenfolge (nicht beobachtbar: Zustand) ist $ T = t_1, t_2 ,. .., t_n $ ist die Wahrscheinlichkeit, die dieses Modell bei der Entscheidung für eine Route angibt
Kann vertreten werden durch. Mit anderen Worten, die Wahrscheinlichkeit, vom vorherigen in den aktuellen Zustand zu wechseln
(Im Beispiel des diesmal verwendeten Satzes ist jedoch $ w_0 $ leer, $ w_1 $ ist 'I', ..., $ w_n $ ist leer.) Dann kann der Teil dieses Satzes vorhergesagt werden, indem die Übergangsmethode mit der höchsten Wahrscheinlichkeit gefunden wird. Übrigens kann im Fall dieses Beispiels gesagt werden, dass es auf leicht verständliche Weise voreingenommen ist. ..
Beiseite Die anfängliche Wahrscheinlichkeit bleibt weiterhin $ \ pi $ Die Übergangswahrscheinlichkeit ist $ a_ {i, j} $ und die Erscheinungswahrscheinlichkeit ist $ b_j (w) $.
Kontinuierliches verstecktes Markov-Modell
Da das zuvor erwähnte versteckte Markov-Modell über vier Wörter nachdachte, kann die Ausgabe als eine Ausgabe mit vier Werten betrachtet werden. Die Ausgabe war also ein diskreter Wert. Daher kann man sagen, dass es sich um ein diskretes verstecktes Markov-Modell handelt. Daher wird das Hidden-Markov-Modell bei Verwendung kontinuierlicher Werte wie Temperatur als kontinuierliches Hidden-Markov-Modell bezeichnet. In dem zuvor erwähnten diskreten Markov-Modell wird die Erscheinungswahrscheinlichkeit jedes Zustands einzeln berechnet. Da es jedoch nicht realistisch ist, alle Werte kontinuierlichen Werten zuzuweisen, definieren wir die Erscheinungswahrscheinlichkeit unter Verwendung einer gemischten Gaußschen Verteilung. .. Bevor wir uns mit der Implementierung des kontinuierlichen Hidden-Markov-Modells befassen, schreiben wir zunächst über die gemischte Gaußsche Verteilung.
Gaußsche Verteilung / gemischte Gaußsche Verteilung
Die gemischte Gaußsche Verteilung ist eine Darstellung der Wahrscheinlichkeitsverteilung durch Hinzufügen einer Anzahl von Gaußschen Verteilungen.
Erstens wird die Gaußsche Verteilung auch als Normalverteilung bezeichnet.
Es wird vertreten durch.
-Bild aus Wikipedia ausgeliehen- Wie Sie in diesem Diagramm sehen können, ist der Wert umso größer, je näher dieses x am Durchschnitt von $ \ mu $ (0 in dieser Abbildung) liegt. Und wenn die Dispersion groß ist, ist sie flach. Mit anderen Worten, wenn dieser Wert als Wahrscheinlichkeit betrachtet wird, ist der Wert umso höher, je näher der Wert von x am Durchschnitt $ \ mu $ liegt, und wenn die Varianz groß ist, kann dieser Wert auch dann angenommen werden, wenn er weit vom Durchschnitt entfernt ist. Ich kann sagen.
Außerdem ist die gemischte Gaußsche Verteilung wie das Addieren mehrerer Normalverteilungen, und wenn K Gaußsche Verteilungen kombiniert werden, lautet die Formel für die gemischte Gaußsche Verteilung
Es gibt jedoch eine Bedingung von $ \ sum_ {k = 1} ^ K \ pi_k = 1 $.
Es wird vertreten durch. Wenn Sie versuchen, dies zu veranschaulichen
Es sieht aus wie das.
Clustering mit gemischter Gaußscher Verteilung
In Zukunft ist es bei der Implementierung des kontinuierlichen Hidden-Markov-Modells erforderlich, die Gaußsche Verteilung anhand der Datenpunkte zu identifizieren. Daher möchte ich den Fluss sehen, indem ich Soft Clustering mit einer gemischten Gaußschen Verteilung implementiere.
Zeichnen wir zunächst die Punkte basierend auf den drei Gaußschen Verteilungen.
Dies ist die Handlung. Diese Punkte befinden sich in einem Zustand, in dem bekannt ist, aus welcher Gaußschen Verteilung sie erzeugt wurden. (Vollständige Daten) Beim eigentlichen Clustering weiß ich jedoch nicht, aus welcher Distribution es generiert wurde. (Unvollständige Daten) Daher werden wir von nun an die Parameter der gemischten Gaußschen Verteilung, die aus den drei Gaußschen Verteilungen besteht, so anpassen, dass sie gemäß der Gaußschen Verteilung am besten in drei geteilt werden können.
Wie können wir es in drei Vertiefungen aufteilen? Dies bedeutet, dass die gleichzeitige Wahrscheinlichkeit aller Punkte gemäß der gewünschten gemischten Gaußschen Verteilung am höchsten sein sollte. Mit anderen Worten
$ X = [x_1, x_2, ..., x_n] $ ・ ・ ・ Jeder Datenpunkt $ K $ ・ ・ ・ Anzahl der Modelle (Anzahl der Cluster)
Stellen Sie die Parameter so ein, dass sie am höchsten sind. Dies wird auch als logarithmische Wahrscheinlichkeitsfunktion bezeichnet. Ich möchte die Parameter ändern, um diese Funktion zu maximieren, aber da die Summe der Summe erscheint und es schwierig ist, sie analytisch zu lösen, werde ich sie mit dem EM-Algorithmus lösen.
Als Methode des EM-Algorithmus E Schritt: Berechnung des erwarteten Wertes M Schritt: Passen Sie die Parameter an, um den erwarteten Wert zu maximieren Es ist eine Methode, um durch Wiederholen die optimale Lösung zu finden.
Der Fluss ist diesmal
- Bestimmen Sie den Anfangswert von $ \ pi, \ mu, \ sigma $
- Berechnen Sie die Belastungsrate (E-Schritt)
- Passen Sie die Werte der Parameter $ \ pi, \ mu, \ sigma $ an.
- Wenn es keinen Unterschied in der Wahrscheinlichkeit gibt, endet der Prozess. Wenn es sich immer noch ändert, gehen Sie zurück zu 2.
Ich werde es im Fluss tun.
Ursprünglicher Wert##
Der anfängliche Anfangswert wird entsprechend festgelegt. (Wenn Sie diesen Wert irgendwie kennen, scheint er natürlich schneller zu konvergieren, wenn Sie ihn anpassen.)
E Schritt
Als nächstes folgt der E-Schritt. Hier kommt der Begriff Belastungsrate auf, der die Wahrscheinlichkeit ist, dass ein Punkt x aus dem Modell k auftritt. Angesichts der obigen Darstellung wusste ich nicht, von welchem Gaußschen Verteilungsmodell stammt. Betrachten Sie die latente Variable z, die angibt, ob diese aus einem Modell dafür geboren wurde. Wenn $ z_ {i, k} = 1 $ ist, dann ist der Punkt $ x_i $ ein Punkt, der aus dem Modell k geboren wurde.
Wenn dies zu einem Ausdruck gemacht wird
Ebenfalls,
weil
Ist. (Bayes-Theorem)
Da $ p (x_i) $ ein in Bezug auf z marginalisierter Wert ist, wird er durch die Formel der gemischten Gaußschen Verteilung ausgedrückt.
Es wird sein.
Hier,
Ich habe nicht verstanden, warum $ p (z_ {i, k} = 1) = \ pi_k $, aber dies ist ein Wert, der nichts mit x zu tun hat und wie viel er zum Modell k gehört. Es ist ein Gefühl, dass es sicherlich der Fall sein kann, wenn es eine Menge ist, um zu vergleichen, ob es einfach ist oder nicht.)
Der Grund, warum dies als Belastungsrate bezeichnet wird, ist, dass es zeigt, wie viel das Modell k den Anteil der gemischten Gaußschen Verteilung einnimmt. Wenn Sie sich die Figur vor zwei Jahren ansehen, können Sie sich das irgendwie vorstellen.
M Schritt
Als nächstes folgt die Aktualisierung der Parameter. Für jede Richtlinie ist für $ \ mu und \ sigma $ die Funktion mit der höchsten Wahrscheinlichkeit ausreichend, sodass wir mit der Richtlinie fortfahren, den Punkt zu differenzieren und nach dem Punkt zu suchen, an dem der Wert 0 wird. Da $ \ pi $ die Bedingung $ \ sum_ {k = 1} ^ K \ pi_k = 1 $ hat, wird es auch durch Lagranges unbestimmte Multiplikatormethode gelöst.
Die Ableitung dieser Unterscheidung ist gründliche Erklärung des EM-Algorithmus, daher erklären andere Leute sie auf sehr leicht verständliche Weise, daher werde ich nur die Ergebnisse zusammenfassen. ~~ Ich persönlich habe es abgeleitet, aber es war sehr schwierig ~~
Ergebnis
Das Ergebnis ist.
Implementierung einer gemischten Gaußschen Verteilung
Da die Ableitung zusammengefasst wurde, habe ich tatsächlich das Clustering der gemischten Gaußschen Verteilung implementiert.
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()
Es sollte so funktionieren. In diesem Programm werden alle Parameter zufällig ausgegeben, sodass möglicherweise keine guten Cluster erstellt werden können. ~~ Wenn solche Daten eingehen, können sie nicht geclustert werden ... ~~
Implementierung des kontinuierlichen Hidden-Markov-Modells
Es war ein Umweg, aber ich würde von nun an gerne die Implementierung des kontinuierlichen Markov-Modells sehen. Wenn die Beobachtungsinformation y = {$ y_1, y_2, ..., y_t $}, die ein kontinuierliches Markov-Modell ist, erhalten wird, wird die Wahrscheinlichkeit $ p (y | M) $, dass dieses Modell diesen Wert zurückruft, durch den Gitteralgorithmus bestimmt. tun können. Der Gitteralgorithmus $ p (y | M) = \ sum_ {i} \ alpha (i, T) $ (Dies ist jedoch der Endzustand)
Wird benötigt bei. Und dieses $ \ alpha (i, t) $ ist
\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}
ist. Und $ a_ {j, i} $ und $ b_ {j, i} (y_t) $, die hier erscheinen, sind vor einiger Zeit erschienen, aber sie zeigen die Übergangswahrscheinlichkeit und Erscheinungswahrscheinlichkeit, die auch im versteckten Markov-Modell erschienen sind. Ich bin. (Diese Geschichte hatte nichts mit Lernen zu tun ...) Dieses $ \ alpha (i, t) $ verfolgt die Wahrscheinlichkeit von vorne (Vorwärtsfunktion), definiert jedoch eine Funktion, die dies vom Ende aus verfolgt (Rückwärtsfunktion).
\begin{eqnarray}
\beta_{i,T}=\left\{ \begin{array}{ll}
\ 1 & (Endzustand) \\
\ 0 & (Nicht im Endzustand) \\
\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)
Die Modellparameter des Hidden-Markov-Modells können mit dem Baum-Welch-Algorithmus gelernt werden, bei dem es sich um eine Art EM-Algorithmus handelt. Der Baumwelch-Algorithmus wird verwendet, um in E-Schritten zu berechnen. $ \ Gamma (i, j, t) = \ frac {\ alpha (i, t-1) ・ a_ {i, j} ・ b_ {j} beta \ beta (j, t)} {p (y | M. )} $
Dieser Wert. Dieser Wert wird als die Wahrscheinlichkeit angesehen, zum Zeitpunkt t vom Zustand i zu j überzugehen und $ y_t $ auszugeben, geteilt durch die Wahrscheinlichkeit, dass y in diesem Modell auftritt.
Und die Parameterschätzungsformel (Aktualisierungsformel) lautet
Nächster, Weil b eine Gaußsche Verteilung hat
Wird benötigt bei.
Dann ist es ein Code. Eine Implementierungsdatei der Gaußschen Distribution, da sie nacheinander angepasst wurde. Es ist in drei Teile unterteilt: eine versteckte Markov-Modelldatei und eine ausführbare Datei. (Zuerst hatte ich vor, es mit GMM-HMM zu machen, daher gibt es viele Teile, deren Implementierung nutzlos ist. Bitte warten Sie eine Weile, da ich es erneut aktualisieren werde, sobald GMM-HMM funktioniert.)
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()
Ergebnis
Die Punkte, die aus dem anfänglichen $ \ mu $ -Wert und der Gaußschen Verteilung ausgegeben werden. Dieses Mal werden vier Zustände angenommen und Punkte werden aus jeder der vier Gaußschen Verteilungen erzeugt. Jeder Zustand wird von Anfang an als Zustand 1 (rot), Zustand 2 (grün), Zustand 3 (blau) und Zustand 4 (grau) bezeichnet.
Der Durchschnittswert $ \ mu $ wurde mit einer Zufallszahl zurückgesetzt.
Wir werden weiter darüber lernen. Dies ist das Ergebnis einer 100-fachen Aktualisierung.
Es ist ein Lernergebnis. Die Zustände liegen von Anfang an in der Reihenfolge von Zustand 2, Zustand 1, Zustand 4 und Zustand 3, aber als Ergebnis habe ich das Gefühl, dass es gut funktioniert.
Schließlich##
Dieses Mal habe ich versucht, das versteckte Markov-Modell auf meine eigene Weise zusammenzufassen, aber es gibt mehr als ich erwartet hatte, und es gibt viele Stellen, an denen ich vergesse zu denken, obwohl ich es vor einigen Monaten implementiert habe, sodass es nicht möglich ist, vorzeitig auszugeben Ich fand es wirklich wichtig. Aufgrund dieses Einflusses kann ich den Bitterbi-Algorithmus und die Funktion der Statusliste nicht erklären. Ich werde es zu einem späteren Zeitpunkt noch einmal zusammenfassen. ~~ (Es gibt auch einen Grund, warum ich wegen mangelnder Schreibfähigkeit nicht eintreten konnte ...) ~~ Außerdem passt die Gaußsche Verteilung möglicherweise nicht richtig in dieses Programm, daher möchte ich etwas genauer überlegen, ob dies auf mein eigenes Programm zurückzuführen ist oder ob es in eine lokal optimale Lösung gefallen ist. Ich habe auch versucht, GMM-HMM zu implementieren (das Ändern von model_num macht es zu GMM), aber es hat nicht gut funktioniert, und wo ist es wieder schief gelaufen (oder ist es wieder in die lokale optimale Lösung gefallen)? Ich würde es gerne überprüfen. Es tut mir leid, dass es nicht richtig abgeschlossen wurde. Außerdem habe ich ein Papier gefunden, das eine Methode zur Lösung eines Problems zu beschreiben scheint, das in eine lokal optimale Lösung fällt, daher werde ich dies ebenfalls untersuchen.
[Referenz]##
https://mieruca-ai.com/ai/markov_model_hmm/ [Technische Erklärung] Markov-Modell und verstecktes Markov-Modell https://qiita.com/kenmatsu4/items/59ea3e5dfa3d4c161efb Ausführliche Erläuterung des EM-Algorithmus
Recommended Posts