# probHMM.py - Hidden Markov Model # AIFCA Python3 code Version 0.7.6 Documentation at http://aipython.org # Artificial Intelligence: Foundations of Computational Agents # http://artint.info # Copyright David L Poole and Alan K Mackworth 2017. # This work is licensed under a Creative Commons # Attribution-NonCommercial-ShareAlike 4.0 International License. # See: http://creativecommons.org/licenses/by-nc-sa/4.0/deed.en import random from probStochSim import sample_one, sample_multiple class HMM(object): def __init__(self, states, obsvars,pobs,trans,indist): """A hidden Markov model. states - set of states obsvars - set of observation variables pobs - probability of observations, pobs[i][s] is P(Obs_i=True | State=s) trans - transition probability - trans[i][j] gives P(State=j | State=i) indist - initial distribution - indist[s] is P(State_0 = s) """ self.states = states self.obsvars = obsvars self.pobs = pobs self.trans = trans self.indist = indist # state # 0=middle, 1,2,3 are corners states1 = {'middle', 'c1', 'c2', 'c3'} # states obs1 = {'m1','m2','m3'} # microphones # pobs gives the observation model: #pobs[mi][state] is P(mi=on | state) closeMic=0.6; farMic=0.1; midMic=0.4 pobs1 = {'m1':{'middle':midMic, 'c1':closeMic, 'c2':farMic, 'c3':farMic}, # mic 1 'm2':{'middle':midMic, 'c1':farMic, 'c2':closeMic, 'c3':farMic}, # mic 2 'm3':{'middle':midMic, 'c1':farMic, 'c2':farMic, 'c3':closeMic}} # mic 3 # trans specifies the dynamics # trans[i] is the distribution over states resulting from state i # trans[i][j] gives P(S=j | S=i) sm=0.7; mmc=0.1 # transition probabilities when in middle sc=0.8; mcm=0.1; mcc=0.05 # transition probabilities when in a corner trans1 = {'middle':{'middle':sm, 'c1':mmc, 'c2':mmc, 'c3':mmc}, # was in middle 'c1':{'middle':mcm, 'c1':sc, 'c2':mcc, 'c3':mcc}, # was in corner 1 'c2':{'middle':mcm, 'c1':mcc, 'c2':sc, 'c3':mcc}, # was in corner 2 'c3':{'middle':mcm, 'c1':mcc, 'c2':mcc, 'c3':sc}} # was in corner 3 # initially we have a uniform distribution over the animal's state indist1 = {st:1.0/len(states1) for st in states1} hmm1 = HMM(states1, obs1, pobs1, trans1, indist1) from display import Displayable class HMM_VE_filter(Displayable): def __init__(self,hmm): self.hmm = hmm self.state_dist = hmm.indist def filter(self, obsseq): """updates and returns the state distribution following the sequence of observations in obsseq using variable elimination. Note that it first advances time. This is what is required if it is called sequentially. If that is not what is wanted initially, do an observe first. """ for obs in obsseq: self.advance() # advance time self.observe(obs) # observe return self.state_dist def observe(self, obs): """updates state conditioned on observations. obs is a list of values for each observation variable""" for i in self.hmm.obsvars: self.state_dist = {st:self.state_dist[st]*(self.hmm.pobs[i][st] if obs[i] else (1-self.hmm.pobs[i][st])) for st in self.hmm.states} norm = sum(self.state_dist.values()) # normalizing constant self.state_dist = {st:self.state_dist[st]/norm for st in self.hmm.states} self.display(2,"After observing",obs,"state distribution:",self.state_dist) def advance(self): """advance to the next time""" nextstate = {st:0.0 for st in self.hmm.states} # distribution over next states for j in self.hmm.states: # j ranges over next states for i in self.hmm.states: # i ranges over previous states nextstate[j] += self.hmm.trans[i][j]*self.state_dist[i] self.state_dist = nextstate hmm1f1 = HMM_VE_filter(hmm1) # hmm1f1.filter([{'m1':0, 'm2':1, 'm3':1}, {'m1':1, 'm2':0, 'm3':1}]) ## HMM_VE_filter.max_display_level = 2 # show more detail in displaying # hmm1f2 = HMM_VE_filter(hmm1) # hmm1f2.filter([{'m1':1, 'm2':0, 'm3':0}, {'m1':0, 'm2':1, 'm3':0}, {'m1':1, 'm2':0, 'm3':0}, # {'m1':0, 'm2':0, 'm3':0}, {'m1':0, 'm2':0, 'm3':0}, {'m1':0, 'm2':0, 'm3':0}, # {'m1':0, 'm2':0, 'm3':0}, {'m1':0, 'm2':0, 'm3':1}, {'m1':0, 'm2':0, 'm3':1}, # {'m1':0, 'm2':0, 'm3':1}]) # hmm1f3 = HMM_VE_filter(hmm1) # hmm1f3.filter([{'m1':1, 'm2':0, 'm3':0}, {'m1':0, 'm2':0, 'm3':0}, {'m1':1, 'm2':0, 'm3':0}, {'m1':1, 'm2':0, 'm3':1}]) # How do the following differ in the resulting state distribution? # Note they start the same, but have different initial observations. ## HMM_VE_filter.max_display_level = 1 # show less detail in displaying # for i in range(100): hmm1f1.advance() # hmm1f1.state_dist # for i in range(100): hmm1f3.advance() # hmm1f3.state_dist from display import Displayable from probStochSim import resample class HMM_particle_filter(Displayable): def __init__(self,hmm,number_particles=1000): self.hmm = hmm self.particles = [sample_one(hmm.indist) for i in range(number_particles)] self.weights = [1 for i in range(number_particles)] def filter(self, obsseq): """returns the state distribution following the sequence of observations in obsseq using particle filtering. Note that it first advances time. This is what is required if it is called after previous filtering. If that is not what is wanted initially, do an observe first. """ for obs in obsseq: self.advance() # advance time self.observe(obs) # observe self.resample_particles() self.display(2,"After observing", str(obs), "state distribution:", self.histogram(self.particles)) self.display(1,"Final state distribution:", self.histogram(self.particles)) return self.histogram(self.particles) def advance(self): """advance to the next time. This assumes that all of the weights are 1.""" self.particles = [sample_one(self.hmm.trans[st]) for st in self.particles] def observe(self, obs): """reweight the particles to incorporate observations obs""" for i in range(len(self.particles)): for obv in obs: if obs[obv]: self.weights[i] *= self.hmm.pobs[obv][self.particles[i]] else: self.weights[i] *= 1-self.hmm.pobs[obv][self.particles[i]] def histogram(self, particles): """returns list of the probability of each state as represented by the particles""" tot=0 hist = {st: 0.0 for st in self.hmm.states} for (st,wt) in zip(self.particles,self.weights): hist[st]+=wt tot += wt return {st:hist[st]/tot for st in hist} def resample_particles(self): """resamples to give a new set of particles.""" self.particles = resample(self.particles, self.weights, len(self.particles)) self.weights = [1] * len(self.particles) hmm1pf1 = HMM_particle_filter(hmm1) # HMM_particle_filter.max_display_level = 2 # show each step # hmm1pf1.filter([{'m1':0, 'm2':1, 'm3':1}, {'m1':1, 'm2':0, 'm3':1}]) # hmm1pf2 = HMM_particle_filter(hmm1) # hmm1pf2.filter([{'m1':1, 'm2':0, 'm3':0}, {'m1':0, 'm2':1, 'm3':0}, {'m1':1, 'm2':0, 'm3':0}, # {'m1':0, 'm2':0, 'm3':0}, {'m1':0, 'm2':0, 'm3':0}, {'m1':0, 'm2':0, 'm3':0}, # {'m1':0, 'm2':0, 'm3':0}, {'m1':0, 'm2':0, 'm3':1}, {'m1':0, 'm2':0, 'm3':1}, # {'m1':0, 'm2':0, 'm3':1}]) # hmm1pf3 = HMM_particle_filter(hmm1) # hmm1pf3.filter([{'m1':1, 'm2':0, 'm3':0}, {'m1':0, 'm2':0, 'm3':0}, {'m1':1, 'm2':0, 'm3':0}, {'m1':1, 'm2':0, 'm3':1}]) def simulate(hmm,horizon): """returns a pair of (state sequence, observation sequence) of length horizon. for each time t, the agent is in state_sequence[t] and observes observation_sequence[t] """ state = sample_one(hmm.indist) obsseq=[] stateseq=[] for time in range(horizon): stateseq.append(state) newobs = {obs:sample_one({0:1-hmm.pobs[obs][state],1:hmm.pobs[obs][state]}) for obs in hmm.obsvars} obsseq.append(newobs) state = sample_one(hmm.trans[state]) return stateseq,obsseq def simobs(hmm,stateseq): """returns observation sequence for the state sequence""" obsseq=[] for state in stateseq: newobs = {obs:sample_one({0:1-hmm.pobs[obs][state],1:hmm.pobs[obs][state]}) for obs in hmm.obsvars} obsseq.append(newobs) return obsseq def create_eg(hmm,n): """Create an annotated example for horizon n""" seq,obs = simulate(hmm,n) print("True state sequence:",seq) print("Sequence of observations:\n",obs) hmmfilter = HMM_VE_filter(hmm) dist = hmmfilter.filter(obs) print("Resulting distribution over states:\n",dist)