Ticket #8547: trac-8547-take2.patch
File trac-8547-take2.patch, 243.4 KB (added by , 12 years ago) |
---|
-
module_list.py
# HG changeset patch # User William Stein <wstein@gmail.com> # Date 1269812759 25200 # Node ID e964f241c1d92393c7c5bc9b24687e72cdb9a218 # Parent a0a95c2cf60608c1bdedef3475e8a515ba62fead trac 8547 -- take 2 -- implement hidden markov models in Cython from scratch diff --git a/module_list.py b/module_list.py
a b 1334 1334 ## 1335 1335 ################################ 1336 1336 1337 Extension('sage.stats.hmm.chmm', 1338 sources = ['sage/stats/hmm/chmm.pyx'], 1339 libraries = ['ghmm'], 1340 include_dirs=numpy_include_dirs), 1337 Extension('sage.stats.hmm.util', 1338 sources = ['sage/stats/hmm/util.pyx']), 1339 1340 Extension('sage.stats.hmm.distributions', 1341 sources = ['sage/stats/hmm/distributions.pyx']), 1341 1342 1342 1343 Extension('sage.stats.hmm.hmm', 1343 sources = ['sage/stats/hmm/hmm.pyx'], 1344 libraries = ['ghmm'], 1345 include_dirs=numpy_include_dirs), 1344 sources = ['sage/stats/hmm/hmm.pyx']), 1345 1346 Extension('sage.stats.hmm.chmm', 1347 sources = ['sage/stats/hmm/chmm.pyx']), 1348 1349 Extension('sage.stats.intlist', 1350 sources = ['sage/stats/intlist.pyx']), 1346 1351 1347 1352 ################################ 1348 1353 ## -
sage/finance/time_series.pxd
diff --git a/sage/finance/time_series.pxd b/sage/finance/time_series.pxd
a b 1 1 cdef class TimeSeries: 2 2 cdef double* _values 3 3 cdef Py_ssize_t _length 4 cpdef rescale(self, double s) 5 cpdef double sum(self) -
sage/finance/time_series.pyx
diff --git a/sage/finance/time_series.pyx b/sage/finance/time_series.pyx
a b 68 68 """ 69 69 self._values = NULL 70 70 71 def __init__(self, values ):71 def __init__(self, values, bint initialize=True): 72 72 """ 73 73 Initialize new time series. 74 74 75 75 INPUT: 76 values -- integer (number of values) or an iterable of floats 76 77 - values -- integer (number of values) or an iterable of 78 floats 79 80 - initialize -- bool (default: True); if False, do not 81 bother to zero out the entries of the new time series. 82 For large series that you are going to just fill in, 83 this can be way faster. 77 84 78 85 EXAMPLES: 79 This implicitly calls init. 86 87 This implicitly calls init.:: 88 80 89 sage: finance.TimeSeries([pi, 3, 18.2]) 81 90 [3.1416, 3.0000, 18.2000] 82 91 83 Conversion from a numpy 1d array, which is very fast. 92 Conversion from a numpy 1d array, which is very fast.:: 93 84 94 sage: v = finance.TimeSeries([1..5]) 85 95 sage: w = v.numpy() 86 96 sage: finance.TimeSeries(w) 87 97 [1.0000, 2.0000, 3.0000, 4.0000, 5.0000] 88 98 89 Conversion from an n-dimensional numpy array also works: 99 Conversion from an n-dimensional numpy array also works:: 100 90 101 sage: import numpy 91 102 sage: v = numpy.array([[1,2], [3,4]], dtype=float); v 92 103 array([[ 1., 2.], … … 98 109 sage: u = numpy.array([[1,2],[3,4]]) 99 110 sage: finance.TimeSeries(u) 100 111 [1.0000, 2.0000, 3.0000, 4.0000] 112 113 For speed purposes we don't initialize (so value is garbage):: 114 115 sage: t = finance.TimeSeries(10, initialize=False) 101 116 """ 102 117 cdef Vector_real_double_dense z 103 118 cdef cnumpy.ndarray np 104 119 cdef double *np_data 105 120 cdef unsigned int j 106 121 if isinstance(values, (int, long, Integer)): 107 values = [0.0]*values 122 self._length = values 123 values = None 108 124 elif PY_TYPE_CHECK(values, Vector_real_double_dense) or PY_TYPE_CHECK(values, cnumpy.ndarray): 109 125 if PY_TYPE_CHECK(values, Vector_real_double_dense): 110 126 np = values._vector_numpy … … 130 146 return 131 147 else: 132 148 values = [float(x) for x in values] 133 self._length = len(values) 149 self._length = len(values) 150 134 151 self._values = <double*> sage_malloc(sizeof(double) * self._length) 135 152 if self._values == NULL: 136 153 raise MemoryError 154 if not initialize: return 137 155 cdef Py_ssize_t i 138 for i from 0 <= i < self._length: 139 self._values[i] = values[i] 156 if values is not None: 157 for i from 0 <= i < self._length: 158 self._values[i] = values[i] 159 else: 160 for i from 0 <= i < self._length: 161 self._values[i] = 0 140 162 141 163 def __reduce__(self): 142 164 """ … … 265 287 """ 266 288 if prec is None: prec = digits 267 289 format = '%.' + str(prec) + 'f' 268 v = self.list() 269 if len(v) > max_print: 270 v0 = v[:max_print//2] 271 v1 = v[-max_print//2:] 290 if len(self) > max_print: 291 v0 = self[:max_print//2] 292 v1 = self[-max_print//2:] 272 293 return '[' + ', '.join([format%x for x in v0]) + ' ... ' + \ 273 294 ', '.join([format%x for x in v1]) + ']' 274 295 else: 275 return '[' + ', '.join([format%x for x in v]) + ']'296 return '[' + ', '.join([format%x for x in self]) + ']' 276 297 277 298 def __len__(self): 278 299 """ … … 363 384 memcpy(t._values, self._values + start, sizeof(double)*t._length) 364 385 return t 365 386 else: 366 if i < 0: 367 i += self._length 368 if i < 0: 387 j = i 388 if j < 0: 389 j += self._length 390 if j < 0: 369 391 raise IndexError, "TimeSeries index out of range" 370 elif i>= self._length:392 elif j >= self._length: 371 393 raise IndexError, "TimeSeries index out of range" 372 return self._values[ i]394 return self._values[j] 373 395 374 396 def __setitem__(self, Py_ssize_t i, double x): 375 397 """ … … 643 665 sage: v.list() 644 666 [1.0, -4.0, 3.0, -2.5, -4.0, 3.0] 645 667 """ 646 v = [0.0]*self._length647 668 cdef Py_ssize_t i 648 for i from 0 <= i < self._length: 649 v[i] = self._values[i] 650 return v 669 return [self._values[i] for i in range(self._length)] 651 670 652 671 def log(self): 653 672 """ … … 666 685 [2.7183, 0.0183, 20.0855, 0.0821, 0.0183, 20.0855] 667 686 sage: v.exp().log() 668 687 [1.0000, -4.0000, 3.0000, -2.5000, -4.0000, 3.0000] 688 689 Log of 0 gives -inf: 690 691 sage: finance.TimeSeries([1,0,3]).log() 692 [0.0000, -inf, 1.0986] 669 693 """ 670 694 cdef Py_ssize_t i 671 for i from 0 <= i < self._length:672 if self._values[i] <= 0:673 raise ValueError, "every entry of self must be positive."674 675 695 cdef TimeSeries t = new_time_series(self._length) 676 696 for i from 0 <= i < self._length: 677 697 t._values[i] = log(self._values[i]) … … 801 821 t._values[i] = self._values[i*k] 802 822 return t 803 823 824 cpdef rescale(self, double s): 825 """ 826 Change self by multiplying every value in the series by s. 827 828 INPUT: 829 s -- float 830 831 EXAMPLES: 832 sage: v = finance.TimeSeries([5,4,1.3,2,8,10,3,-5]); v 833 [5.0000, 4.0000, 1.3000, 2.0000, 8.0000, 10.0000, 3.0000, -5.0000] 834 sage: v.rescale(0.5) 835 sage: v 836 [2.5000, 2.0000, 0.6500, 1.0000, 4.0000, 5.0000, 1.5000, -2.5000] 837 """ 838 for i from 0 <= i < self._length: 839 self._values[i] = self._values[i] * s 840 804 841 def scale(self, double s): 805 842 """ 806 843 Return new time series obtained by multiplying every value in the series by s. … … 1067 1104 t._values[i] = s 1068 1105 return t 1069 1106 1070 defsum(self):1107 cpdef double sum(self): 1071 1108 """ 1072 1109 Return the sum of all the entries of self. If self has 1073 1110 length 0, returns 0. … … 1555 1592 """ 1556 1593 if self._length == 0: 1557 1594 raise ValueError, "min() arg is an empty sequence" 1558 return 0.0, -11559 1595 cdef Py_ssize_t i, j 1560 1596 cdef double s = self._values[0] 1561 1597 j = 0 … … 1589 1625 """ 1590 1626 if self._length == 0: 1591 1627 raise ValueError, "max() arg is an empty sequence" 1592 cdef Py_ssize_t i, j 1628 cdef Py_ssize_t i, j = 0 1593 1629 cdef double s = self._values[0] 1594 1630 for i from 1 <= i < self._length: 1595 1631 if self._values[i] > s: … … 1743 1779 bins -- positive integer (default: 50) 1744 1780 normalize -- (default: True) whether to normalize so the total 1745 1781 area in the bars of the histogram is 1. 1746 **kwds -- passed to the bar_chart function1747 1782 OUTPUT: 1748 1783 a histogram plot 1749 1784 … … 1751 1786 sage: v = finance.TimeSeries([1..50]) 1752 1787 sage: v.plot_histogram(bins=10) 1753 1788 """ 1754 from sage.plot.all import bar_chart,polygon1789 from sage.plot.all import polygon 1755 1790 counts, intervals = self.histogram(bins, normalize=normalize) 1756 1791 s = 0 1757 1792 for i, (x0,x1) in enumerate(intervals): … … 2239 2274 2240 2275 def unpickle_time_series_v1(v, Py_ssize_t n): 2241 2276 """ 2242 Version 0unpickle method.2277 Version 1 unpickle method. 2243 2278 2244 2279 INPUT: 2245 2280 v -- a raw char buffer -
sage/matrix/matrix_complex_double_dense.pxd
diff --git a/sage/matrix/matrix_complex_double_dense.pxd b/sage/matrix/matrix_complex_double_dense.pxd
a b 1 import matrix_double_dense 2 1 3 cimport matrix_double_dense 2 4 3 5 cdef class Matrix_complex_double_dense(matrix_double_dense.Matrix_double_dense): -
sage/matrix/matrix_complex_double_dense.pyx
diff --git a/sage/matrix/matrix_complex_double_dense.pyx b/sage/matrix/matrix_complex_double_dense.pyx
a b 36 36 # The full text of the GPL is available at: 37 37 # http://www.gnu.org/licenses/ 38 38 ############################################################################## 39 40 import matrix_double_dense 41 39 42 from sage.rings.complex_double import CDF 40 43 41 44 cimport numpy as cnumpy -
sage/matrix/matrix_double_dense.pyx
diff --git a/sage/matrix/matrix_double_dense.pyx b/sage/matrix/matrix_double_dense.pyx
a b 52 52 cnumpy.import_array() 53 53 54 54 55 from matrix_complex_double_dense import Matrix_complex_double_dense56 57 55 cdef class Matrix_double_dense(matrix_dense.Matrix_dense): 58 56 """ 59 57 Base class for matrices over the Real Double Field and the Complex -
sage/stats/all.py
diff --git a/sage/stats/all.py b/sage/stats/all.py
a b 3 3 4 4 import hmm.all as hmm 5 5 6 from sage.finance.time_series import TimeSeries 7 from intlist import IntList 8 -
sage/stats/hmm/__init__.py
diff --git a/sage/stats/hmm/__init__.py b/sage/stats/hmm/__init__.py
a b 1 # Hidden Markov Models1 # hmm -
sage/stats/hmm/all.py
diff --git a/sage/stats/hmm/all.py b/sage/stats/hmm/all.py
a b 1 1 ############################################################################# 2 # Copyright (C) 2008 William Stein <wstein@gmail.com> 3 # Distributed under the terms of the GNU General Public License (GPL), 4 # version 2 or any later version at your option. 2 # Copyright (C) 2010 William Stein <wstein@gmail.com> 3 # Distributed under the terms of the GNU General Public License (GPL), v2+. 5 4 # The full text of the GPL is available at: 6 5 # http://www.gnu.org/licenses/ 7 6 ############################################################################# 8 7 9 8 from hmm import DiscreteHiddenMarkovModel 10 from chmm import GaussianHiddenMarkovModel, GaussianMixtureHiddenMarkovModel 9 10 from chmm import (GaussianHiddenMarkovModel, 11 GaussianMixtureHiddenMarkovModel) 12 13 from distributions import GaussianMixtureDistribution 14 -
deleted file sage/stats/hmm/chmm.pxd
diff --git a/sage/stats/hmm/chmm.pxd b/sage/stats/hmm/chmm.pxd deleted file mode 100644
+ - 1 #############################################################################2 # Copyright (C) 2008 William Stein <wstein@gmail.com>3 # Distributed under the terms of the GNU General Public License (GPL),4 # version 2 or any later version at your option.5 # The full text of the GPL is available at:6 # http://www.gnu.org/licenses/7 #############################################################################8 9 from hmm cimport HiddenMarkovModel10 11 cdef extern from "ghmm/ghmm.h":12 cdef int GHMM_kContinuousHMM13 14 cdef extern from "ghmm/smodel.h":15 #############################################################16 # Continuous emission density types17 #############################################################18 cdef enum ghmm_density_t:19 normal,20 normal_right, # right tail21 normal_approx,22 normal_left, # left tail23 uniform,24 binormal,25 multinormal,26 density_number27 28 #############################################################29 # Continuous emission30 #############################################################31 # Mean and variance types32 cdef union mean_t:33 double val34 double *vec35 36 cdef union variance_t:37 double val38 double *mat39 40 cdef struct ghmm_c_emission:41 # Flag for density function for each component of the mixture42 # 0: normal density43 # 1: truncated normal (right side) density44 # 2: approximated normal density45 # 3: truncated normal (left side)46 # 4: uniform distribution47 # 6: multivariate normal48 int type49 # dimension > 1 for multivariate normals50 int dimension51 # mean for output functions (pointer to mean vector for multivariate)52 mean_t mean53 variance_t variance54 # pointer to inverse of covariance matrix if multivariate normal;55 # else NULL56 double *sigmainv57 # determinant of covariance matrix for multivariate normal58 double det59 # Cholesky decomposition of covariance matrix A,60 # if A = G*G' sigmacd only holds G61 double *sigmacd62 # minimum of uniform distribution or left boundary for63 # right-tail gaussians64 double min65 # maximum of uniform distribution or right boundary for66 # left-tail gaussians67 double max68 # if fixed != 0 the parameters of the density are fixed69 int fixed70 71 #############################################################72 # Continous Emission States73 #############################################################74 ctypedef struct ghmm_cstate:75 # Number of output densities per state76 int M77 # initial prob.78 double pi79 # IDs of successor states80 int *out_id81 # IDs of predecessor states82 int *in_id83 # transition probs to successor states. It is a84 # matrix in case of mult. transition matrices (cos > 1)85 double **out_a86 # transition probs from predecessor states. It is a87 # matrix in case of mult. transition matrices (cos > 1)88 double **in_a89 # number of successor states90 int out_states91 # number of predecessor states92 int in_states93 # weight vector for output function components94 double *c95 # flag for fixation of parameter. If fix = 1 do not change parameters of96 # output functions, and if fix = 0 do normal training. Default is 0.97 int fix98 # vector of ghmm_c_emission99 # (type and parameters of output function components)100 ghmm_c_emission *e101 # contains a description of the state (null terminated utf-8)102 char *desc103 # x coordinate position for graph representation plotting104 int xPosition105 # y coordinate position for graph representation plotting106 int yPosition107 108 109 double **ighmm_cmatrix_alloc (int nrows, int ncols)110 111 #############################################################112 # Continous Hidden Markov Model113 #############################################################114 ctypedef struct ghmm_cmodel115 ctypedef struct ghmm_cmodel_class_change_context116 117 ctypedef struct ghmm_cmodel:118 # Number of states119 int N120 # Maximun number of components in the states121 int M122 # Number of dimensions of the emission components.123 # NOTE: All emissions must have the same number of dimensions.124 int dim125 # The ghmm_cmodel includes two continuous models126 # cos=1: one transition matrix127 # cos>1: (integer) extension for models with several matrices.128 int cos129 # Prior for a priori probability of the model.130 # -1 means no prior specified (all models have equal prob. a priori.)131 double prior132 # Contains a arbitrary name for the model (null terminated utf-8)133 char * name134 # Contains bit flags for various model extensions such as135 # kSilentStates (see ghmm.h for a complete list).136 int model_type137 # All states of the model. Transition probabilities are part of the states.138 ghmm_cstate *s139 # pointer to a ghmm_cmodel_class_change_context struct140 # necessary for multiple transition classes141 ghmm_cmodel_class_change_context *class_change142 143 int ghmm_cmodel_class_change_alloc(ghmm_cmodel * smo)144 145 #############################################################146 # Continous Sequence of states147 #148 # Sequence structure for double sequences. Contains an array of149 # sequences and corresponding data like sequence labels, sequence150 # weights, etc. Sequences may have different length.151 # multi-dimension sequences are linearized152 #############################################################153 ctypedef struct ghmm_cseq:154 # sequence array. sequence[i][j] = j-th symbol of i-th seq.155 # sequence[i][D * j] = first dimension of j-th observation of i-th sequence156 double **seq157 # array of sequence lengths158 int *seq_len159 # array of sequence IDs160 double *seq_id161 # positive! sequence weights. default is 1 = no weight162 double *seq_w163 # total number of sequences164 long seq_number165 # reserved space for sequences is always >= seq_number166 long capacity167 # sum of sequence weights168 double total_w169 # total number of dimensions170 int dim171 # flags (internal)172 unsigned int flags173 # obsolete174 int* seq_label175 176 int ghmm_cseq_free (ghmm_cseq ** sq)177 178 #############################################################179 # Continous Baum-Welch algorithm context180 #############################################################181 ctypedef struct ghmm_cmodel_baum_welch_context:182 # pointer of continuous model183 ghmm_cmodel *smo184 # sequence pointer185 ghmm_cseq *sqd186 # calculated log likelihood187 double *logp188 # leave reestimation loop if diff. between successive logp values189 # is smaller than eps190 double eps191 # max. no of iterations192 int max_iter193 194 195 int ghmm_cmodel_free (ghmm_cmodel ** smo)196 197 # Normalizes initial and transition probabilities and mixture weights.198 # 0 on success / -1 on error199 int ghmm_cmodel_normalize(ghmm_cmodel *smo)200 201 202 # Produces sequences computed using a given model. All memory that203 # is needed for the sequences is allocated inside the function. It204 # is possible to define the length of the sequences global205 # (global_len > 0) or it can be set inside the function, when a206 # final state in the model is reached (a state with no output). If207 # the model has no final state, the sequences will have length208 # MAX_SEQ_LEN.209 #210 # INPUT:211 # smo -- model212 # seed -- initial parameter for the random number generator;213 # if seed == 0, don't reinitialize random number generator214 # global_en -- length of sequence. If 0, sequence ended automatically215 # when the final state is reached216 # seq_number-- number of sequences217 # Tmax -- maximal sequence length (or -1)218 #219 # OUTPUT:220 # pointer to an array of sequences221 #222 223 ghmm_cseq *ghmm_cmodel_generate_sequences (ghmm_cmodel * smo, int seed,224 int global_len, long seq_number,225 int Tmax)226 227 # Computes sum over all sequence of seq_w * log( P ( O|lambda )).228 # If a sequence can't be generated by smo error cost of seq_w *229 # PENALTY_LOGP are imposed.230 # OUTPUT: number of evaluated sequences; -1: error231 int ghmm_cmodel_likelihood (ghmm_cmodel * smo, ghmm_cseq * sqd, double *log_p)232 233 234 # Viterbi algorithm: calculation of the viterbi path (best possible235 # state sequence for a given sequence for a given model (smo)). Also236 # calculates log(p) according to this path, the matrices in the local_store237 # struct are allocated using stat_matrix_d_alloc.238 int *ghmm_cmodel_viterbi (ghmm_cmodel * smo, double *o, int n, double *log_p)239 240 cdef extern from "ghmm/sreestimate.h":241 ctypedef struct ghmm_cmodel_class_change_context:242 # Names of class change module/function (for python callback)243 char *python_module244 char *python_function245 # index of current sequence246 int k247 # pointer to class function248 int (*get_class) (ghmm_cmodel*, double*, int, int)249 # space for any data necessary for class switch, USER is RESPONSIBLE250 void *user_data251 252 # Baum-Welch Algorithm for continuous HMMs253 # Training of model parameter with multiple double sequences (incl. scaling).254 # New parameters set directly in hmm (no storage of previous values!). Matrices255 # are allocated with stat_matrix_d_alloc.256 int ghmm_cmodel_baum_welch (ghmm_cmodel_baum_welch_context * cs)257 258 259 cdef class ContinuousHiddenMarkovModel(HiddenMarkovModel):260 cdef ghmm_cmodel* m261 cdef bint initialized262 263 cdef class GaussianHiddenMarkovModel(ContinuousHiddenMarkovModel):264 pass265 266 -
sage/stats/hmm/chmm.pyx
diff --git a/sage/stats/hmm/chmm.pyx b/sage/stats/hmm/chmm.pyx
a b 1 """ nodoctest2 Continuous Hidden Markov Models1 """ 2 Continuous Emission Hidden Markov Models 3 3 4 AUTHOR S:5 -- William Stein (2008) 6 -- The authors of GHMM http://ghmm.sourceforge.net/4 AUTHOR: 5 6 - William Stein, 2010-03 7 7 """ 8 8 9 9 ############################################################################# 10 # Copyright (C) 2008 William Stein <wstein@gmail.com> 11 # Distributed under the terms of the GNU General Public License (GPL), 12 # version 2 or any later version at your option. 10 # Copyright (C) 2010 William Stein <wstein@gmail.com> 11 # Distributed under the terms of the GNU General Public License (GPL) 13 12 # The full text of the GPL is available at: 14 13 # http://www.gnu.org/licenses/ 15 14 ############################################################################# 16 15 17 16 include "../../ext/stdsage.pxi" 17 include "../../ext/cdefs.pxi" 18 include "../../ext/interrupt.pxi" 18 19 19 include "misc.pxi" 20 cdef extern from "math.h": 21 double log(double) 22 double sqrt(double) 23 double exp(double) 24 int isnormal(double) 25 int isfinite(double) 20 26 21 from sage.misc.randstate import random 22 from sage.misc.flatten import flatten 27 import math 28 29 from sage.misc.flatten import flatten 30 from sage.matrix.all import is_Matrix 31 32 cdef double sqrt2pi = sqrt(2*math.pi) 33 cdef double NaN = 1.0/0.0 23 34 24 35 from sage.finance.time_series cimport TimeSeries 36 from sage.stats.intlist cimport IntList 25 37 26 cdef extern from "math.h": 27 double sqrt(double) 28 29 cdef class ContinuousHiddenMarkovModel(HiddenMarkovModel): 38 from hmm cimport HiddenMarkovModel 39 from util cimport HMM_Util 40 from distributions cimport GaussianMixtureDistribution 41 42 cdef HMM_Util util = HMM_Util() 43 44 from sage.misc.randstate cimport current_randstate, randstate 45 46 47 # DELETE THIS FUNCTION WHEN MOVE Gaussian stuff to distributions.pyx!!! 48 cdef double random_normal(double mean, double std, randstate rstate): 49 # Ported from http://users.tkk.fi/~nbeijar/soft/terrain/source_o2/boxmuller.c 50 # This the box muller algorithm. 51 # Client code can get the current random state from: 52 # cdef randstate rstate = current_randstate() 53 54 cdef double x1, x2, w, y1, y2 55 while True: 56 x1 = 2*rstate.c_rand_double() - 1 57 x2 = 2*rstate.c_rand_double() - 1 58 w = x1*x1 + x2*x2 59 if w < 1: break 60 w = sqrt( (-2*log(w))/w ) 61 y1 = x1 * w 62 return mean + y1*std 63 64 cdef class GaussianHiddenMarkovModel(HiddenMarkovModel): 30 65 """ 31 Abstract base class for continuous hidden Markov models.66 GaussianHiddenMarkovModel(A, B, pi) 32 67 33 68 Gaussian emissions Hidden Markov Model. 69 34 70 INPUT: 35 A -- square matrix (or list of lists) 36 B -- list or matrix with numbers of rows equal to number of rows of A; 37 each row defines the emissions 38 pi -- list 39 name -- (default: None); a string 71 - A -- matrix; the N x N transition matrix 72 - B -- list of pairs (mu,sigma) that define the distributions 73 - pi -- initial state probabilities 40 74 41 EXAMPLES: 42 sage: sage.stats.hmm.chmm.ContinuousHiddenMarkovModel([[1.0]], [(-0.0,10.0)], [1], "model") 43 <sage.stats.hmm.chmm.ContinuousHiddenMarkovModel object at ...> 75 EXAMPLES:: 76 77 We illustrate the primary functions with an example 2-state Gaussian HMM:: 78 79 sage: m = hmm.GaussianHiddenMarkovModel([[.1,.9],[.5,.5]], [(1,1), (-1,1)], [.5,.5]); m 80 Gaussian Hidden Markov Model with 2 States 81 Transition matrix: 82 [0.1 0.9] 83 [0.5 0.5] 84 Emission parameters: 85 [(1.0, 1.0), (-1.0, 1.0)] 86 Initial probabilities: [0.5000, 0.5000] 87 88 We query the defining transition matrix, emission parameters, and 89 initial state probabilities:: 90 91 sage: m.transition_matrix() 92 [0.1 0.9] 93 [0.5 0.5] 94 sage: m.emission_parameters() 95 [(1.0, 1.0), (-1.0, 1.0)] 96 sage: m.initial_probabilities() 97 [0.5000, 0.5000] 98 99 We obtain a sample sequence with 10 entries in it, and compute the 100 logarithm of the probability of obtaining his sequence, given the 101 model:: 102 103 sage: obs = m.sample(10); obs 104 [-1.6835, 0.0635, -2.1688, 0.3043, -0.3188, -0.7835, 1.0398, -1.3558, 1.0882, 0.4050] 105 sage: m.log_likelihood(obs) 106 -15.2262338077988... 107 108 We compute the Viterbi path, and probability that the given path 109 of states produced obs:: 110 111 sage: m.viterbi(obs) 112 ([1, 0, 1, 0, 1, 1, 0, 1, 0, 1], -16.677382701707881) 113 114 We use the Baum-Welch iterative algorithm to find another model 115 for which our observation sequence is more likely:: 116 117 sage: m.baum_welch(obs) 118 (-10.610333495739708, 14) 119 sage: m.log_likelihood(obs) 120 -10.610333495739708 121 122 Notice that running Baum-Welch changed our model:: 123 124 sage: m 125 Gaussian Hidden Markov Model with 2 States 126 Transition matrix: 127 [ 0.415498136619 0.584501863381] 128 [ 0.999999317425 6.82574625899e-07] 129 Emission parameters: 130 [(0.417888242712, 0.517310966436), (-1.50252086313, 0.508551283606)] 131 Initial probabilities: [0.0000, 1.0000] 44 132 """ 45 def __init__(self, A, B, pi, name): 133 cdef TimeSeries B, prob 134 cdef int n_out 135 136 def __init__(self, A, B, pi, bint normalize=True): 46 137 """ 47 Constructor for continuous Hidden Markov model abstract base class. 138 Create a Gaussian emissions HMM with transition probability 139 matrix A, normal emissions given by B, and initial state 140 probability distribution pi. 48 141 49 EXAMPLES: 50 This class is an abstract base class, so shouldn't ever be constructed by users. 51 sage: sage.stats.hmm.chmm.ContinuousHiddenMarkovModel([[1.0]], [(0.0,1.0)], [1], None) 52 <sage.stats.hmm.chmm.ContinuousHiddenMarkovModel object at ...> 142 INPUT:: 143 144 - A -- a list of lists or a square N x N matrix, whose 145 (i,j) entry gives the probability of transitioning from 146 state i to state j. 147 148 - B -- a list of N pairs (mu,std), where if B[i]=(mu,std), 149 then the probability distribution associated with state i 150 normal with mean mu and standard deviation std. 151 152 - pi -- the probabilities of starting in each initial 153 state, i.e,. pi[i] is the probability of starting in 154 state i. 155 156 - normalize --bool (default: True); if given, input is 157 normalized to define valid probability distributions, 158 e.g., the entries of A are made nonnegative and the rows 159 sum to 1. 160 161 EXAMPLES:: 162 163 164 sage: hmm.GaussianHiddenMarkovModel([[.1,.9],[.5,.5]], [(1,1), (-1,1)], [.5,.5]) 165 Gaussian Hidden Markov Model with 2 States 166 Transition matrix: 167 [0.1 0.9] 168 [0.5 0.5] 169 Emission parameters: 170 [(1.0, 1.0), (-1.0, 1.0)] 171 Initial probabilities: [0.5000, 0.5000] 172 173 We input a model in which both A and pi have to be 174 renormalized to define valid probability distributions:: 175 176 sage: hmm.GaussianHiddenMarkovModel([[-1,.7],[.3,.4]], [(1,1), (-1,1)], [-1,.3]) 177 Gaussian Hidden Markov Model with 2 States 178 Transition matrix: 179 [ 0.0 1.0] 180 [0.428571428571 0.571428571429] 181 Emission parameters: 182 [(1.0, 1.0), (-1.0, 1.0)] 183 Initial probabilities: [0.0000, 1.0000] 184 185 Bad things can happen:: 186 187 sage: hmm.GaussianHiddenMarkovModel([[-1,.7],[.3,.4]], [(1,1), (-1,1)], [-1,.3], normalize=False) 188 Gaussian Hidden Markov Model with 2 States 189 Transition matrix: 190 [-1.0 0.7] 191 [ 0.3 0.4] 192 ... 53 193 """ 54 self.initialized = False 55 HiddenMarkovModel.__init__(self, A, B, pi) 56 self.m = <ghmm_cmodel*> safe_malloc(sizeof(ghmm_cmodel)) 57 # Set number of states 58 self.m.N = self.A.nrows() 59 # Assign model identifier (the name) if specified 60 if name is not None: 61 name = str(name) 62 self.m.name = <char*> safe_malloc(len(name)+1) 63 strcpy(self.m.name, name) 194 self.pi = util.initial_probs_to_TimeSeries(pi, normalize) 195 self.N = len(self.pi) 196 self.A = util.state_matrix_to_TimeSeries(A, self.N, normalize) 197 198 # B should be a matrix of N rows, with column 0 the mean and 1 199 # the standard deviation. 200 if is_Matrix(B): 201 B = B.list() 64 202 else: 65 self.m.name = NULL 66 203 B = flatten(B) 204 self.B = TimeSeries(B) 205 self.probability_init() 67 206 68 def name(self):207 def __cmp__(self, other): 69 208 """ 70 Return the name of this model. 209 Compare self and other, which must both be GaussianHiddenMarkovModel's. 210 211 EXAMPLES:: 212 213 sage: m = hmm.GaussianHiddenMarkovModel([[1]], [(0,1)], [1]) 214 sage: n = hmm.GaussianHiddenMarkovModel([[1]], [(1,1)], [1]) 215 sage: m < n 216 True 217 sage: m == m 218 True 219 sage: n > m 220 True 221 sage: n < m 222 False 223 """ 224 if not isinstance(other, GaussianHiddenMarkovModel): 225 raise ValueError 226 return cmp(self.__reduce__()[1], other.__reduce__()[1]) 227 228 def __getitem__(self, Py_ssize_t i): 229 """ 230 Return the mean and standard distribution for the i-th state. 231 232 INPUT: 233 234 - i -- integer 71 235 72 236 OUTPUT: 73 string or None74 237 75 EXAMPLES: 76 sage: m = hmm.GaussianHiddenMarkovModel([[0.4,0.6],[0.1,0.9]], [(0.0,1.0),(1,1)], [1,0], 'Test Model') 77 sage: m.name() 78 'Test Model' 238 - 2 floats 79 239 80 If the model is not explicitly named then this function returns None: 81 sage: m = hmm.GaussianHiddenMarkovModel([[0.4,0.6],[0.1,0.9]], [(0.0,1.0),(1,1)], [1,0]) 82 sage: m.name() 83 sage: m.name() is None 240 EXAMPLES:: 241 242 sage: m = hmm.GaussianHiddenMarkovModel([[.1,.9],[.5,.5]], [(1,.5), (-2,.3)], [.5,.5]) 243 sage: m[0] 244 (1.0, 0.5) 245 sage: m[1] 246 (-2.0, 0.29999999999999999) 247 sage: m[-1] 248 (-2.0, 0.29999999999999999) 249 sage: m[3] 250 Traceback (most recent call last): 251 ... 252 IndexError: index out of range 253 sage: m[-3] 254 Traceback (most recent call last): 255 ... 256 IndexError: index out of range 257 """ 258 if i < 0: 259 i += self.N 260 if i < 0 or i >= self.N: 261 raise IndexError, 'index out of range' 262 263 # TODO: change to be a normal distribution class 264 return self.B[2*i], self.B[2*i+1] 265 266 def __reduce__(self): 267 """ 268 Used in pickling. 269 270 EXAMPLES:: 271 272 sage: G = hmm.GaussianHiddenMarkovModel([[1]], [(0,1)], [1]) 273 sage: loads(dumps(G)) == G 84 274 True 85 275 """ 86 if self.m.name: 87 s = str(self.m.name) 88 return s 276 return unpickle_gaussian_hmm_v1, \ 277 (self.A, self.B, self.pi, self.prob, self.n_out) 278 279 def emission_parameters(self): 280 """ 281 Return the parameters that define the normal distributions 282 associated to all of the states. 283 284 OUTPUT: 285 286 - a list B of pairs B[i] = (mu, std), such that the 287 distribution associated to state i is normal with mean 288 mu and standard deviation std. 289 290 EXAMPLES:: 291 292 sage: hmm.GaussianHiddenMarkovModel([[.1,.9],[.5,.5]], [(1,.5), (-1,3)], [.1,.9]).emission_parameters() 293 [(1.0, 0.5), (-1.0, 3.0)] 294 """ 295 cdef Py_ssize_t i 296 from sage.rings.all import RDF 297 return [(RDF(self.B[2*i]),RDF(self.B[2*i+1])) for i in range(self.N)] 298 299 def __repr__(self): 300 """ 301 Return string representation. 302 303 EXAMPLES:: 304 305 sage: hmm.GaussianHiddenMarkovModel([[.1,.9],[.5,.5]], [(1,.5), (-1,3)], [.1,.9]).__repr__() 306 'Gaussian Hidden Markov Model with 2 States\nTransition matrix:\n[0.1 0.9]\n[0.5 0.5]\nEmission parameters:\n[(1.0, 0.5), (-1.0, 3.0)]\nInitial probabilities: [0.1000, 0.9000]' 307 """ 308 s = "Gaussian Hidden Markov Model with %s States"%self.N 309 s += '\nTransition matrix:\n%s'%self.transition_matrix() 310 s += '\nEmission parameters:\n%s'%self.emission_parameters() 311 s += '\nInitial probabilities: %s'%self.initial_probabilities() 312 return s 313 314 315 def generate_sequence(self, Py_ssize_t length): 316 """ 317 Return a sample of the given length from this HMM. 318 319 INPUT: 320 321 - length -- positive integer 322 323 OUTPUT: 324 325 - an IntList or list of emission symbols 326 - TimeSeries of emissions 327 328 EXAMPLES:: 329 330 sage: m =hmm.GaussianHiddenMarkovModel([[.1,.9],[.5,.5]], [(1,.5), (-1,3)], [.1,.9]) 331 sage: m.generate_sequence(5) 332 ([-3.0505, 0.5317, -4.5065, 0.6521, 1.0435], [1, 0, 1, 0, 1]) 333 sage: m.generate_sequence(0) 334 ([], []) 335 sage: m.generate_sequence(-1) 336 Traceback (most recent call last): 337 ... 338 ValueError: length must be nonnegative 339 """ 340 if length < 0: 341 raise ValueError, "length must be nonnegative" 342 343 # Create Integer lists for states and observations 344 cdef IntList states = IntList(length) 345 cdef TimeSeries obs = TimeSeries(length) 346 if length == 0: 347 return states, obs 348 349 # Setup variables, including random state. 350 cdef Py_ssize_t i, j 351 cdef randstate rstate = current_randstate() 352 cdef int q = 0 353 cdef double r, accum 354 355 # Choose the starting state. 356 # See the remark in hmm.pyx about how this should get 357 # replaced by some general fast discrete distribution code. 358 r = rstate.c_rand_double() 359 accum = 0 360 for i in range(self.N): 361 if r < self.pi._values[i] + accum: 362 q = i 363 else: 364 accum += self.pi._values[i] 365 366 states._values[0] = q 367 obs._values[0] = self.random_sample(q, rstate) 368 369 cdef double* row 370 cdef int O 371 _sig_on 372 for i in range(1, length): 373 accum = 0 374 row = self.A._values + q*self.N 375 r = rstate.c_rand_double() 376 for j in range(self.N): 377 if r < row[j] + accum: 378 q = j 379 break 380 else: 381 accum += row[j] 382 states._values[i] = q 383 obs._values[i] = self.random_sample(q, rstate) 384 _sig_off 385 386 return obs, states 387 388 cdef probability_init(self): 389 """ 390 Used internally to compute caching information that makes 391 certain computations in the Baum-Welch algorithm faster. 392 """ 393 self.prob = TimeSeries(2*self.N) 394 cdef int i 395 for i in range(self.N): 396 self.prob[2*i] = 1.0/(sqrt2pi*self.B[2*i+1]) 397 self.prob[2*i+1] = -1.0/(2*self.B[2*i+1]*self.B[2*i+1]) 398 399 cdef double random_sample(self, int state, randstate rstate): 400 """ 401 Return a random sample from the normal distribution associated 402 to the given state. 403 404 This is only used internally, and no bounds or other error 405 checking is done, so calling this improperly can lead to seg 406 faults. 407 408 INPUT: 409 410 - state -- integer 411 - rstate -- randstate instance 412 413 OUTPUT: 414 415 - double 416 """ 417 return random_normal(self.B._values[state*2], self.B._values[state*2+1], rstate) 418 419 cdef double probability_of(self, int state, double observation): 420 """ 421 Return the probability b_j(o) of seeing the given observation 422 given that we're in the given state j (=state). 423 424 This is a continuous probability, so this really returns a 425 number p such that the probability of a value in the interval 426 [o,o+d] is p*d. 427 428 INPUT: 429 430 - state -- integer 431 - observation -- double 432 433 OUTPUT: 434 435 - double 436 """ 437 # The code below is an optimized version of the following code: 438 # cdef double mean = self.B._values[2*state], \ 439 # std = self.B._values[2*state+1] 440 # return 1/(sqrt2pi*std) * \ 441 # exp(-(observation-mean)*(observation-mean)/(2*std*std)) 442 443 cdef double x = observation - self.B._values[2*state] # observation - mean 444 return self.prob._values[2*state] * exp(x*x*self.prob._values[2*state+1]) 445 446 def log_likelihood(self, obs): 447 """ 448 Return the logarithm of the probability that this model 449 produced the given observation sequence. 450 451 INPUT: 452 453 - obs -- sequence of observations 454 455 OUTPUT: 456 457 - float 458 459 EXAMPLES:: 460 461 sage: m = hmm.GaussianHiddenMarkovModel([[.1,.9],[.5,.5]], [(1,.5), (-1,3)], [.1,.9]) 462 sage: m.log_likelihood([1,1,1]) 463 -4.2978807660724856 464 sage: set_random_seed(0); s = m.sample(20) 465 sage: m.log_likelihood(s) 466 -40.115714129484... 467 """ 468 if len(obs) == 0: 469 return 1.0 470 if not isinstance(obs, TimeSeries): 471 obs = TimeSeries(obs) 472 return self._forward_scale(obs) 473 474 def _forward_scale(self, TimeSeries obs): 475 """ 476 Memory-efficient implementation of the forward algorithm (with scaling). 477 478 INPUT: 479 480 - obs -- an integer list of observation states. 481 482 OUTPUT: 483 484 - float -- the log of the probability that the model 485 produced this sequence 486 487 EXAMPLES:: 488 489 sage: m = hmm.GaussianHiddenMarkovModel([[.1,.9],[.5,.5]], [(1,.5), (-1,3)], [.1,.9]) 490 sage: m._forward_scale(stats.TimeSeries([1,-1,-1,1])) 491 -7.641988207069133 492 """ 493 cdef Py_ssize_t i, j, t, T = len(obs) 494 495 # The running some of the log probabilities 496 cdef double log_probability = 0, sum, a 497 498 cdef TimeSeries alpha = TimeSeries(self.N), \ 499 alpha2 = TimeSeries(self.N) 500 501 # Initialization 502 sum = 0 503 for i in range(self.N): 504 a = self.pi[i] * self.probability_of(i, obs._values[0]) 505 alpha[i] = a 506 sum += a 507 508 log_probability = log(sum) 509 alpha.rescale(1/sum) 510 511 # Induction 512 cdef double s 513 for t in range(1, T): 514 sum = 0 515 for j in range(self.N): 516 s = 0 517 for i in range(self.N): 518 s += alpha._values[i] * self.A._values[i*self.N + j] 519 a = s * self.probability_of(j, obs._values[t]) 520 alpha2._values[j] = a 521 sum += a 522 523 log_probability += log(sum) 524 for j in range(self.N): 525 alpha._values[j] = alpha2._values[j] / sum 526 527 # Termination 528 return log_probability 529 530 def viterbi(self, obs): 531 """ 532 Determine "the" hidden sequence of states that is most likely 533 to produce the given sequence seq of observations, along with 534 the probability that this hidden sequence actually produced 535 the observation. 536 537 INPUT: 538 539 - seq -- sequence of emitted ints or symbols 540 541 OUTPUT: 542 543 - list -- "the" most probable sequence of hidden states, i.e., 544 the Viterbi path. 545 546 - float -- log of probability that the observed sequence 547 was produced by the Viterbi sequence of states. 548 549 EXAMPLES:: 550 551 We find the optimal state sequence for a given model:: 552 553 sage: m = hmm.GaussianHiddenMarkovModel([[0.5,0.5],[0.5,0.5]], [(0,1),(10,1)], [0.5,0.5]) 554 sage: m.viterbi([0,1,10,10,1]) 555 ([0, 0, 1, 1, 0], -9.0604285688230899) 556 557 Another example in which the most likely states change based 558 on the last observation:: 559 560 sage: m = hmm.GaussianHiddenMarkovModel([[.1,.9],[.5,.5]], [(1,.5), (-1,3)], [.1,.9]) 561 sage: m.viterbi([-2,-1,.1,0.1]) 562 ([1, 1, 0, 1], -9.61823698847639...) 563 sage: m.viterbi([-2,-1,.1,0.3]) 564 ([1, 1, 1, 0], -9.5660236533785135) 565 """ 566 cdef TimeSeries _obs 567 if not isinstance(obs, TimeSeries): 568 _obs = TimeSeries(obs) 89 569 else: 90 return None 91 570 _obs = obs 92 571 93 cdef class GaussianHiddenMarkovModel(ContinuousHiddenMarkovModel): 94 r""" 95 Create a Gaussian Hidden Markov Model. The probability 96 distribution associated with each state is a Gaussian 97 distribution. 572 # The algorithm is the same as _viterbi above, except 573 # we take the logarithms of everything first, and add 574 # instead of multiply. 575 cdef Py_ssize_t t, T = _obs._length 576 cdef IntList state_sequence = IntList(T) 577 if T == 0: 578 return state_sequence, 0.0 98 579 99 GaussianHiddenMarkovModel(A, B, pi, name) 580 cdef int i, j, N = self.N 581 582 # delta[i] is the maximum of the probabilities over all 583 # paths ending in state i. 584 cdef TimeSeries delta = TimeSeries(N), delta_prev = TimeSeries(N) 585 586 # We view psi as an N x T matrix of ints. The quantity 587 # psi[N*t + j] 588 # is a most probable hidden state at time t-1, given 589 # that j is a most probable state at time j. 590 cdef IntList psi = IntList(N * T) # initialized to 0 by default 591 592 # Log Preprocessing 593 cdef TimeSeries A = self.A.log() 594 cdef TimeSeries pi = self.pi.log() 595 596 # Initialization: 597 for i in range(N): 598 delta._values[i] = pi._values[i] + log(self.probability_of(i, _obs._values[0])) 599 600 # Recursion: 601 cdef double mx, tmp, minus_inf = float('-inf') 602 cdef int index 603 604 for t in range(1, T): 605 delta_prev, delta = delta, delta_prev 606 for j in range(N): 607 # Compute delta_t[j] = max_i(delta_{t-1}(i) a_{i,j}) * b_j(_obs[t]) 608 mx = minus_inf 609 index = -1 610 for i in range(N): 611 tmp = delta_prev._values[i] + A._values[i*N+j] 612 if tmp > mx: 613 mx = tmp 614 index = i 615 delta._values[j] = mx + log(self.probability_of(j, _obs._values[t])) 616 psi._values[N*t + j] = index 617 618 # Termination: 619 mx, index = delta.max(index=True) 620 621 # Path (state sequence) backtracking: 622 state_sequence._values[T-1] = index 623 t = T-2 624 while t >= 0: 625 state_sequence._values[t] = psi._values[N*(t+1) + state_sequence._values[t+1]] 626 t -= 1 627 628 return state_sequence, mx 629 630 cdef TimeSeries _backward_scale_all(self, TimeSeries obs, TimeSeries scale): 631 """ 632 This function returns the matrix beta_t(i), and is used 633 internally as part of the Baum-Welch algorithm. 634 635 The quantity beta_t(i) is the probability of observing the 636 sequence obs[t+1:] assuming that the model is in state i at 637 time t. 638 639 INPUT: 640 641 - obs -- TimeSeries 642 - scale -- TimeSeries 643 644 OUTPUT: 645 646 - TimeSeries beta such that beta_t(i) = beta[t*N + i] 647 - scale is also changed by this function 648 """ 649 cdef Py_ssize_t t, T = obs._length 650 cdef int N = self.N, i, j 651 cdef double s 652 cdef TimeSeries beta = TimeSeries(N*T, initialize=False) 653 654 # 1. Initialization 655 for i in range(N): 656 beta._values[(T-1)*N + i] = 1 / scale._values[T-1] 657 658 # 2. Induction 659 t = T-2 660 while t >= 0: 661 for i in range(N): 662 s = 0 663 for j in range(N): 664 s += self.A._values[i*N+j] * \ 665 self.probability_of(j, obs._values[t+1]) * beta._values[(t+1)*N+j] 666 beta._values[t*N + i] = s/scale._values[t] 667 t -= 1 668 return beta 669 670 cdef _forward_scale_all(self, TimeSeries obs): 671 """ 672 Return scaled values alpha_t(i), the sequence of scalings, and 673 the log probability. 674 675 The quantity alpha_t(i) is the probability of observing the 676 sequence obs[:t+1] assuming that the model is in state i at 677 time t. 678 679 INPUT: 680 681 - obs -- TimeSeries 682 683 OUTPUT: 684 685 - TimeSeries alpha with alpha_t(i) = alpha[t*N + i] 686 - TimeSeries scale with scale[t] the scaling at step t 687 - float -- log_probability of the observation sequence 688 being produced by the model. 689 """ 690 cdef Py_ssize_t i, j, t, T = len(obs) 691 cdef int N = self.N 692 693 # The running some of the log probabilities 694 cdef double log_probability = 0, sum, a 695 696 cdef TimeSeries alpha = TimeSeries(N*T, initialize=False) 697 cdef TimeSeries scale = TimeSeries(T, initialize=False) 698 699 # Initialization 700 sum = 0 701 for i in range(self.N): 702 a = self.pi._values[i] * self.probability_of(i, obs._values[0]) 703 alpha._values[0*N + i] = a 704 sum += a 705 706 scale._values[0] = sum 707 log_probability = log(sum) 708 for i in range(self.N): 709 alpha._values[0*N + i] /= sum 710 711 # Induction 712 # The code below is just an optimized version of: 713 # alpha = (alpha * A).pairwise_product(B[O[t+1]]) 714 # alpha = alpha.scale(1/alpha.sum()) 715 # along with keeping track of the log of the scaling factor. 716 cdef double s 717 for t in range(1,T): 718 sum = 0 719 for j in range(N): 720 s = 0 721 for i in range(N): 722 s += alpha._values[(t-1)*N + i] * self.A._values[i*N + j] 723 a = s * self.probability_of(j, obs._values[t]) 724 alpha._values[t*N + j] = a 725 sum += a 726 727 log_probability += log(sum) 728 scale._values[t] = sum 729 for j in range(self.N): 730 alpha._values[t*N + j] /= sum 731 732 # Termination 733 return alpha, scale, log_probability 734 735 cdef TimeSeries _baum_welch_gamma(self, TimeSeries alpha, TimeSeries beta): 736 # TODO: factor out to hmm.pyx 737 cdef int j, N = self.N 738 cdef Py_ssize_t t, T = alpha._length//N 739 cdef TimeSeries gamma = TimeSeries(alpha._length, initialize=False) 740 cdef double denominator 741 for t in range(T): 742 denominator = 0 743 for j in range(N): 744 gamma._values[t*N + j] = alpha._values[t*N + j] * beta._values[t*N + j] 745 denominator += gamma._values[t*N + j] 746 for j in range(N): 747 gamma._values[t*N + j] /= denominator 748 return gamma 749 750 cdef TimeSeries _baum_welch_xi(self, TimeSeries alpha, TimeSeries beta, TimeSeries obs): 751 # TODO: factor out to hmm.pyx 752 """ 753 Return 3-dimensional array of xi values. 754 We have x[t,i,j] = x[t*N*N + i*N + j]. 755 """ 756 cdef int i, j, N = self.N 757 cdef double sum 758 cdef Py_ssize_t t, T = alpha._length//N 759 cdef TimeSeries xi = TimeSeries(T*N*N, initialize=False) 760 for t in range(T-1): 761 sum = 0.0 762 for i in range(N): 763 for j in range(N): 764 xi._values[t*N*N+i*N+j] = alpha._values[t*N+i]*beta._values[(t+1)*N+j]*\ 765 self.A._values[i*N+j] * self.probability_of(j, obs._values[t+1]) 766 sum += xi._values[t*N*N+i*N+j] 767 for i in range(N): 768 for j in range(N): 769 xi._values[t*N*N+i*N+j] /= sum 770 return xi 771 772 def baum_welch(self, obs, int max_iter=500, double log_likelihood_cutoff=1e-4, double eps=1e-12, bint fix_emissions=False): 773 """ 774 Given an observation sequence obs, improve this HMM using the 775 Baum-Welch algorithm to increase the probability of observing obs. 776 777 INPUT: 778 779 - obs -- a time series of emissions 780 781 - max_iter -- integer (default: 500) maximum number 782 of Baum-Welch steps to take 783 784 - log_likehood_cutoff -- positive float (default: 1e-4); 785 the minimal improvement in likelihood with respect to 786 the last iteration required to continue. Relative value 787 to log likelihood. 788 789 - eps -- very small positive float (default: 1e-12); used 790 during the iteration to replace numbers, e.g., standard 791 deviations that are 0 but are not allowed to be 0. 792 793 - fix_emissions -- bool (default: False); if True, do not 794 change emissions when updating 795 796 OUTPUT: 797 798 - changes the model in places, and returns the log 799 likelihood and number of iterations. 800 801 EXAMPLES:: 802 803 sage: m = hmm.GaussianHiddenMarkovModel([[.1,.9],[.5,.5]], [(1,.5), (-1,3)], [.1,.9]) 804 sage: m.log_likelihood([-2,-1,.1,0.1]) 805 -8.8582822159862751 806 sage: m.baum_welch([-2,-1,.1,0.1]) 807 (22.164539478647512, 8) 808 sage: m.log_likelihood([-2,-1,.1,0.1]) 809 22.164539478647512 810 sage: m 811 Gaussian Hidden Markov Model with 2 States 812 Transition matrix: 813 [ 1.0 1.99514384486e-21] 814 [ 0.499999997782 0.500000002218] 815 Emission parameters: 816 [(0.09999999999..., 1.48492424379e-06), (-1.4999999929, 0.500000010249)] 817 Initial probabilities: [0.0000, 1.0000] 818 819 We watch the log likelihoods of the model converge, step by step:: 820 821 sage: m = hmm.GaussianHiddenMarkovModel([[.1,.9],[.5,.5]], [(1,.5), (-1,3)], [.1,.9]) 822 sage: v = m.sample(10) 823 sage: stats.TimeSeries([m.baum_welch(v,max_iter=1)[0] for _ in range(len(v))]) 824 [-20.1167, -17.7611, -16.9814, -16.9364, -16.9314, -16.9309, -16.9309, -16.9309, -16.9309, -16.9309] 825 826 We illustrate fixing emissions:: 827 828 sage: m = hmm.GaussianHiddenMarkovModel([[.1,.9],[.9,.1]], [(1,2),(-1,.5)], [.3,.7]) 829 sage: set_random_seed(0); v = m.sample(100) 830 sage: m.baum_welch(v,fix_emissions=True) 831 (-164.72944548204..., 23) 832 sage: m.emission_parameters() 833 [(1.0, 2.0), (-1.0, 0.5)] 834 sage: m = hmm.GaussianHiddenMarkovModel([[.1,.9],[.9,.1]], [(1,2),(-1,.5)], [.3,.7]) 835 sage: m.baum_welch(v) 836 (-162.85437039799817, 49) 837 sage: m.emission_parameters() 838 [(1.27224191726, 2.37136875176), (-0.948617467518, 0.576236038512)] 839 """ 840 if not isinstance(obs, TimeSeries): 841 obs = TimeSeries(obs) 842 cdef TimeSeries _obs = obs 843 cdef TimeSeries alpha, beta, scale, gamma, xi 844 cdef double log_probability, log_probability0, log_probability_prev, delta 845 cdef int i, j, k, N, n_iterations 846 cdef Py_ssize_t t, T 847 cdef double denominator_A, numerator_A, denominator_B, numerator_mean, numerator_std 848 849 # Initialization 850 alpha, scale, log_probability0 = self._forward_scale_all(_obs) 851 if not isfinite(log_probability0): 852 return (0.0, 0) 853 log_probability = log_probability0 854 beta = self._backward_scale_all(_obs, scale) 855 gamma = self._baum_welch_gamma(alpha, beta) 856 xi = self._baum_welch_xi(alpha, beta, _obs) 857 log_probability_prev = log_probability 858 N = self.N 859 n_iterations = 0 860 T = len(_obs) 861 862 # Re-estimation 863 while True: 864 865 # Reestimate 866 for i in range(N): 867 if not gamma._values[0*N+i]: 868 gamma._values[0*N+i] = eps 869 elif not isfinite(gamma._values[0*N+i]): 870 util.normalize_probability_TimeSeries(self.pi, 0, self.pi._length) 871 raise RuntimeError, "impossible to compute gamma during reestimation" 872 self.pi._values[i] = gamma._values[0*N+i] 873 874 # Update the probabilities pi to define a valid discrete distribution 875 util.normalize_probability_TimeSeries(self.pi, 0, self.pi._length) 876 877 # Reestimate transition matrix and emission probabilities in 878 # each state. 879 for i in range(N): 880 # Compute the updated transition matrix 881 denominator_A = 0.0 882 for t in range(T-1): 883 denominator_A += gamma._values[t*N+i] 884 if not isnormal(denominator_A): 885 raise RuntimeError, "unable to re-estimate transition matrix" 886 for j in range(N): 887 numerator_A = 0.0 888 for t in range(T-1): 889 numerator_A += xi._values[t*N*N+i*N+j] 890 self.A._values[i*N+j] = numerator_A / denominator_A 891 892 # Rescale the i-th row of the transition matrix to be 893 # a valid stochastic matrix: 894 util.normalize_probability_TimeSeries(self.A, i*N, (i+1)*N) 895 896 if not fix_emissions: 897 denominator_B = denominator_A + gamma._values[(T-1)*N + i] 898 if not isnormal(denominator_B): 899 raise RuntimeError, "unable to re-estimate emission probabilities" 900 901 numerator_mean = 0.0 902 numerator_std = 0.0 903 for t in range(T): 904 numerator_mean += gamma._values[t*N + i] * _obs._values[t] 905 numerator_std += gamma._values[t*N + i] * \ 906 (_obs._values[t] - self.B._values[2*i])*(_obs._values[t] - self.B._values[2*i]) 907 # restimated mean 908 self.B._values[2*i] = numerator_mean / denominator_B 909 # restimated standard deviation 910 self.B._values[2*i+1] = sqrt(numerator_std / denominator_B) 911 if not self.B._values[2*i+1]: 912 self.B._values[2*i+1] = eps 913 self.probability_init() 914 915 n_iterations += 1 916 if n_iterations >= max_iter: break 917 918 # Initialization for next iteration 919 alpha, scale, log_probability0 = self._forward_scale_all(_obs) 920 if not isfinite(log_probability0): break 921 log_probability = log_probability0 922 beta = self._backward_scale_all(_obs, scale) 923 gamma = self._baum_welch_gamma(alpha, beta) 924 xi = self._baum_welch_xi(alpha, beta, _obs) 925 926 # Compute the difference between the log probability of 927 # two iterations. 928 delta = log_probability - log_probability_prev 929 log_probability_prev = log_probability 930 931 # If the log probability does not change by more than delta, 932 # then terminate 933 if delta >= 0 and delta <= log_likelihood_cutoff: 934 break 935 936 return log_probability, n_iterations 937 938 939 cdef class GaussianMixtureHiddenMarkovModel(GaussianHiddenMarkovModel): 940 """ 941 GaussianMixtureHiddenMarkovModel(A, B, pi) 942 943 Gaussian mixture Hidden Markov Model. 100 944 101 945 INPUT: 102 A -- matrix; the transition matrix (n x n) 103 B -- list of n pairs (mu, sigma) that define the 104 Gaussian distributions associated to each state, 105 where mu is the mean and sigma the standard deviation. 106 pi -- list of floats that sums to 1.0; these are 107 the initial probabilities of each hidden state 108 name -- (default: None); a string 109 normalize -- (optional; default=True) whether or not to 110 normalize the model so the probabilities add to 1 946 - A -- matrix; the N x N transition matrix 947 - B -- list of pairs (mu,sigma) that define the distributions 948 - pi -- initial state probabilities 111 949 112 EXAMPLES:113 Define the transition matrix:114 sage: A = [[0,1,0],[0.5,0,0.5],[0.3,0.3,0.4]]115 950 116 Parameters of the normal emission distributions in pairs of (mu, sigma): 117 sage: B = [(0,1), (-1,0.5), (1,0.2)] 951 EXAMPLES:: 118 952 119 The initial probabilities per state: 120 sage: pi = [1,0,0] 953 sage: A = [[0.5,0.5],[0.5,0.5]] 954 sage: B = [[(0.9,(0.0,1.0)), (0.1,(1,10000))],[(1,(1,1)), (0,(0,0.1))]] 955 sage: hmm.GaussianMixtureHiddenMarkovModel(A, B, [1,0]) 956 Gaussian Mixture Hidden Markov Model with 2 States 957 Transition matrix: 958 [0.5 0.5] 959 [0.5 0.5] 960 Emission parameters: 961 [0.9*N(0.0,1.0) + 0.1*N(1.0,10000.0), 1.0*N(1.0,1.0) + 0.0*N(0.0,0.1)] 962 Initial probabilities: [1.0000, 0.0000] 121 963 122 We create the continuous Gaussian hidden Markov model defined by $A,B,\pi$: 123 sage: hmm.GaussianHiddenMarkovModel(A, B, pi) 124 Gaussian Hidden Markov Model with 3 States 964 TESTS: 965 966 If a standard deviation is 0, it is normalized to be slightly bigger than 0.:: 967 968 sage: hmm.GaussianMixtureHiddenMarkovModel([[1]], [[(1,(0,0))]], [1]) 969 Gaussian Mixture Hidden Markov Model with 1 States 125 970 Transition matrix: 126 [0.0 1.0 0.0] 127 [0.5 0.0 0.5] 128 [0.3 0.3 0.4] 971 [1.0] 129 972 Emission parameters: 130 [(0.0, 1.0), (-1.0, 0.5), (1.0, 0.20000000000000001)] 131 Initial probabilities: [1.0, 0.0, 0.0] 973 [1.0*N(0.0,1e-08)] 974 Initial probabilities: [1.0000] 975 976 We test that number of emission distributions must be the same as the number of states:: 977 978 sage: hmm.GaussianMixtureHiddenMarkovModel([[1]], [], [1]) 979 Traceback (most recent call last): 980 ... 981 ValueError: number of GaussianMixtures must be the same as number of entries of pi 982 983 sage: hmm.GaussianMixtureHiddenMarkovModel([[1]], [[]], [1]) 984 Traceback (most recent call last): 985 ... 986 ValueError: must specify at least one component of the mixture model 987 988 We test that the number of initial probabilities must equal the number of states:: 989 990 sage: hmm.GaussianMixtureHiddenMarkovModel([[1]], [[]], [1,2]) 991 Traceback (most recent call last): 992 ... 993 ValueError: number of entries of transition matrix A must be the square of the number of entries of pi 132 994 """ 133 def __init__(self, A, B, pi=None, name=None, normalize=True): 995 996 cdef object mixture # mixture 997 998 def __init__(self, A, B, pi=None, bint normalize=True): 134 999 """ 135 EXAMPLES: 136 We make a very simple model: 137 sage: hmm.Gaussian HiddenMarkovModel([[1]], [(0,1)], [1], 'simple')138 Gaussian Hidden Markov Model simple with 1States1000 EXAMPLES:: 1001 1002 sage: hmm.GaussianMixtureHiddenMarkovModel([[.9,.1],[.4,.6]], [[(.4,(0,1)), (.6,(1,0.1))],[(1,(0,1))]], [.7,.3]) 1003 Gaussian Mixture Hidden Markov Model with 2 States 139 1004 Transition matrix: 140 [1.0] 1005 [0.9 0.1] 1006 [0.4 0.6] 141 1007 Emission parameters: 142 [(0.0, 1.0)] 143 Initial probabilities: [1.0] 1008 [0.4*N(0.0,1.0) + 0.6*N(1.0,0.1), 1.0*N(0.0,1.0)] 1009 Initial probabilities: [0.7000, 0.3000] 1010 """ 1011 self.pi = util.initial_probs_to_TimeSeries(pi, normalize) 1012 self.N = len(self.pi) 1013 self.A = util.state_matrix_to_TimeSeries(A, self.N, normalize) 1014 if self.N*self.N != len(self.A): 1015 raise ValueError, "number of entries of transition matrix A must be the square of the number of entries of pi" 144 1016 145 We test a degenerate case: 146 sage: hmm.GaussianHiddenMarkovModel([], [], [], 'simple') 147 Gaussian Hidden Markov Model simple with 0 States 148 Transition matrix: 149 [] 150 Emission parameters: 151 [] 152 Initial probabilities: [] 1017 self.mixture = [b if isinstance(b, GaussianMixtureDistribution) else \ 1018 GaussianMixtureDistribution([flatten(x) for x in b]) for b in B] 1019 if len(self.mixture) != self.N: 1020 raise ValueError, "number of GaussianMixtures must be the same as number of entries of pi" 153 1021 154 TESTS: 155 We test normalize=False: 156 sage: hmm.GaussianHiddenMarkovModel([[1,100],[0,1]], [(0,1),(0,2)], [1/2,1/2], normalize=False).transition_matrix() 157 [ 1.0 100.0] 158 [ 0.0 1.0] 1022 def __repr__(self): 159 1023 """ 160 ContinuousHiddenMarkovModel.__init__(self, A, B, pi, name=name)1024 Return string representation. 161 1025 162 # Set number of outputs. 163 self.m.M = 1 164 # Set the model type to continuous 165 self.m.model_type = GHMM_kContinuousHMM 166 # 1 transition matrix 167 self.m.cos = 1 168 # Set that no a prior model probabilities are set. 169 self.m.prior = -1 170 # Dimension is 1 171 self.m.dim = 1 172 self._initialize_state(pi) 173 self.m.class_change = NULL 174 self.initialized = True 175 if normalize: 176 self.normalize() 1026 EXAMPLES:: 177 1027 178 def _initialize_state(self, pi): 1028 sage: hmm.GaussianMixtureHiddenMarkovModel([[.9,.1],[.4,.6]], [[(.4,(0,1)), (.6,(1,0.1))],[(1,(0,1))]], [.7,.3]).__repr__() 1029 'Gaussian Mixture Hidden Markov Model with 2 States\nTransition matrix:\n[0.9 0.1]\n[0.4 0.6]\nEmission parameters:\n[0.4*N(0.0,1.0) + 0.6*N(1.0,0.1), 1.0*N(0.0,1.0)]\nInitial probabilities: [0.7000, 0.3000]' 179 1030 """ 180 Allocate and initialize states. 181 182 INPUT: 183 pi -- initial probabilities 184 185 All other inputs are set as self.A and self.B. 186 187 EXAMPLES: 188 This function is called implicitly during object creation. It should 189 never be called directly by the user, unless they want to LEAK MEMORY. 190 191 sage: m = hmm.GaussianHiddenMarkovModel([[1]], [(1,10)], [1]) # indirect test 192 """ 193 cdef ghmm_cstate* states = <ghmm_cstate*> safe_malloc(sizeof(ghmm_cstate) * self.m.N) 194 cdef ghmm_cstate* state 195 cdef ghmm_c_emission* e 196 cdef Py_ssize_t i, j, k 197 198 for i in range(self.m.N): 199 # Parameters of normal distribution 200 mu, sigma = self.B[i] 201 # Get a reference to the i-th state for convenience of the notation below. 202 state = &(states[i]) 203 state.M = 1 204 state.pi = pi[i] 205 state.desc = NULL 206 state.fix = 0 207 e = <ghmm_c_emission*> safe_malloc(sizeof(ghmm_c_emission)) 208 e.type = 0 # normal 209 e.dimension = 1 210 e.mean.val = mu 211 e.variance.val = sigma*sigma # variance! not standard deviation 212 # fixing of emissions is deactivated by default 213 e.fixed = 0 214 e.sigmacd = NULL 215 e.sigmainv = NULL 216 state.e = e 217 state.c = to_double_array([1.0]) 218 219 ######################################################### 220 # Initialize state transition data. 221 # NOTE: This code is similar to a block of code in hmm.pyx. 222 223 # Set "out" probabilities, i.e., the probabilities to 224 # transition to another hidden state from this state. 225 v = self.A[i] 226 k = self.m.N 227 state.out_states = k 228 state.out_id = <int*> safe_malloc(sizeof(int)*k) 229 state.out_a = ighmm_cmatrix_alloc(1, k) 230 for j in range(k): 231 state.out_id[j] = j 232 state.out_a[0][j] = v[j] 233 234 # Set "in" probabilities 235 v = self.A.column(i) 236 state.in_states = k 237 state.in_id = <int*> safe_malloc(sizeof(int)*k) 238 state.in_a = ighmm_cmatrix_alloc(1, k) 239 for j in range(k): 240 state.in_id[j] = j 241 state.in_a[0][j] = v[j] 242 243 ######################################################### 244 # Set states 245 self.m.s = states 1031 s = "Gaussian Mixture Hidden Markov Model with %s States"%self.N 1032 s += '\nTransition matrix:\n%s'%self.transition_matrix() 1033 s += '\nEmission parameters:\n%s'%self.emission_parameters() 1034 s += '\nInitial probabilities: %s'%self.initial_probabilities() 1035 return s 246 1036 247 1037 def __reduce__(self): 248 1038 """ 249 1039 Used in pickling. 250 1040 251 1041 EXAMPLES: 252 sage: m = hmm.GaussianHiddenMarkovModel([[1]], [(0,1)], [1], 'test') 253 sage: f,g = m.__reduce__() 254 sage: f(*g) == m 1042 sage: m = hmm.GaussianMixtureHiddenMarkovModel([[1]], [[(.4,(0,1)), (.6,(1,0.1))]], [1]) 1043 sage: loads(dumps(m)) == m 255 1044 True 256 1045 """ 257 return unpickle_gaussian_ hmm_v0, (self.transition_matrix(), self.emission_parameters(),258 self.initial_probabilities(), self.name())1046 return unpickle_gaussian_mixture_hmm_v1, \ 1047 (self.A, self.B, self.pi, self.mixture) 259 1048 260 def __dealloc__(self):261 """262 Dealloc the memory used by this Gaussian HMM, but only if the263 HMM was successfully initialized.264 265 EXAMPLES:266 sage: m = hmm.GaussianHiddenMarkovModel([[1.0]], [(0.0,1.0)], [1]) # implicit doctest267 sage: del m268 """269 if self.initialized:270 ghmm_cmodel_free(&self.m)271 272 def __copy__(self):273 """274 Return a copy of this Gaussian HMM.275 276 EXAMPLES:277 sage: m = hmm.GaussianHiddenMarkovModel([[0.4,0.6],[0.1,0.9]], [(0.0,1.0),(1,1)], [1,0], 'NAME')278 sage: copy(m)279 Gaussian Hidden Markov Model NAME with 2 States280 Transition matrix:281 [0.4 0.6]282 [0.1 0.9]283 Emission parameters:284 [(0.0, 1.0), (1.0, 1.0)]285 Initial probabilities: [1.0, 0.0]286 """287 288 return GaussianHiddenMarkovModel(self.transition_matrix(), self.emission_parameters(),289 self.initial_probabilities(), self.name())290 1049 291 1050 def __cmp__(self, other): 292 1051 """ 293 Compare two Gaussian HMM's. 1052 Compare self and other, which must both be GaussianMixtureHiddenMarkovModel's. 1053 1054 EXAMPLES:: 1055 1056 sage: m = hmm.GaussianMixtureHiddenMarkovModel([[1]], [[(.4,(0,1)), (.6,(1,0.1))]], [1]) 1057 sage: n = hmm.GaussianMixtureHiddenMarkovModel([[1]], [[(.5,(0,1)), (.5,(1,0.1))]], [1]) 1058 sage: m < n 1059 True 1060 sage: m == m 1061 True 1062 sage: n > m 1063 True 1064 sage: n < m 1065 False 1066 """ 1067 if not isinstance(other, GaussianMixtureHiddenMarkovModel): 1068 raise ValueError 1069 return cmp(self.__reduce__()[1], other.__reduce__()[1]) 1070 1071 def __getitem__(self, Py_ssize_t i): 1072 """ 1073 Return the Gaussian mixture distribution associated to the 1074 i-th state. 294 1075 295 1076 INPUT: 296 self, other -- objects; if other is not a Gaussian HMM compare types.297 OUTPUT:298 -1,0,1299 1077 300 The transition matrices are compared, then the emission 301 parameters, and the initial probabilities. The names are not 302 compared, so two GHMM's with the same parameters but different 303 names compare to be equal. 304 305 EXAMPLES: 306 sage: m = hmm.GaussianHiddenMarkovModel([[0.4,0.6],[0.1,0.9]], [(0.0,1.0),(1,1)], [1,2], "Test 1", normalize=False) 307 sage: m.__cmp__(m) 308 0 309 310 Note that the name doesn't matter: 311 sage: n = hmm.GaussianHiddenMarkovModel([[0.4,0.6],[0.1,0.9]], [(0.0,1.0),(1,1)], [1,2], "Test 2", normalize=False) 312 sage: m.__cmp__(n) 313 0 314 315 Normalizing fixes the initial probabilities, hence m and n are no longer equal. 316 sage: n.normalize() 317 sage: m.__cmp__(n) 318 1 319 """ 320 if not isinstance(other, GaussianHiddenMarkovModel): 321 return cmp(type(self), type(other)) 322 323 if self is other: return 0 # easy special case 324 325 cdef GaussianHiddenMarkovModel o = other 326 if self.m.N < o.m.N: 327 return -1 328 elif self.m.N > o.m.N: 329 return 1 330 cdef Py_ssize_t i, j 331 332 # The code below is somewhat long and tedious because I want it to be 333 # very fast. All it does is explicitly loop through the transition 334 # matrix, emission parameters and initial state probabilities checking 335 # that they agree, and if not returning -1 or 1. 336 # Compare transition matrices 337 for i from 0 <= i < self.m.N: 338 for j from 0 <= j < self.m.s[i].out_states: 339 if self.m.s[i].out_a[0][j] < o.m.s[i].out_a[0][j]: 340 return -1 341 elif self.m.s[i].out_a[0][j] > o.m.s[i].out_a[0][j]: 342 return 1 343 344 # Compare emissions parameters 345 for i from 0 <= i < self.m.N: 346 if self.m.s[i].e.mean.val < o.m.s[i].e.mean.val: 347 return -1 348 elif self.m.s[i].e.mean.val > o.m.s[i].e.mean.val: 349 return 1 350 if self.m.s[i].e.variance.val < o.m.s[i].e.variance.val: 351 return -1 352 elif self.m.s[i].e.variance.val > o.m.s[i].e.variance.val: 353 return 1 354 355 # Compare initial state probabilities 356 for 0 <= i < self.m.N: 357 if self.m.s[i].pi < o.m.s[i].pi: 358 return -1 359 elif self.m.s[i].pi > o.m.s[i].pi: 360 return 1 361 362 return 0 363 364 def fix_emissions(self, Py_ssize_t i, bint fixed=True): 365 """ 366 Sets the Gaussian emission parameters for the i-th state to be 367 either fixed or not fixed. If it is fixed, then running the 368 Baum-Welch algorithm will not change the emission parameters 369 for the i-th state. 370 371 INPUT: 372 i -- nonnegative integer < self.m.N 373 fixed -- bool 374 375 EXAMPLES: 376 We run Baum-Welch once without fixing the emission states: 377 sage: m = hmm.GaussianHiddenMarkovModel([[0.4,0.6],[0.1,0.9]], [(0.0,1.0),(1,1)], [1,0]) 378 sage: m.baum_welch([0,1]) 379 sage: m 380 Gaussian Hidden Markov Model with 2 States 381 Transition matrix: 382 [0.0 1.0] 383 [0.1 0.9] 384 Emission parameters: 385 [(0.0, 0.01), (1.0, 0.01)] 386 Initial probabilities: [1.0, 0.0] 387 388 Now we run Baum-Welch with the emission states fixed. Notice that they don't change. 389 sage: m = hmm.GaussianHiddenMarkovModel([[0.4,0.6],[0.1,0.9]], [(0.0,1.0),(1,1)], [1,0]) 390 sage: m.fix_emissions(0); m.fix_emissions(1) 391 sage: m.baum_welch([0,1]) 392 sage: m 393 Gaussian Hidden Markov Model with 2 States 394 Transition matrix: 395 [0.000368587006957 0.999631412993] 396 [ 0.1 0.9] 397 Emission parameters: 398 [(0.0, 1.0), (1.0, 1.0)] 399 Initial probabilities: [..., 0.0] 400 """ 401 if i < 0 or i >= self.m.N: 402 raise IndexError, "index out of range" 403 self.m.s[i].e.fixed = fixed 404 405 def __repr__(self): 406 """ 407 Return string representation of this Continuous HMM. 1078 - i -- integer 408 1079 409 1080 OUTPUT: 410 string411 1081 412 EXAMPLES: 413 sage: m = hmm.GaussianHiddenMarkovModel([[0.0,1.0,0],[0.5,0.0,0.5],[0.3,0.3,0.4]], [(0.0,1.0), (-1.0,0.5), (1.0,0.2)], [1,0,0]) 414 sage: m.__repr__() 415 'Gaussian Hidden Markov Model with 3 States\nTransition matrix:\n[0.0 1.0 0.0]\n[0.5 0.0 0.5]\n[0.3 0.3 0.4]\nEmission parameters:\n[(0.0, 1.0), (-1.0, 0.5), (1.0, 0.20000000000000001)]\nInitial probabilities: [1.0, 0.0, 0.0]' 1082 - a Gaussian mixture distribution object 1083 1084 EXAMPLES:: 1085 1086 sage: m = hmm.GaussianMixtureHiddenMarkovModel([[.9,.1],[.4,.6]], [[(.4,(0,1)), (.6,(1,0.1))],[(1,(0,1))]], [.7,.3]) 1087 sage: m[0] 1088 0.4*N(0.0,1.0) + 0.6*N(1.0,0.1) 1089 sage: m[1] 1090 1.0*N(0.0,1.0) 1091 1092 Negative indexing works:: 1093 1094 sage: m[-1] 1095 1.0*N(0.0,1.0) 1096 1097 Bounds are checked:: 1098 1099 sage: m[2] 1100 Traceback (most recent call last): 1101 ... 1102 IndexError: index out of range 1103 sage: m[-3] 1104 Traceback (most recent call last): 1105 ... 1106 IndexError: index out of range 416 1107 """ 417 s = "Gaussian Hidden Markov Model%s with %s States"%( 418 ' ' + self.m.name if self.m.name else '', 419 self.m.N) 420 s += '\nTransition matrix:\n%s'%self.transition_matrix() 421 s += '\nEmission parameters:\n%s'%self.emission_parameters() 422 s += '\nInitial probabilities: %s'%self.initial_probabilities() 423 return s 424 425 def initial_probabilities(self): 426 """ 427 Return the list of initial state probabilities. 428 429 OUTPUT: 430 list of floats 431 432 EXAMPLES: 433 sage: m = hmm.GaussianHiddenMarkovModel([[0.0,1.0,0],[0.5,0.0,0.5],[0.3,0.3,0.4]], [(0.0,1.0), (-1.0,0.5), (1.0,0.2)], [0.4,0.3,0.3]) 434 sage: m.initial_probabilities() 435 [0.40000000000000002, 0.29999999999999999, 0.29999999999999999] 436 """ 437 cdef Py_ssize_t i 438 return [self.m.s[i].pi for i in range(self.m.N)] 439 440 def transition_matrix(self): 441 """ 442 Return the hidden state transition matrix. 443 444 OUTPUT: 445 matrix whose rows give the transition probabilities of the 446 Hidden Markov Model states. 447 448 EXAMPLES: 449 sage: m = hmm.GaussianHiddenMarkovModel([[0.0,1.0,0],[0.5,0.0,0.5],[0.3,0.3,0.4]], [(0.0,1.0), (-1.0,0.5), (1.0,0.2)], [1,0,0]) 450 sage: m.transition_matrix() 451 [0.0 1.0 0.0] 452 [0.5 0.0 0.5] 453 [0.3 0.3 0.4] 454 """ 455 cdef Py_ssize_t i, j 456 # Update the state of the "immutable" matrix A, then return a reference to it. 457 for i from 0 <= i < self.m.N: 458 for j from 0 <= j < self.m.s[i].out_states: 459 self.A.set_unsafe_double(i,j,self.m.s[i].out_a[0][j]) 460 return self.A 1108 if i < 0: 1109 i += self.N 1110 if i < 0 or i >= self.N: 1111 raise IndexError, 'index out of range' 1112 return self.mixture[i] 461 1113 462 1114 def emission_parameters(self): 463 1115 """ 464 Return the emission parameters list.1116 Returns a list of all the emission distributions. 465 1117 466 1118 OUTPUT: 467 list of tuples (mu, sigma) that define Gaussian distributions associated to each state. 468 Here mu is the mean and sigma the standard deviation. 469 470 EXAMPLES: 471 sage: m = hmm.GaussianHiddenMarkovModel([[0.4,0.6],[0.1,0.9]], [(1.5,2),(-1,3)], [1,0], 'NAME') 1119 1120 - list of Gaussian mixtures 1121 1122 EXAMPLES:: 1123 1124 sage: m = hmm.GaussianMixtureHiddenMarkovModel([[.9,.1],[.4,.6]], [[(.4,(0,1)), (.6,(1,0.1))],[(1,(0,1))]], [.7,.3]) 472 1125 sage: m.emission_parameters() 473 [ (1.5, 2.0), (-1.0, 3.0)]1126 [0.4*N(0.0,1.0) + 0.6*N(1.0,0.1), 1.0*N(0.0,1.0)] 474 1127 """ 475 cdef Py_ssize_t i 476 return [(self.m.s[i].e.mean.val, sqrt(self.m.s[i].e.variance.val)) for i in range(self.m.N)] 1128 return list(self.mixture) 477 1129 478 def normalize(self):1130 cdef double random_sample(self, int state, randstate rstate): 479 1131 """ 480 Normalize the transition and emission probabilities, if481 applicable. This changes self in place.1132 Return a random sample from the normal distribution associated 1133 to the given state. 482 1134 483 EXAMPLES: 484 sage: m = hmm.GaussianHiddenMarkovModel([[1.0,1.2],[0,0.1]], [(0.0,1.0),(1,1)], [1,2]) 485 sage: m.normalize() 486 sage: m 487 Gaussian Hidden Markov Model with 2 States 488 Transition matrix: 489 [0.454545454545 0.545454545455] 490 [ 0.0 1.0] 491 Emission parameters: 492 [(0.0, 1.0), (1.0, 1.0)] 493 Initial probabilities: [0.33333333333333331, 0.66666666666666663] 494 """ 495 if ghmm_cmodel_normalize(self.m): 496 raise RuntimeError, "error normalizing model (note that model may still have been partly changed)" 497 498 def sample(self, long length, number=None): 499 """ 500 Return number samples from this HMM of given length. 1135 This is only used internally, and no bounds or other error 1136 checking is done, so calling this improperly can lead to seg 1137 faults. 501 1138 502 1139 INPUT: 503 length -- positive integer 504 1140 1141 - state -- integer 1142 - rstate -- randstate instance 1143 505 1144 OUTPUT: 506 if number is not given, return a single TimeSeries.507 if number is given, return a list of TimeSeries.508 1145 509 EXAMPLES: 510 sage: m = hmm.GaussianHiddenMarkovModel([[1]], [(0,1)], [1]) 511 sage: set_random_seed(0) 512 sage: m.sample(5,2) 513 [[-2.2808, -0.0699, 0.1858, 1.3624, 1.8252], 514 [0.0080, 0.1244, 0.5098, 0.9961, 0.4710]] 1146 - double 1147 """ 1148 cdef GaussianMixtureDistribution G = self.mixture[state] 1149 return G._sample(rstate) 515 1150 516 sage: m = hmm.GaussianHiddenMarkovModel([[1]], [(0,1)], [1]) 517 sage: set_random_seed(0) 518 sage: m.sample(5) 519 [-2.2808, -0.0699, 0.1858, 1.3624, 1.8252] 1151 cdef double probability_of(self, int state, double observation): 520 1152 """ 521 if number is None: 522 single = True 523 number = 1 524 else: 525 single = False 526 seed = random() 527 cdef ghmm_cseq *d = ghmm_cmodel_generate_sequences(self.m, seed, length, number, length) 528 cdef Py_ssize_t i, j 529 cdef TimeSeries T 530 ans = [] 531 for j from 0 <= j < number: 532 T = TimeSeries(length) 533 for i from 0 <= i < length: 534 T._values[i] = d.seq[j][i] 535 ans.append(T) 536 if single: 537 return ans[0] 538 return ans 539 # The line below would replace the above by something that returns a list of lists. 540 #return [[d.seq[j][i] for i in range(length)] for j in range(number)] 1153 Return the probability b_j(o) of see the given observation o 1154 (=observation) given that we're in the given state j (=state). 541 1155 542 543 544 545 #################################################################### 546 # HMM Problem 1 -- Computing likelihood: Given the parameter set 547 # lambda of an HMM model and an observation sequence O, determine 548 # the likelihood P(O | lambda). 549 #################################################################### 550 def log_likelihood(self, seq): 551 r""" 552 HMM Problem 1: Likelihood. 553 Computes sum over all sequence of seq_w * log( P ( O|lambda )). 1156 This is a continuous probability, so this really returns a 1157 number p such that the probability of a value in the interval 1158 [o,o+d] is p*d. 554 1159 555 1160 INPUT: 556 seq -- a single TimeSeries, or 557 a list of object z where z is either a TimeSeries 558 or a pair (TimeSeries, float), where float is a positive 559 weight. When the weight isn't given it defaults to 1. 1161 1162 - state -- integer 1163 - observation -- double 560 1164 561 1165 OUTPUT: 562 float -- the weight sum of logs of the probability of563 observing these sequences given the model564 1166 565 EXAMPLES:566 We create a very simple GHMM that generates random numbers with mean 0567 and standard deviation 1.568 sage: m = hmm.GaussianHiddenMarkovModel([[1]], [(0,1)], [1])1167 - double 1168 """ 1169 cdef GaussianMixtureDistribution G = self.mixture[state] 1170 return G.prob(observation) 569 1171 570 We compute the log probability of a certain series: 571 sage: m.log_likelihood(finance.TimeSeries([1,0,1,1])) 572 -5.1757541328186907 1172 # TODO: factor out to hmm.pyx 1173 cdef TimeSeries _baum_welch_gamma(self, TimeSeries alpha, TimeSeries beta): 1174 """ 1175 Returns gamma and gamma_all, where gamma_all contains the 1176 (double-precision float) numbers gamma_t(j,m), and gamma 1177 contains just the gamma_t(j). 1178 """ 1179 cdef int j, N = self.N 1180 cdef Py_ssize_t t, T = alpha._length//N 1181 cdef TimeSeries gamma = TimeSeries(alpha._length, initialize=False) 1182 cdef double denominator 1183 for t in range(T): 1184 denominator = 0 1185 for j in range(N): 1186 gamma._values[t*N + j] = alpha._values[t*N + j] * beta._values[t*N + j] 1187 denominator += gamma._values[t*N + j] 1188 for j in range(N): 1189 gamma._values[t*N + j] /= denominator 573 1190 574 We compute the log probability of another much less likely sequence: 575 sage: m.log_likelihood(finance.TimeSeries([1,0,1,20])) 576 -204.6757541328187 1191 return gamma 577 1192 578 We compute weighted sum of log probabilities of two sequences: 579 sage: m.log_likelihood([ ([1,0,1,1], 10), ([1,0,1,20], 0.1) ]) 580 -72.2251167414687... 581 582 We verify that the above weight computation is right. 583 sage: 10 * m.log_likelihood([1,0,1,1]) + 0.1 * m.log_likelihood([1,0,1,20]) 584 -72.2251167414688 585 586 We generate two normally distributed sequences and see that the one 587 with mean 0 and standard deviation 1 is vastly more likely given 588 the model than the one with mean 10 and standard deviation 1. 589 sage: set_random_seed(0) 590 sage: m.log_likelihood(finance.TimeSeries(100).randomize('normal',0,1)) 591 -129.87209711900121 592 sage: m.log_likelihood(finance.TimeSeries(100).randomize('normal',10,1)) 593 -5010.151947016132 594 595 However, if we change the model then the situation reverses: 596 sage: m = hmm.GaussianHiddenMarkovModel([[1]], [(10,1)], [1]) 597 sage: m.log_likelihood(finance.TimeSeries(100).randomize('normal',10,1)) 598 -149.68737352182211 599 sage: m.log_likelihood(finance.TimeSeries(100).randomize('normal',0,1)) 600 -5275.30829407876... 1193 cdef TimeSeries _baum_welch_mixed_gamma(self, TimeSeries alpha, TimeSeries beta, 1194 TimeSeries obs, int j): 601 1195 """ 602 603 cdef double log_p 604 cdef ghmm_cseq* sqd 605 try: 606 sqd = to_cseq(seq) 607 except RuntimeError: 608 # no sequences 609 return float(0) 610 cdef int ret = ghmm_cmodel_likelihood(self.m, sqd, &log_p) 611 ghmm_cseq_free(&sqd) 612 if ret == -1: 613 raise RuntimeError, "unable to compute log likelihood" 614 return log_p 615 616 617 #################################################################### 618 # HMM Problem 2 -- Decoding: Given the complete parameter set that 619 # defines the model and an observation sequence seq, determine the 620 # best hidden sequence Q. Use the Viterbi algorithm. 621 #################################################################### 622 def viterbi(self, seq): 623 """ 624 HMM Problem 2: Decoding. Determine a hidden sequence of 625 states that is most likely to produce the given sequence seq 626 of observations. 1196 Let gamma_t(j,m) be the m-component (in the mixture) of the 1197 probability of being in state j at time t, given the 1198 observation sequence. This function outputs a TimeSeries v 1199 such that v[m*T + t] gives gamma_t(j, m) where T is the number 1200 of time steps. 627 1201 628 1202 INPUT: 629 seq -- TimeSeries 630 1203 1204 - alpha -- TimeSeries 1205 - beta -- TimeSeries 1206 - obs -- TimeSeries 1207 - j -- int 1208 631 1209 OUTPUT: 632 list -- a most probable sequence of hidden states, i.e., the633 Viterbi path.634 float -- log of the probability that the sequence of hidden635 states actually produced the given sequence seq.636 1210 637 EXAMPLES: 638 We find the optimal state sequence for a given model. 639 sage: m = hmm.GaussianHiddenMarkovModel([[0.5,0.5],[0.5,0.5]], [(0,1),(10,1)], [0.5,0.5]) 640 sage: m.viterbi([0,1,10,10,1]) 641 ([0, 0, 1, 1, 0], -9.0604285688230899) 1211 - TimeSeries 642 1212 """ 643 cdef double log_p 644 if not isinstance(seq, TimeSeries): 645 seq = TimeSeries(seq) 646 cdef TimeSeries T = seq 647 cdef int* path = ghmm_cmodel_viterbi(self.m, T._values, T._length, &log_p) 648 if not path: 649 raise RuntimeError, "sequence can't be built from model" 650 cdef Py_ssize_t i 651 v = [path[i] for i in range(T._length)] 652 sage_free(path) 653 return v, log_p 1213 cdef int i, k, m, N = self.N 1214 cdef Py_ssize_t t, T = alpha._length//N 654 1215 655 #################################################################### 656 # HMM Problem 3 -- Learning: Given an observation sequence O and 657 # the set of states in the HMM, improve the HMM to increase the 658 # probability of observing O. 659 #################################################################### 660 def baum_welch(self, training_seqs, max_iter=10000, log_likelihood_cutoff=0.00001): 1216 cdef double numer, alpha_minus, P, s, prob 1217 cdef GaussianMixtureDistribution G = self.mixture[j] 1218 cdef int M = len(G) 1219 cdef TimeSeries mixed_gamma = TimeSeries(T*M) 1220 1221 for t in range(T): 1222 prob = self.probability_of(j, obs._values[t]) 1223 if prob == 0: 1224 # If the probability of observing obs[t] in state j is 0, then 1225 # all of the m-mixture components *have* to automatically be 0, 1226 # since prob is the sum of those and they are all nonnegative. 1227 for m in range(M): 1228 mixed_gamma._values[m*T + t] = 0 1229 else: 1230 # Compute the denominator we used when scaling gamma. 1231 # The key thing is that this is consistent between 1232 # gamma and mixed_gamma. 1233 P = 0 1234 for k in range(N): 1235 P += alpha._values[t*N+k]*beta._values[t*N+k] 1236 1237 # Divide out the total probability, so we can multiply back in 1238 # the m-components of the probability. 1239 alpha_minus = alpha._values[t*N + j] / prob 1240 for m in range(M): 1241 numer = alpha_minus * G.prob_m(obs._values[t], m) * beta._values[t*N + j] 1242 mixed_gamma._values[m*T + t] = numer / P 1243 1244 return mixed_gamma 1245 1246 def baum_welch(self, obs, int max_iter=1000, double log_likelihood_cutoff=1e-12, double eps=1e-12, bint fix_emissions=False): 661 1247 """ 662 HMM Problem 3: Learning. Given an observation sequence O and 663 the set of states in the HMM, improve the HMM using the 664 Baum-Welch algorithm to increase the probability of observing O. 1248 Given an observation sequence obs, improve this HMM using the 1249 Baum-Welch algorithm to increase the probability of observing obs. 665 1250 666 1251 INPUT: 667 training_seqs -- a list of lists of emission symbols, where all sequences of 668 length 0 are ignored; or, a list of pairs 669 (sample_sequence, weight), 670 where sample_sequence is a list or TimeSeries, and weight is 671 a positive real number. 672 max_iter -- integer or None (default: 10000) maximum number 673 of Baum-Welch steps to take 674 log_likehood_cutoff -- positive float or None (default: 0.00001); 675 the minimal improvement in likelihood 676 with respect to the last iteration required to 677 continue. Relative value to log likelihood 1252 1253 - obs -- a time series of emissions 1254 - max_iter -- integer (default: 1000) maximum number 1255 of Baum-Welch steps to take 1256 - log_likehood_cutoff -- positive float (default: 1e-12); 1257 the minimal improvement in likelihood with respect to 1258 the last iteration required to continue. Relative value 1259 to log likelihood. 1260 - eps -- very small positive float (default: 1e-12); used 1261 during the iteration to replace numbers, e.g., standard 1262 deviations that are 0 but are not allowed to be 0. 1263 - fix_emissions -- bool (default: False); if True, do not 1264 change emissions when updating 678 1265 679 1266 OUTPUT: 680 changes the model in places, or raises a RuntimError681 exception on error682 1267 683 EXAMPLES: 684 We train a very simple model: 685 sage: m = hmm.GaussianHiddenMarkovModel([[1]], [(0,1)], [1]) 686 sage: m.baum_welch([1,1,1,1]) 1268 - changes the model in places, and returns the log 1269 likelihood and number of iterations. 687 1270 688 Notice that after training the mean in the emission parameters changes 689 form 0 to 1. 1271 EXAMPLES:: 1272 1273 sage: m = hmm.GaussianMixtureHiddenMarkovModel([[.9,.1],[.4,.6]], [[(.4,(0,1)), (.6,(1,0.1))],[(1,(0,1))]], [.7,.3]) 1274 sage: set_random_seed(0); v = m.sample(10); v 1275 [0.3576, -0.9365, 0.9449, -0.6957, 1.0217, 0.9644, 0.9987, -0.5950, -1.0219, 0.6477] 1276 sage: m.log_likelihood(v) 1277 -8.31408655939536... 1278 sage: m.baum_welch(v) 1279 (2.18905068682301..., 15) 1280 sage: m.log_likelihood(v) 1281 2.18905068682301... 690 1282 sage: m 691 Gaussian Hidden Markov Model with 1States1283 Gaussian Mixture Hidden Markov Model with 2 States 692 1284 Transition matrix: 693 [1.0] 1285 [ 0.874636333977 0.125363666023] 1286 [ 1.0 1.45168520229e-40] 694 1287 Emission parameters: 695 [ (1.0, 0.01)]696 Initial probabilities: [ 1.0]1288 [0.500161629343*N(-0.812298726239,0.173329026744) + 0.499838370657*N(0.982433690378,0.029719932009), 1.0*N(0.503260056832,0.145881515324)] 1289 Initial probabilities: [0.0000, 1.0000] 697 1290 698 We train using a list of lists: 699 sage: m = hmm.GaussianHiddenMarkovModel([[1,0],[0,1]], [(0,1),(0,2)], [1/2,1/2]) 700 sage: m.baum_welch([[1,2], [3,2]]) 701 sage: m 702 Gaussian Hidden Markov Model with 2 States 703 Transition matrix: 704 [1.0 0.0] 705 [0.0 1.0] 706 Emission parameters: 707 [(1.94653953598434..., 0.705082962992410...), (2.02081569132933..., 0.70680033099099593)] 708 Initial probabilities: [0.28024729110782..., 0.71975270889217...] 1291 We illustrate fixing all emissions:: 709 1292 710 We train the same model, but waiting one of the lists more than the other. 711 sage: m = hmm.GaussianHiddenMarkovModel([[1,0],[0,1]], [(0,1),(0,2)], [1/2,1/2]) 712 sage: m.baum_welch([([1,2,],10), ([3,2],1)]) 713 sage: m 714 Gaussian Hidden Markov Model with 2 States 715 Transition matrix: 716 [1.0 0.0] 717 [0.0 1.0] 718 Emission parameters: 719 [(1.58517861517798..., 0.572645807401051...), (1.594503566606473..., 0.57928632238916...)] 720 Initial probabilities: [0.385468570528119..., 0.61453142947188...] 721 1293 sage: m = hmm.GaussianMixtureHiddenMarkovModel([[.9,.1],[.4,.6]], [[(.4,(0,1)), (.6,(1,0.1))],[(1,(0,1))]], [.7,.3]) 1294 sage: set_random_seed(0); v = m.sample(10) 1295 sage: m.baum_welch(v, fix_emissions=True) 1296 (-7.5865685899788922, 36) 1297 sage: m.emission_parameters() 1298 [0.4*N(0.0,1.0) + 0.6*N(1.0,0.1), 1.0*N(0.0,1.0)] 1299 """ 1300 if not isinstance(obs, TimeSeries): 1301 obs = TimeSeries(obs) 1302 cdef TimeSeries _obs = obs 1303 cdef TimeSeries alpha, beta, scale, gamma, mixed_gamma, mixed_gamma_m, xi 1304 cdef double log_probability, log_probability0, log_probability_prev, delta 1305 cdef int i, j, k, m, N, n_iterations 1306 cdef Py_ssize_t t, T 1307 cdef double denominator_A, numerator_A, denominator_B, numerator_mean, numerator_std, \ 1308 numerator_c, c, mu, std, numer, denom, new_mu, new_std, new_c, s 1309 cdef GaussianMixtureDistribution G 722 1310 723 Training sequences of length 0 are gracefully ignored: 724 sage: m.baum_welch([]) 725 sage: m.baum_welch([([],1)]) 726 """ 727 cdef ghmm_cmodel_baum_welch_context cs 1311 # Initialization 1312 alpha, scale, log_probability0 = self._forward_scale_all(_obs) 1313 if not isfinite(log_probability0): 1314 return (0.0, 0) 1315 log_probability = log_probability0 1316 beta = self._backward_scale_all(_obs, scale) 1317 gamma = self._baum_welch_gamma(alpha, beta) 1318 xi = self._baum_welch_xi(alpha, beta, _obs) 1319 log_probability_prev = log_probability 1320 N = self.N 1321 n_iterations = 0 1322 T = len(_obs) 728 1323 729 cs.smo = self.m 730 try: 731 cs.sqd = to_cseq(training_seqs) 732 except RuntimeError: 733 # No nonempty sequences 734 return 735 cs.logp = <double*> safe_malloc(sizeof(double)) 736 cs.eps = log_likelihood_cutoff 737 cs.max_iter = max_iter 738 739 if ghmm_cmodel_baum_welch(&cs): 740 raise RuntimeError, "error running Baum-Welch algorithm" 741 ghmm_cseq_free(&cs.sqd) 742 743 744 cdef ghmm_cseq* to_cseq(seq) except NULL: 745 """ 746 Return a pointer to a ghmm_cseq C struct. 1324 # Re-estimation 1325 while True: 747 1326 748 All empty sequences are ignored. If there are no nonempty 749 sequences a RuntimeError is raised, since GHMM doesn't treat 750 this degenerate case well. If there are any nonpositive 751 weights, then a ValueError is raised. 752 """ 753 if isinstance(seq, list) and len(seq) > 0 and not isinstance(seq[0], (list, tuple, TimeSeries)): 754 seq = TimeSeries(seq) 755 if isinstance(seq, TimeSeries): 756 seq = [(seq,float(1))] 757 cdef Py_ssize_t i, n 758 for i in range(len(seq)): 759 z = seq[i] 760 if isinstance(z, tuple) and len(z) == 2: 761 if isinstance(z[0],TimeSeries): 762 z = (z[0], float(z[1])) 763 else: 764 z = (TimeSeries(z[0]), float(z[1])) 765 else: 766 if isinstance(z, TimeSeries): 767 z = (z, float(1)) 768 else: 769 z = (TimeSeries(z), float(1)) 770 seq[i] = z 771 seq = [x for x in seq if len(x[0]) > 0] 1327 # Reestimate frequency of state i in time t=0 1328 for i in range(N): 1329 if not gamma._values[0*N+i]: 1330 gamma._values[0*N+i] = eps 1331 elif not isnormal(gamma._values[0*N+i]): 1332 util.normalize_probability_TimeSeries(self.pi, 0, self.pi._length) 1333 raise RuntimeError, "impossible to compute gamma during reestimation" 1334 self.pi._values[i] = gamma._values[0*N+i] 772 1335 773 n = len(seq) 774 if n == 0: 775 raise RuntimeError, "there must be at least one nonempty sequence" 1336 # Update the probabilities pi to define a valid discrete distribution 1337 util.normalize_probability_TimeSeries(self.pi, 0, self.pi._length) 776 1338 777 for _, w in seq: 778 if w <= 0: 779 raise ValueError, "each weight must be positive" 780 781 cdef ghmm_cseq* sqd = <ghmm_cseq*>safe_malloc(sizeof(ghmm_cseq)) 782 sqd.seq = <double**>safe_malloc(sizeof(double*) * n) 783 sqd.seq_len = to_int_array([len(v) for v,_ in seq]) 784 sqd.seq_id = to_double_array([0]*n) 785 sqd.seq_label = to_int_array([]) # obsolete but no choice 786 weights = [w for _, w in seq] 787 sqd.seq_w = to_double_array(weights) 788 sqd.seq_number = n 789 sqd.total_w = sum(weights) 790 sqd.dim = 1 791 sqd.flags = 0 792 sqd.capacity = n 1339 # Reestimate transition matrix and emission probabilities in 1340 # each state. 1341 for i in range(N): 1342 # Reestimate the state transition matrix 1343 denominator_A = 0.0 1344 for t in range(T-1): 1345 denominator_A += gamma._values[t*N+i] 1346 if not isnormal(denominator_A): 1347 raise RuntimeError, "unable to re-estimate pi (1)" 1348 for j in range(N): 1349 numerator_A = 0.0 1350 for t in range(T-1): 1351 numerator_A += xi._values[t*N*N+i*N+j] 1352 self.A._values[i*N+j] = numerator_A / denominator_A 793 1353 794 cdef TimeSeries T 795 for i from 0 <= i < len(seq): 796 T = seq[i][0] 797 sqd.seq[i] = <double*>safe_malloc(sizeof(double) * T._length) 798 memcpy(sqd.seq[i], T._values , sizeof(double)*T._length) 1354 # Rescale the i-th row of the transition matrix to be 1355 # a valid stochastic matrix: 1356 util.normalize_probability_TimeSeries(self.A, i*N, (i+1)*N) 799 1357 800 return sqd 1358 ######################################################################## 1359 # Re-estimate the emission probabilities 1360 ######################################################################## 1361 G = self.mixture[i] 1362 if not fix_emissions and not G.is_fixed(): 1363 mixed_gamma = self._baum_welch_mixed_gamma(alpha, beta, _obs, i) 1364 new_G = [] 1365 for m in range(len(G)): 1366 if G.fixed._values[m]: 1367 new_G.append(G[m]) 1368 continue 801 1369 1370 # Compute re-estimated mu_{j,m} 1371 numer = 0 1372 denom = 0 1373 for t in range(T): 1374 numer += mixed_gamma._values[m*T + t] * _obs._values[t] 1375 denom += mixed_gamma._values[m*T + t] 1376 new_mu = numer / denom 1377 1378 # Compute re-estimated standard deviation 1379 numer = 0 1380 mu = G[m][1] 1381 for t in range(T): 1382 numer += mixed_gamma._values[m*T + t] * \ 1383 (_obs._values[t] - mu)*(_obs._values[t] - mu) 1384 1385 new_std = sqrt(numer / denom) 1386 1387 # Compute re-estimated weighting coefficient 1388 new_c = denom 1389 s = 0 1390 for t in range(T): 1391 s += gamma._values[t*N + i] 1392 new_c /= s 1393 1394 new_G.append((new_c,new_mu,new_std)) 1395 1396 self.mixture[i] = GaussianMixtureDistribution(new_G) 1397 1398 n_iterations += 1 1399 if n_iterations >= max_iter: break 1400 1401 ######################################################################## 1402 # Initialization for next iteration 1403 ######################################################################## 1404 alpha, scale, log_probability0 = self._forward_scale_all(_obs) 1405 if not isfinite(log_probability0): break 1406 log_probability = log_probability0 1407 beta = self._backward_scale_all(_obs, scale) 1408 gamma = self._baum_welch_gamma(alpha, beta) 1409 xi = self._baum_welch_xi(alpha, beta, _obs) 1410 1411 # Compute the difference between the log probability of 1412 # two iterations. 1413 delta = log_probability - log_probability_prev 1414 log_probability_prev = log_probability 1415 1416 # If the log probability does not improve by more than 1417 # delta, then terminate 1418 if delta >= 0 and delta <= log_likelihood_cutoff: 1419 break 1420 1421 return log_probability, n_iterations 1422 1423 1424 ################################################## 1425 # For Unpickling 1426 ################################################## 1427 1428 # We keep the _v0 function for backwards compatible. 802 1429 def unpickle_gaussian_hmm_v0(A, B, pi, name): 803 1430 """ 804 EXAMPLES: 805 sage: m = hmm.GaussianHiddenMarkovModel([[1]], [(0,1)], [1], 'test') 806 sage: loads(dumps(m)) == m 807 True 1431 EXAMPLES:: 1432 1433 sage: m = hmm.GaussianHiddenMarkovModel([[1]], [(0,1)], [1]) 808 1434 sage: sage.stats.hmm.chmm.unpickle_gaussian_hmm_v0(m.transition_matrix(), m.emission_parameters(), m.initial_probabilities(), 'test') 809 Gaussian Hidden Markov Model testwith 1 States1435 Gaussian Hidden Markov Model with 1 States 810 1436 Transition matrix: 811 1437 [1.0] 812 1438 Emission parameters: 813 1439 [(0.0, 1.0)] 814 Initial probabilities: [1.0 ]1440 Initial probabilities: [1.0000] 815 1441 """ 816 return GaussianHiddenMarkovModel(A,B,pi ,name)1442 return GaussianHiddenMarkovModel(A,B,pi) 817 1443 818 1444 819 cdef class GaussianMixtureHiddenMarkovModel(GaussianHiddenMarkovModel):1445 def unpickle_gaussian_hmm_v1(A, B, pi, prob, n_out): 820 1446 """ 821 GaussianMixtureHiddenMarkovModel(A, B, pi, name)1447 EXAMPLES:: 822 1448 823 INPUT:824 A -- matrix; the transition matrix (n x n)825 B -- list of lists of pairs (w, (mu, sigma)) that define the826 Gaussian mixture associated to each state, where w is827 the weight, mu is the mean and sigma the standard828 deviation.829 pi -- list of floats that sums to 1.0; these are830 the initial probabilities of each hidden state831 name -- (default: None); a string832 normalize -- (optional; default=True) whether or not to normalize833 the model so the probabilities add to 11449 sage: m = hmm.GaussianHiddenMarkovModel([[1]], [(0,1)], [1]) 1450 sage: loads(dumps(m)) == m # indirect test 1451 True 1452 """ 1453 cdef GaussianHiddenMarkovModel m = PY_NEW(GaussianHiddenMarkovModel) 1454 m.A = A 1455 m.B = B 1456 m.pi = pi 1457 m.prob = prob 1458 m.n_out = n_out 1459 return m 834 1460 835 EXAMPLES: 836 sage: A = [[0.5,0.5],[0.5,0.5]] 837 sage: B = [[(0.5,(0.0,1.0)), (0.1,(1,10000))],[(1,(1,1)), (0,(0,0.1))]] 838 sage: pi = [1,0] 839 sage: hmm.GaussianMixtureHiddenMarkovModel(A, B, pi) 840 Gaussian Hidden Markov Model with 2 States 841 Transition matrix: 842 [0.5 0.5] 843 [0.5 0.5] 844 Emission parameters: 845 [[(1.0, (0.0, 1.0)), (0.10000000000000001, (1.0, 10000.0))], [(1.0, (1.0, 1.0)), (0.0, (0.0, 0.10000000000000001))]] 846 Initial probabilities: [1.0, 0.0] 847 1461 def unpickle_gaussian_mixture_hmm_v1(A, B, pi, mixture): 1462 """ 1463 EXAMPLES:: 848 1464 849 TESTS: 850 We test that standard deviations must be positive: 851 sage: m = hmm.GaussianMixtureHiddenMarkovModel([[1]], [[(1,(0,0))]], [1]) 852 Traceback (most recent call last): 853 ... 854 ValueError: sigma must be positive (if weight is nonzero) 1465 sage: m = hmm.GaussianMixtureHiddenMarkovModel([[1]], [[(.4,(0,1)), (.6,(1,0.1))]], [1]) 1466 sage: loads(dumps(m)) == m # indirect test 1467 True 1468 """ 1469 cdef GaussianMixtureHiddenMarkovModel m = PY_NEW(GaussianMixtureHiddenMarkovModel) 1470 m.A = A 1471 m.B = B 1472 m.pi = pi 1473 m.mixture = mixture 1474 return m 855 1475 856 We test that number of mixtures must be positive:857 sage: m = hmm.GaussianMixtureHiddenMarkovModel([[1]], [[]], [1])858 Traceback (most recent call last):859 ...860 ValueError: number of Gaussian mixtures must be positive861 """862 def __init__(self, A, B, pi=None, name=None, normalize=True):863 """864 EXAMPLES:865 sage: hmm.GaussianMixtureHiddenMarkovModel([[1]], [[(1,(0,1))]], [1], name='simple')866 Gaussian Hidden Markov Model simple with 1 States867 Transition matrix:868 [1.0]869 Emission parameters:870 [[(1.0, (0.0, 1.0))]]871 Initial probabilities: [1.0]872 """873 # Turn B into a list of lists874 B = [flatten(x) for x in B]875 m = max([len(x) for x in B])876 if m == 0:877 raise ValueError, "number of Gaussian mixtures must be positive"878 B = [x + [0]*(m-len(x)) for x in B]879 GaussianHiddenMarkovModel.__init__(self, A, B, pi, name=name)880 self.m.M = m//3881 # Set number of outputs.882 883 def _initialize_state(self, pi):884 """885 Allocate and initialize states.886 887 INPUT:888 pi -- initial probabilities889 890 All other inputs are set as self.A and self.B.891 892 EXAMPLES:893 This function is called implicitly during object creation. It should894 never be called directly by the user, unless they want to LEAK MEMORY.895 896 sage: m = hmm.GaussianMixtureHiddenMarkovModel([[1]], [[(1,(1,10))]], [1]) # indirect test897 """898 cdef ghmm_cstate* states = <ghmm_cstate*> safe_malloc(sizeof(ghmm_cstate) * self.m.N)899 cdef ghmm_cstate* state900 cdef ghmm_c_emission* e901 cdef Py_ssize_t i, j, k, M, n902 903 for i in range(self.m.N):904 # Parameters of Gaussian distributions905 v = self.B[i]906 M = len(v)//3 # number of distinct Gaussians907 908 # Get a reference to the i-th state for convenience of the notation below.909 state = &(states[i])910 state.M = M911 state.pi = pi[i]912 state.desc = NULL913 state.fix = 0914 e = <ghmm_c_emission*> safe_malloc(sizeof(ghmm_c_emission)*M)915 weights = []916 917 for n in range(M):918 e[n].type = 0 # Gaussian919 e[n].dimension = 1920 mu = v[n*3+1]921 sigma = v[n*3+2]922 if sigma <= 0 and v[n*3]:923 raise ValueError, "sigma must be positive (if weight is nonzero)"924 weights.append( v[n*3] )925 e[n].mean.val = mu926 e[n].variance.val = sigma*sigma # variance! not standard deviation927 928 # fixing of emissions is deactivated by default929 e[n].fixed = 0930 e[n].sigmacd = NULL931 e[n].sigmainv = NULL932 933 state.e = e934 state.c = to_double_array(weights)935 936 #########################################################937 # Initialize state transition data.938 # NOTE: This code is similar to a block of code in hmm.pyx.939 940 # Set "out" probabilities, i.e., the probabilities to941 # transition to another hidden state from this state.942 v = self.A[i]943 k = self.m.N944 state.out_states = k945 state.out_id = <int*> safe_malloc(sizeof(int)*k)946 state.out_a = ighmm_cmatrix_alloc(1, k)947 for j in range(k):948 state.out_id[j] = j949 state.out_a[0][j] = v[j]950 951 # Set "in" probabilities952 v = self.A.column(i)953 state.in_states = k954 state.in_id = <int*> safe_malloc(sizeof(int)*k)955 state.in_a = ighmm_cmatrix_alloc(1, k)956 for j in range(k):957 state.in_id[j] = j958 state.in_a[0][j] = v[j]959 960 #########################################################961 # Set states962 self.m.s = states963 964 def emission_parameters(self):965 """966 Return the emission parameters list.967 968 OUTPUT:969 list of lists of tuples (weight, (mu, sigma))970 971 EXAMPLES:972 sage: m = hmm.GaussianMixtureHiddenMarkovModel([[0.5,0.5],[0.5,0.5]], [[(0.5,(0.0,1.0)), (0.1,(1,10000))],[(1,(1,1))]], [1,0])973 sage: m.emission_parameters()974 [[(1.0, (0.0, 1.0)), (0.10000000000000001, (1.0, 10000.0))], [(1.0, (1.0, 1.0)), (0.0, (0.0, 0.0))]]975 """976 cdef Py_ssize_t i,j977 978 return [[(self.m.s[i].c[j], (self.m.s[i].e[j].mean.val, sqrt(self.m.s[i].e[j].variance.val)))979 for j in range(self.m.s[i].M)] for i in range(self.m.N)]980 981 982 983 def fix_emissions(self, Py_ssize_t i, Py_ssize_t j, bint fixed=True):984 """985 Sets the j-th Gaussian of the emission parameters for the i-th986 state to be either fixed or not fixed. If it is fixed, then987 running the Baum-Welch algorithm will not change the emission988 parameters for the i-th state.989 990 INPUT:991 i -- nonnegative integer < self.m.N992 j -- nonnegative integer < self.m.M993 fixed -- bool994 995 EXAMPLES:996 We run Baum-Welch once without fixing the emission states:997 sage: m = hmm.GaussianHiddenMarkovModel([[0.4,0.6],[0.1,0.9]], [(0.0,1.0),(1,1)], [1,0])998 sage: m.baum_welch([0,1])999 sage: m1000 Gaussian Hidden Markov Model with 2 States1001 Transition matrix:1002 [0.0 1.0]1003 [0.1 0.9]1004 Emission parameters:1005 [(0.0, 0.01), (1.0, 0.01)]1006 Initial probabilities: [1.0, 0.0]1007 1008 Now we run Baum-Welch with the emission states fixed. Notice that they don't change.1009 sage: m = hmm.GaussianHiddenMarkovModel([[0.4,0.6],[0.1,0.9]], [(0.0,1.0),(1,1)], [1,0])1010 sage: m.fix_emissions(0); m.fix_emissions(1)1011 sage: m.baum_welch([0,1])1012 sage: m1013 Gaussian Hidden Markov Model with 2 States1014 Transition matrix:1015 [0.000368587006957 0.999631412993]1016 [ 0.1 0.9]1017 Emission parameters:1018 [(0.0, 1.0), (1.0, 1.0)]1019 Initial probabilities: [..., 0.0]1020 """1021 if i < 0 or i >= self.m.N:1022 raise IndexError, "index i out of range"1023 if j < 0 or j >= self.m.M:1024 raise IndexError, "index j out of range"1025 self.m.s[i].e[j].fixed = fixed -
new file sage/stats/hmm/distributions.pxd
diff --git a/sage/stats/hmm/distributions.pxd b/sage/stats/hmm/distributions.pxd new file mode 100644
- + 1 ############################################################################# 2 # Copyright (C) 2010 William Stein <wstein@gmail.com> 3 # Distributed under the terms of the GNU General Public License (GPL) 4 # The full text of the GPL is available at: 5 # http://www.gnu.org/licenses/ 6 ############################################################################# 7 8 from sage.finance.time_series cimport TimeSeries 9 from sage.stats.intlist cimport IntList 10 from sage.misc.randstate cimport randstate 11 12 cdef class Distribution: 13 pass 14 15 cdef class DiscreteDistribution(Distribution): 16 cdef object v 17 18 cdef class GaussianDistribution(Distribution): 19 pass 20 21 cdef class GaussianMixtureDistribution(Distribution): 22 cdef TimeSeries c0, c1, param 23 cdef IntList fixed 24 25 cdef double _sample(self, randstate rstate) 26 cpdef double prob(self, double x) 27 cpdef double prob_m(self, double x, int m) 28 cpdef is_fixed(self, i=?) 29 30 31 -
new file sage/stats/hmm/distributions.pyx
diff --git a/sage/stats/hmm/distributions.pyx b/sage/stats/hmm/distributions.pyx new file mode 100644
- + 1 """ 2 Distributions used in implementing Hidden Markov Models 3 4 These distribution classes are designed specifically for HMM's and not 5 for general use in statistics. For example, they have fixed or 6 non-fixed status, which only make sense relative to being used in a 7 hidden Markov model. 8 9 AUTHOR: 10 11 - William Stein, 2010-03 12 """ 13 14 ############################################################################# 15 # Copyright (C) 2010 William Stein <wstein@gmail.com> 16 # Distributed under the terms of the GNU General Public License (GPL) 17 # The full text of the GPL is available at: 18 # http://www.gnu.org/licenses/ 19 ############################################################################# 20 21 include "../../ext/stdsage.pxi" 22 23 cdef extern from "math.h": 24 double exp(double) 25 double log(double) 26 double sqrt(double) 27 28 import math 29 cdef double sqrt2pi = sqrt(2*math.pi) 30 31 from sage.misc.randstate cimport current_randstate, randstate 32 from sage.finance.time_series cimport TimeSeries 33 34 35 36 cdef double random_normal(double mean, double std, randstate rstate): 37 """ 38 Return a floating point number chosen from the normal distribution 39 with given mean and standard deviation, using the given randstate. 40 The computation uses the box muller algorithm. 41 42 INPUT: 43 44 - mean -- float; the mean 45 - std -- float; the standard deviation 46 - rstate -- randstate; the random number generator state 47 48 OUTPUT: 49 50 - double 51 """ 52 # Ported from http://users.tkk.fi/~nbeijar/soft/terrain/source_o2/boxmuller.c 53 # This the box muller algorithm. 54 # Client code can get the current random state from: 55 # cdef randstate rstate = current_randstate() 56 cdef double x1, x2, w, y1, y2 57 while True: 58 x1 = 2*rstate.c_rand_double() - 1 59 x2 = 2*rstate.c_rand_double() - 1 60 w = x1*x1 + x2*x2 61 if w < 1: break 62 w = sqrt( (-2*log(w))/w ) 63 y1 = x1 * w 64 return mean + y1*std 65 66 # Abstract base class for distributions used for hidden Markov models. 67 68 cdef class Distribution: 69 """ 70 A distribution. 71 """ 72 def sample(self, n=None): 73 """ 74 Return either a single sample (the default) or n samples from 75 this probability distribution. 76 77 INPUT: 78 79 - n -- None or a positive integer 80 81 OUTPUT: 82 83 - a single sample if n is 1; otherwise many samples 84 85 EXAMPLES:: 86 87 This method must be defined in a derived class:: 88 89 sage: sage.stats.hmm.distributions.Distribution().sample() 90 Traceback (most recent call last): 91 ... 92 NotImplementedError 93 """ 94 raise NotImplementedError 95 96 def prob(self, x): 97 """ 98 The probability density function evaluated at x. 99 100 INPUT: 101 102 - x -- object 103 104 OUTPUT: 105 106 - float 107 108 EXAMPLES:: 109 110 This method must be defined in a derived class:: 111 112 sage: sage.stats.hmm.distributions.Distribution().prob(0) 113 Traceback (most recent call last): 114 ... 115 NotImplementedError 116 """ 117 raise NotImplementedError 118 119 def plot(self, *args, **kwds): 120 """ 121 Return a plot of the probability density function. 122 123 INPUT: 124 125 - args and kwds, passed to the Sage plot function 126 127 OUTPUT: 128 129 - a Graphics object 130 131 EXAMPLES:: 132 133 sage: P = hmm.GaussianMixtureDistribution([(.2,-10,.5),(.6,1,1),(.2,20,.5)]) 134 sage: P.plot(-10,30) 135 """ 136 from sage.plot.all import plot 137 return plot(self.prob, *args, **kwds) 138 139 cdef class GaussianMixtureDistribution(Distribution): 140 """ 141 A probability distribution defined by taking a weighted linear 142 combination of Gaussian distributions. 143 144 EXAMPLES:: 145 146 sage: P = hmm.GaussianMixtureDistribution([(.3,1,2),(.7,-1,1)]); P 147 0.3*N(1.0,2.0) + 0.7*N(-1.0,1.0) 148 sage: P[0] 149 (0.29999999999999999, 1.0, 2.0) 150 sage: P.is_fixed() 151 False 152 sage: P.fix(1) 153 sage: P.is_fixed(0) 154 False 155 sage: P.is_fixed(1) 156 True 157 sage: P.unfix(1) 158 sage: P.is_fixed(1) 159 False 160 """ 161 def __init__(self, B, eps=1e-8, bint normalize=True): 162 """ 163 INPUT: 164 165 - `B` -- a list of triples `(c_i, mean_i, std_i)`, where 166 the `c_i` and `std_i` are positive and the sum of the 167 `c_i` is `1`. 168 169 - eps -- positive real number; any standard deviation in B 170 less than eps is replaced by eps. 171 172 - normalize -- if True, ensure that the c_i are nonnegative 173 174 EXAMPLES:: 175 176 sage: hmm.GaussianMixtureDistribution([(.3,1,2),(.7,-1,1)]) 177 0.3*N(1.0,2.0) + 0.7*N(-1.0,1.0) 178 sage: hmm.GaussianMixtureDistribution([(1,-1,0)], eps=1e-3) 179 1.0*N(-1.0,0.001) 180 """ 181 B = [[c if c>=0 else 0, mu, std if std>0 else eps] for c,mu,std in B] 182 if len(B) == 0: 183 raise ValueError, "must specify at least one component of the mixture model" 184 cdef double s 185 if normalize: 186 s = sum([a[0] for a in B]) 187 if s != 1: 188 if s == 0: 189 s = 1.0/len(B) 190 for a in B: 191 a[0] = s 192 else: 193 for a in B: 194 a[0] /= s 195 self.c0 = TimeSeries([c/(sqrt2pi*std) for c,_,std in B]) 196 self.c1 = TimeSeries([-1.0/(2*std*std) for _,_,std in B]) 197 self.param = TimeSeries(sum([list(x) for x in B],[])) 198 self.fixed = IntList(self.c0._length) 199 200 def __getitem__(self, Py_ssize_t i): 201 """ 202 Returns triple (coefficient, mu, std). 203 204 INPUT: 205 206 - i -- integer 207 208 OUTPUT: 209 210 - triple of floats 211 212 EXAMPLES:: 213 214 sage: P = hmm.GaussianMixtureDistribution([(.2,-10,.5),(.6,1,1),(.2,20,.5)]) 215 sage: P[0] 216 (0.20000000000000001, -10.0, 0.5) 217 sage: P[2] 218 (0.20000000000000001, 20.0, 0.5) 219 sage: [-1] 220 [-1] 221 sage: P[-1] 222 (0.20000000000000001, 20.0, 0.5) 223 sage: P[3] 224 Traceback (most recent call last): 225 ... 226 IndexError: index out of range 227 sage: P[-4] 228 Traceback (most recent call last): 229 ... 230 IndexError: index out of range 231 """ 232 if i < 0: i += self.param._length//3 233 if i < 0 or i >= self.param._length//3: raise IndexError, "index out of range" 234 return self.param._values[3*i], self.param._values[3*i+1], self.param._values[3*i+2] 235 236 def __reduce__(self): 237 """ 238 Used in pickling. 239 240 EXAMPLES:: 241 242 sage: G = sage.stats.hmm.distributions.GaussianMixtureDistribution([(.1,1,2), (.9,0,1)]) 243 sage: loads(dumps(G)) == G 244 True 245 """ 246 return unpickle_gaussian_mixture_distribution_v1, ( 247 self.c0, self.c1, self.param, self.fixed) 248 249 def __cmp__(self, other): 250 """ 251 EXAMPLES:: 252 253 sage: G = sage.stats.hmm.distributions.GaussianMixtureDistribution([(.1,1,2), (.9,0,1)]) 254 sage: H = sage.stats.hmm.distributions.GaussianMixtureDistribution([(.3,1,2), (.7,1,5)]) 255 sage: G < H 256 True 257 sage: H > G 258 True 259 sage: G == H 260 False 261 sage: G == G 262 True 263 """ 264 if not isinstance(other, GaussianMixtureDistribution): 265 raise ValueError 266 return cmp(self.__reduce__()[1], other.__reduce__()[1]) 267 268 def __len__(self): 269 """ 270 Return the number of components of this GaussianMixtureDistribution. 271 272 EXAMPLES:: 273 274 sage: len(hmm.GaussianMixtureDistribution([(.2,-10,.5),(.6,1,1),(.2,20,.5)])) 275 3 276 """ 277 return self.c0._length 278 279 cpdef is_fixed(self, i=None): 280 """ 281 Return whether or not this GaussianMixtureDistribution is 282 fixed when using Baum-Welch to update the corresponding HMM. 283 284 INPUT: 285 286 - i - None (default) or integer; if given, only return 287 whether the i-th component is fixed 288 289 EXAMPLES:: 290 291 sage: P = hmm.GaussianMixtureDistribution([(.2,-10,.5),(.6,1,1),(.2,20,.5)]) 292 sage: P.is_fixed() 293 False 294 sage: P.is_fixed(0) 295 False 296 sage: P.fix(0); P.is_fixed() 297 False 298 sage: P.is_fixed(0) 299 True 300 sage: P.fix(); P.is_fixed() 301 True 302 """ 303 if i is None: 304 return bool(self.fixed.prod()) 305 else: 306 return bool(self.fixed[i]) 307 308 def fix(self, i=None): 309 """ 310 Set that this GaussianMixtureDistribution (or its ith 311 component) is fixed when using Baum-Welch to update 312 the corresponding HMM. 313 314 INPUT: 315 316 - i - None (default) or integer; if given, only fix the 317 i-th component 318 319 EXAMPLES:: 320 321 sage: P = hmm.GaussianMixtureDistribution([(.2,-10,.5),(.6,1,1),(.2,20,.5)]) 322 sage: P.fix(1); P.is_fixed() 323 False 324 sage: P.is_fixed(1) 325 True 326 sage: P.fix(); P.is_fixed() 327 True 328 """ 329 cdef int j 330 if i is None: 331 for j in range(self.c0._length): 332 self.fixed[j] = 1 333 else: 334 self.fixed[i] = 1 335 336 def unfix(self, i=None): 337 """ 338 Set that this GaussianMixtureDistribution (or its ith 339 component) is not fixed when using Baum-Welch to update the 340 corresponding HMM. 341 342 INPUT: 343 344 - i - None (default) or integer; if given, only fix the 345 i-th component 346 347 EXAMPLES:: 348 349 sage: P = hmm.GaussianMixtureDistribution([(.2,-10,.5),(.6,1,1),(.2,20,.5)]) 350 sage: P.fix(1); P.is_fixed(1) 351 True 352 sage: P.unfix(1); P.is_fixed(1) 353 False 354 sage: P.fix(); P.is_fixed() 355 True 356 sage: P.unfix(); P.is_fixed() 357 False 358 359 """ 360 cdef int j 361 if i is None: 362 for j in range(self.c0._length): 363 self.fixed[j] = 0 364 else: 365 self.fixed[i] = 0 366 367 368 def __repr__(self): 369 """ 370 Return string representation of this mixed Gaussian distribution. 371 372 EXAMPLES:: 373 374 sage: hmm.GaussianMixtureDistribution([(.2,-10,.5),(.6,1,1),(.2,20,.5)]).__repr__() 375 '0.2*N(-10.0,0.5) + 0.6*N(1.0,1.0) + 0.2*N(20.0,0.5)' 376 """ 377 return ' + '.join(["%s*N(%s,%s)"%x for x in self]) 378 379 def sample(self, n=None): 380 """ 381 Return a single sample from this distribution (by default), or 382 if n>1, return a TimeSeries of samples. 383 384 INPUT: 385 386 - n -- integer or None (default: None) 387 388 OUTPUT: 389 390 - float if n is None (default); otherwise a TimeSeries 391 392 EXAMPLES:: 393 394 sage: P = hmm.GaussianMixtureDistribution([(.2,-10,.5),(.6,1,1),(.2,20,.5)]) 395 sage: P.sample() 396 19.658243610875129 397 sage: P.sample(1) 398 [-10.4683] 399 sage: P.sample(5) 400 [-0.1688, -10.3479, 1.6812, 20.1083, -9.9801] 401 sage: P.sample(0) 402 [] 403 sage: P.sample(-3) 404 Traceback (most recent call last): 405 ... 406 ValueError: n must be nonnegative 407 """ 408 cdef randstate rstate = current_randstate() 409 cdef Py_ssize_t i 410 cdef TimeSeries T 411 if n is None: 412 return self._sample(rstate) 413 else: 414 _n = n 415 if _n < 0: 416 raise ValueError, "n must be nonnegative" 417 T = TimeSeries(_n) 418 for i in range(_n): 419 T._values[i] = self._sample(rstate) 420 return T 421 422 cdef double _sample(self, randstate rstate): 423 """ 424 Used internally to compute a sample from this distribution quickly. 425 426 INPUT: 427 428 - rstate -- a randstate object 429 430 OUTPUT: 431 432 - double 433 """ 434 cdef double accum, r 435 cdef int n 436 accum = 0 437 r = rstate.c_rand_double() 438 439 # See the remark in hmm.pyx about using GSL to remove this 440 # silly way of sampling from a discrete distribution. 441 for n in range(self.c0._length): 442 accum += self.param._values[3*n] 443 if r <= accum: 444 return random_normal(self.param._values[3*n+1], self.param._values[3*n+2], rstate) 445 raise RuntimeError, "invalid probability distribution" 446 447 cpdef double prob(self, double x): 448 """ 449 Return the probability of x. 450 451 Since this is a continuous distribution, this is defined to be 452 the limit of the p's such that the probability of [x,x+h] is p*h. 453 454 INPUT: 455 456 - x -- float 457 458 OUTPUT: 459 460 - float 461 462 EXAMPLES:: 463 464 sage: P = hmm.GaussianMixtureDistribution([(.2,-10,.5),(.6,1,1),(.2,20,.5)]) 465 sage: P.prob(.5) 466 0.21123919605857971 467 sage: P.prob(-100) 468 0.0 469 sage: P.prob(20) 470 0.1595769121605731 471 """ 472 # The tricky-looking code below is a fast version of this: 473 # return sum([c/(sqrt(2*math.pi)*std) * \ 474 # exp(-(x-mean)*(x-mean)/(2*std*std)) for 475 # c, mean, std in self.B]) 476 cdef double s=0, mu 477 cdef int n 478 for n in range(self.c0._length): 479 mu = self.param._values[3*n+1] 480 s += self.c0._values[n]*exp((x-mu)*(x-mu)*self.c1._values[n]) 481 return s 482 483 cpdef double prob_m(self, double x, int m): 484 """ 485 Return the probability of x using just the m-th summand. 486 487 INPUT: 488 489 - x -- float 490 - m -- integer 491 492 OUTPUT: 493 494 - float 495 496 EXAMPLES:: 497 498 sage: P = hmm.GaussianMixtureDistribution([(.2,-10,.5),(.6,1,1),(.2,20,.5)]) 499 sage: P.prob_m(.5, 0) 500 2.760811768050888...e-97 501 sage: P.prob_m(.5, 1) 502 0.21123919605857971 503 sage: P.prob_m(.5, 2) 504 0.0 505 """ 506 cdef double s, mu 507 if m < 0 or m >= self.param._length//3: 508 raise IndexError, "index out of range" 509 mu = self.param._values[3*m+1] 510 return self.c0._values[m]*exp((x-mu)*(x-mu)*self.c1._values[m]) 511 512 def unpickle_gaussian_mixture_distribution_v1(TimeSeries c0, TimeSeries c1, 513 TimeSeries param, IntList fixed): 514 """ 515 Used in unpickling GaussianMixtureDistribution's. 516 517 EXAMPLES:: 518 519 sage: P = hmm.GaussianMixtureDistribution([(.2,-10,.5),(.6,1,1),(.2,20,.5)]) 520 sage: loads(dumps(P)) == P # indirect doctest 521 True 522 """ 523 cdef GaussianMixtureDistribution G = PY_NEW(GaussianMixtureDistribution) 524 G.c0 = c0 525 G.c1 = c1 526 G.param = param 527 G.fixed = fixed 528 return G -
sage/stats/hmm/hmm.pxd
diff --git a/sage/stats/hmm/hmm.pxd b/sage/stats/hmm/hmm.pxd
a b 1 1 ############################################################################# 2 # Copyright (C) 2008 William Stein <wstein@gmail.com> 3 # Distributed under the terms of the GNU General Public License (GPL), 4 # version 2 or any later version at your option. 2 # Copyright (C) 2010 William Stein <wstein@gmail.com> 3 # Distributed under the terms of the GNU General Public License (GPL) v2+. 5 4 # The full text of the GPL is available at: 6 5 # http://www.gnu.org/licenses/ 7 6 ############################################################################# 8 7 9 from sage.matrix.matrix_real_double_dense cimport Matrix_real_double_dense10 8 11 cdef extern from "ghmm/ghmm.h": 12 cdef int GHMM_kSilentStates 13 cdef int GHMM_kHigherOrderEmissions 14 cdef int GHMM_kDiscreteHMM 15 cdef int GHMM_EPS_ITER_BW 16 cdef int GHMM_MAX_ITER_BW 17 18 cdef extern from "ghmm/model.h": 19 ctypedef struct ghmm_dstate: 20 # Initial probability 21 double pi 22 # Output probability 23 double *b 24 25 # ID's of the following states 26 int *out_id 27 # ID's of the previous states 28 int *in_id 29 30 # Transition probabilities to successor states. 31 double *out_a 32 # Transition probabilities from predecessor states. 33 double *in_a 34 35 # Number of successor states 36 int out_states 37 # Number of precursor states 38 int in_states 39 # if fix == 1 then output probabilities b stay fixed during the training 40 int fix 41 # Contains a description of the state (null terminated utf-8) 42 char * desc 43 # x coordinate position for graph representation plotting 44 int xPosition 45 # y coordinate position for graph representation plotting 46 int yPosition 47 48 ctypedef struct ghmm_dbackground: 49 pass 50 51 ctypedef struct ghmm_alphabet: 52 pass 53 54 # Discrete HMM's. 55 ctypedef struct ghmm_dmodel: 56 # Number of states 57 int N 58 # Number of outputs 59 int M 60 # Vector of states 61 ghmm_dstate *s 62 # Contains bit flags for varios model extensions such as 63 # kSilentStates, kTiedEmissions (see ghmm.h for a complete list) 64 int model_type 65 # The a priori probability for the model. 66 # A value of -1 indicates that no prior is defined. 67 # Note: this is not to be confused with priors on emission distributions. 68 double prior 69 # An arbitrary name for the model (null terminated utf-8) 70 char * name 71 72 # Flag variables for each state indicating whether it is emitting or not. 73 # Note: silent != NULL iff (model_type & kSilentStates) == 1 74 int *silent 75 76 # Int variable for the maximum level of higher order emissions. 77 int maxorder 78 79 # saves the history of emissions as int, 80 # the nth-last emission is (emission_history * |alphabet|^n+1) % |alphabet| 81 int emission_history 82 83 # Flag variables for each state indicating whether the states emissions 84 # are tied to another state. Groups of tied states are represented 85 # by their tie group leader (the lowest numbered member of the group). 86 # tied_to[s] == kUntied : s is not a tied state 87 # tied_to[s] == s : s is a tie group leader 88 # tied_to[t] == s : t is tied to state s (t>s) 89 # Note: tied_to != NULL iff (model_type & kTiedEmissions) != 0 */ 90 int *tied_to 91 92 # Note: State store order information of the emissions. 93 # Classical HMMS have emission order 0; that is, the emission probability 94 # is conditioned only on the state emitting the symbol. 95 # For higher order emissions, the emission are conditioned on the state s 96 # as well as the previous emission_order observed symbols. 97 # The emissions are stored in the state's usual double* b. 98 # Note: order != NULL iff (model_type & kHigherOrderEmissions) != 0 */ 99 int * order 100 101 # ghmm_dbackground is a pointer to a 102 # ghmm_dbackground structure, which holds (essentially) an 103 # array of background distributions (which are just vectors of floating 104 # point numbers like state.b). 105 # For each state the array background_id indicates which of the background 106 # distributions to use in parameter estimation. A value of kNoBackgroundDistribution 107 # indicates that none should be used. 108 # Note: background_id != NULL iff (model_type & kHasBackgroundDistributions) != 0 109 int *background_id 110 ghmm_dbackground *bp 111 112 # Topological ordering of silent states 113 # Condition: topo_order != NULL iff (model_type & kSilentStates) != 0 114 int *topo_order 115 int topo_order_length 116 117 # pow_lookup is a array of precomputed powers 118 # It contains in the i-th entry M (alphabet size) to the power of i 119 # The last entry is maxorder+1. 120 int *pow_lookup 121 122 # Store for each state a class label. Limits the possibly state sequence 123 # Note: label != NULL iff (model_type & kLabeledStates) != 0 124 int* label 125 ghmm_alphabet* label_alphabet 126 ghmm_alphabet* alphabet 127 128 129 # Discrete sequences 130 ctypedef struct ghmm_dseq: 131 # sequence array. sequence[i] [j] = j-th symbol of i-th seq 132 int **seq 133 # matrix of state ids, can be used to save the viterbi path during sequence generation. 134 # ATTENTION: is NOT allocated by ghmm_dseq_calloc 135 int **states 136 # array of sequence length 137 int *seq_len 138 # array of state path lengths 139 int *states_len 140 # array of sequence IDs 141 double *seq_id 142 # positive sequence weights. default is 1 = no weight 143 double *seq_w 144 # total number of sequences 145 long seq_number 146 # reserved space for sequences is always >= seq_number 147 long capacity 148 # sum of sequence weights 149 double total_w 150 # matrix of state labels corresponding to seq 151 int **state_labels 152 # number of labels for each sequence 153 int *state_labels_len 154 # internal flags 155 unsigned int flags 156 157 ghmm_dseq *ghmm_dmodel_label_generate_sequences( 158 ghmm_dmodel * mo, int seed, int global_len, long seq_number, int Tmax) 159 ghmm_dseq *ghmm_dmodel_generate_sequences(ghmm_dmodel* mo, int seed, int global_len, 160 long seq_number, int Tmax) 161 int ghmm_dseq_free (ghmm_dseq ** sq) 162 ghmm_dseq *ghmm_dseq_calloc (long seq_number) 163 164 int ghmm_dmodel_normalize(ghmm_dmodel *m) 165 int ghmm_dmodel_free(ghmm_dmodel **m) 166 167 # Problem 1: Likelihood 168 int ghmm_dmodel_logp(ghmm_dmodel *m, int *O, int len, double *log_p) 169 170 # Problem 2: Decoding 171 int *ghmm_dmodel_viterbi (ghmm_dmodel *m, int *O, int len, int *pathlen, double *log_p) 172 173 # Problem 3: Learning 174 int ghmm_dmodel_baum_welch (ghmm_dmodel *m, ghmm_dseq *sq) 175 int ghmm_dmodel_baum_welch_nstep (ghmm_dmodel *m, ghmm_dseq *sq, 176 int nsteps, double log_likelihood_cutoff) 9 from sage.finance.time_series cimport TimeSeries 177 10 178 11 cdef class HiddenMarkovModel: 179 cdef Matrix_real_double_dense A, B180 cdef listpi12 cdef int N 13 cdef TimeSeries A, pi 181 14 182 cdef class DiscreteHiddenMarkovModel(HiddenMarkovModel):183 cdef ghmm_dmodel* m184 cdef bint initialized185 cdef object _emission_symbols, _emission_symbols_dict -
sage/stats/hmm/hmm.pyx
diff --git a/sage/stats/hmm/hmm.pyx b/sage/stats/hmm/hmm.pyx
a b 1 r"""nodoctest 1 """ 2 2 Hidden Markov Models 3 3 4 AUTHOR: William Stein 4 This is a complete pure-Cython optimized implementation of Hidden 5 Markov Models. It fully supports Discrete, Gaussian, and Mixed 6 Gaussian emissions. 7 8 The best references for the basic HMM algorithms implemented here are: 9 10 - Tapas Kanungo's "Hidden Markov Models" 11 12 - Jackson's HMM tutorial: 13 http://personal.ee.surrey.ac.uk/Personal/P.Jackson/tutorial/ 14 15 LICENSE: Some of the code in this file is based on reading Kanungo's 16 GPLv2+ implementation of discrete HMM's, hence the present code must 17 be licensed with a GPLv2+ compatible license. 18 19 AUTHOR: 20 21 - William Stein, 2010-03 5 22 """ 6 23 7 24 ############################################################################# 8 # Copyright (C) 2008 William Stein <wstein@gmail.com> 9 # Distributed under the terms of the GNU General Public License (GPL), 10 # version 2 or any later version at your option. 25 # Copyright (C) 2010 William Stein <wstein@gmail.com> 26 # Distributed under the terms of the GNU General Public License (GPL) v2+. 11 27 # The full text of the GPL is available at: 12 28 # http://www.gnu.org/licenses/ 13 29 ############################################################################# 14 30 15 import math 31 include "../../ext/stdsage.pxi" 32 include "../../ext/cdefs.pxi" 33 include "../../ext/interrupt.pxi" 16 34 35 cdef extern from "math.h": 36 double log(double) 37 38 from sage.finance.time_series cimport TimeSeries 39 from sage.stats.intlist cimport IntList 17 40 from sage.matrix.all import is_Matrix, matrix 18 from sage.rings.all import RDF 19 from sage.misc.randstate import random 41 from sage.misc.randstate cimport current_randstate, randstate 20 42 21 include "../../ext/stdsage.pxi" 43 from util cimport HMM_Util 22 44 23 include "misc.pxi" 45 cdef HMM_Util util = HMM_Util() 46 47 ########################################### 24 48 25 49 cdef class HiddenMarkovModel: 26 50 """ 27 51 Abstract base class for all Hidden Markov Models. 52 """ 53 def initial_probabilities(self): 54 """ 55 Return the initial probabilities, which as a TimeSeries of 56 length N, where N is the number of states of the Markov model. 28 57 29 INPUT: 30 A -- matrix or list 31 B -- matrix or list 32 pi -- list of floats 58 EXAMPLES:: 33 59 34 EXAMPLES: 35 One shouldn't directly call this constructor since this class is an abstract 36 base class. 37 sage: sage.stats.hmm.hmm.HiddenMarkovModel([[0.4,0.6],[1,0]], [[1,0],[0.5,0.5]], [0.5,0.5]) 38 <sage.stats.hmm.hmm.HiddenMarkovModel object at ...> 39 """ 40 def __init__(self, A, B, pi=None): 60 sage: m = hmm.DiscreteHiddenMarkovModel([[0.4,0.6],[0.1,0.9]], [[0.1,0.9],[0.5,0.5]], [.2,.8]) 61 sage: pi = m.initial_probabilities(); pi 62 [0.2000, 0.8000] 63 sage: type(pi) 64 <type 'sage.finance.time_series.TimeSeries'> 65 66 The returned time series is a copy, so changing it does not 67 change the model. 68 69 sage: pi[0] = .1; pi[1] = .9 70 sage: m.initial_probabilities() 71 [0.2000, 0.8000] 72 73 Some other models:: 74 75 sage: hmm.GaussianHiddenMarkovModel([[.1,.9],[.5,.5]], [(1,1), (-1,1)], [.1,.9]).initial_probabilities() 76 [0.1000, 0.9000] 77 sage: hmm.GaussianMixtureHiddenMarkovModel([[.9,.1],[.4,.6]], [[(.4,(0,1)), (.6,(1,0.1))],[(1,(0,1))]], [.7,.3]).initial_probabilities() 78 [0.7000, 0.3000] 41 79 """ 42 INPUT: 43 A -- matrix or list 44 B -- matrix or list 45 pi -- list of floats 80 return TimeSeries(self.pi) 46 81 47 EXAMPLES: 48 sage: hmm.DiscreteHiddenMarkovModel([[1]], [[0.0,1.0]], [1]) 49 Discrete Hidden Markov Model with 1 States and 2 Emissions 50 Transition matrix: 51 [1.0] 52 Emission matrix: 53 [0.0 1.0] 54 Initial probabilities: [1.0] 82 def transition_matrix(self): 55 83 """ 56 cdef Py_ssize_t n 57 58 # Convert A to a matrix 59 if not is_Matrix(A): 60 n = len(A) 61 A = matrix(RDF, n, len(A[0]) if n > 0 else 0, A) 84 Return the state transition matrix. 62 85 63 # Convert B to a matrix 64 if not is_Matrix(B): 65 n = len(B) 66 B = matrix(RDF, n, len(B[0]) if n > 0 else 0, B) 86 OUTPUT: 67 87 68 # Some consistency checks 69 if not A.is_square(): 70 print A.parent() 71 raise ValueError, "A must be square" 72 if A.nrows() != B.nrows(): 73 raise ValueError, "number of rows of A and B must be the same" 88 - a Sage matrix with real double precision (RDF) entries. 74 89 75 # Make sure A and B are over RDF. 76 if A.base_ring() != RDF: 77 A = A.change_ring(RDF) 78 if B.base_ring() != RDF: 79 B = B.change_ring(RDF) 90 EXAMPLES:: 80 91 81 # Make sure the initial probabilities are all floats. 82 if pi is None: 83 if A.nrows() == 0: 84 self.pi = [] 85 else: 86 self.pi = [1.0/A.nrows()]*A.nrows() 87 else: 88 self.pi = [float(x) for x in pi] 89 if len(self.pi) != A.nrows(): 90 raise ValueError, "length of pi must equal number of rows of A" 92 sage: M = hmm.DiscreteHiddenMarkovModel([[0.7,0.3],[0.9,0.1]], [[0.5,.5],[.1,.9]], [0.3,0.7]) 93 sage: T = M.transition_matrix(); T 94 [0.7 0.3] 95 [0.9 0.1] 91 96 92 # Record the now validated matrices A and B as attributes. 93 # They get used later as attributes in the constructors for 94 # derived classes. 95 self.A = A 96 self.B = B 97 The returned matrix is mutable, but changing it does not 98 change the transition matrix for the model:: 97 99 98 # Check that all entries of A are nonnegative and raise a 99 # ValueError otherwise. Note that we don't check B since it 100 # has entries that are potentially negative in the continuous 101 # case. But GHMM prints clear warnings when the emission 102 # probabilities are negative, i.e., it does not silently give 103 # wrong results like it does for negative transition 104 # probabilities. 105 cdef Py_ssize_t i, j 106 for i from 0 <= i < self.A._nrows: 107 for j from 0 <= j < self.A._ncols: 108 if self.A.get_unsafe_double(i,j) < 0: 109 raise ValueError, "each transition probability must be nonnegative" 110 100 sage: T[0,0] = .1; T[0,1] = .9 101 sage: M.transition_matrix() 102 [0.7 0.3] 103 [0.9 0.1] 111 104 105 Transition matrices for other types of models:: 112 106 113 cdef class DiscreteHiddenMarkovModel(HiddenMarkovModel): 114 """ 115 Create a discrete hidden Markov model. 107 sage: hmm.GaussianHiddenMarkovModel([[.1,.9],[.5,.5]], [(1,1), (-1,1)], [.5,.5]).transition_matrix() 108 [0.1 0.9] 109 [0.5 0.5] 110 sage: hmm.GaussianMixtureHiddenMarkovModel([[.9,.1],[.4,.6]], [[(.4,(0,1)), (.6,(1,0.1))],[(1,(0,1))]], [.7,.3]).transition_matrix() 111 [0.9 0.1] 112 [0.4 0.6] 113 """ 114 from sage.matrix.constructor import matrix 115 from sage.rings.all import RDF 116 return matrix(RDF, self.N, self.A.list()) 116 117 117 hmm.DiscreteHiddenMarkovModel(A, B, pi=None, emission_symbols=None, name=None, normalize=True) 118 n 119 INPUTS: 120 A -- square matrix of doubles; the state change probabilities 121 B -- matrix of doubles; emission probabilities 122 pi -- list of floats; probabilities for each initial state 123 emission_symbols -- list of B.ncols() symbols (just used for printing) 124 name -- (optional) name of the model 125 normalize -- (optional; default=True) whether or not to normalize 126 the model so the probabilities add to 1 127 128 EXAMPLES: 129 We create a discrete HMM with 2 internal states on an alphabet of size 2. 130 sage: hmm.DiscreteHiddenMarkovModel([[0.2,0.8],[0.5,0.5]], [[1,0],[0,1]], [0,1]) 131 Discrete Hidden Markov Model with 2 States and 2 Emissions 132 Transition matrix: 133 [0.2 0.8] 134 [0.5 0.5] 135 Emission matrix: 136 [1.0 0.0] 137 [0.0 1.0] 138 Initial probabilities: [0.0, 1.0] 139 140 The transition probabilities must be nonnegative: 141 sage: hmm.DiscreteHiddenMarkovModel([[-0.2,0.8],[0.5,0.5]], [[1,0],[0,1]], [0,1]) 142 Traceback (most recent call last): 143 ... 144 ValueError: each transition probability must be nonnegative 145 146 The transition probabilities are by default automatically normalized: 147 sage: a = hmm.DiscreteHiddenMarkovModel([[0.2,0.3],[0.5,0.5]], [[1,0],[0,1]], [0,1]) 148 sage: a.transition_matrix() 149 [0.4 0.6] 150 [0.5 0.5] 151 """ 152 def __init__(self, A, B, pi=None, emission_symbols=None, name=None, normalize=True): 118 def graph(self, eps=1e-3): 153 119 """ 154 EXAMPLES: 155 sage: hmm.DiscreteHiddenMarkovModel([[1]], [[0.0,1.0]], [1]) 156 Discrete Hidden Markov Model with 1 States and 2 Emissions 157 Transition matrix: 158 [1.0] 159 Emission matrix: 160 [0.0 1.0] 161 Initial probabilities: [1.0] 162 """ 163 # memory has not all been setup yet. 164 self.initialized = False 165 166 # This sets self.A, self.B and pi after doing appropriate coercions, etc. 167 HiddenMarkovModel.__init__(self, A, B, pi) 168 169 self.set_emission_symbols(emission_symbols) 170 171 self.m = <ghmm_dmodel*> safe_malloc(sizeof(ghmm_dmodel)) 172 173 self.m.label = to_int_array(range(len(self._emission_symbols))) 174 175 # Set all pointers to NULL 176 self.m.s = NULL; self.m.name = NULL; self.m.silent = NULL 177 self.m.tied_to = NULL; self.m.order = NULL; self.m.background_id = NULL 178 self.m.bp = NULL; self.m.topo_order = NULL; self.m.pow_lookup = NULL; 179 self.m.label_alphabet = NULL; self.m.alphabet = NULL 180 181 # Set number of states and number of outputs 182 self.m.N = self.A.nrows() 183 self.m.M = len(self._emission_symbols) 184 # Set the model type to discrete 185 self.m.model_type = GHMM_kDiscreteHMM 186 187 # Set that no a prior model probabilities are set. 188 self.m.prior = -1 189 190 # Assign model identifier if specified 191 if name is not None: 192 name = str(name) 193 self.m.name = <char*> safe_malloc(len(name)+1) 194 strcpy(self.m.name, name) 195 else: 196 self.m.name = NULL 197 198 # Allocate and initialize states 199 cdef ghmm_dstate* states = <ghmm_dstate*> safe_malloc(sizeof(ghmm_dstate) * self.m.N) 200 cdef ghmm_dstate* state 201 202 silent_states = [] 203 tmp_order = [] 204 205 cdef Py_ssize_t i, j, k 206 207 for i from 0 <= i < self.m.N: 208 v = self.B[i] 209 210 # Get a reference to the i-th state for convenience of the notation below. 211 state = &(states[i]) 212 state.desc = NULL 213 214 # Compute state order 215 if self.m.M > 1: 216 order = math.log( len(v), self.m.M ) - 1 217 else: 218 order = len(v) - 1 219 220 # Check for valid number of emission parameters 221 order = int(order) 222 if self.m.M**(order+1) == len(v): 223 tmp_order.append(order) 224 else: 225 raise ValueError, "number of columns (= %s) of B must be a power of the number of emission symbols (= %s)"%( 226 self.B.ncols(), len(emission_symbols)) 227 228 state.b = to_double_array(v) 229 state.pi = self.pi[i] 230 231 silent_states.append( 1 if sum(v) == 0 else 0) 232 233 # Set "out" probabilities, i.e., the probabilities to 234 # transition to another hidden state from this state. 235 v = self.A[i] 236 k = self.m.N 237 state.out_states = k 238 state.out_id = <int*> safe_malloc(sizeof(int)*k) 239 state.out_a = <double*> safe_malloc(sizeof(double)*k) 240 for j from 0 <= j < k: 241 state.out_id[j] = j 242 state.out_a[j] = v[j] 243 244 # Set "in" probabilities 245 v = self.A.column(i) 246 state.in_states = k 247 state.in_id = <int*> safe_malloc(sizeof(int)*k) 248 state.in_a = <double*> safe_malloc(sizeof(double)*k) 249 for j from 0 <= j < k: 250 state.in_id[j] = j 251 state.in_a[j] = v[j] 252 253 state.fix = 0 254 255 self.m.s = states 256 self.initialized = True 257 if normalize: 258 self.normalize() 259 260 def __cmp__(self, other): 261 """ 262 Compare two Discrete HMM's. 120 Create a weighted directed graph from the transition matrix, 121 not including any edge with a probability less than eps. 263 122 264 123 INPUT: 265 self, other -- objects; if other is not a Discrete HMM compare types.266 OUTPUT:267 -1,0,1268 124 269 The transition matrices are compared, then the emission 270 parameters, and the initial probabilities. The names are not 271 compared, so two GHMM's with the same parameters but different 272 names compare to be equal. 273 274 EXAMPLES: 275 sage: m = hmm.DiscreteHiddenMarkovModel([[0.4,0.6],[0.1,0.9]], [[0.0,1.0],[1,1]], [1,2], normalize=False) 276 sage: m.__cmp__(m) 277 0 278 279 Note that the name doesn't matter: 280 sage: n = hmm.DiscreteHiddenMarkovModel([[0.4,0.6],[0.1,0.9]], [[0.0,1.0],[1,1]], [1,2], normalize=False) 281 sage: m.__cmp__(n) 282 0 283 284 Normalizing fixes the initial probabilities, hence m and n are no longer equal. 285 sage: n.normalize() 286 sage: m.__cmp__(n) 287 1 288 """ 289 if not isinstance(other, DiscreteHiddenMarkovModel): 290 return cmp(type(self), type(other)) 291 292 if self is other: return 0 # easy special case 293 294 cdef DiscreteHiddenMarkovModel o = other 295 if self.m.N < o.m.N: 296 return -1 297 elif self.m.N > o.m.N: 298 return 1 299 cdef Py_ssize_t i, j 300 301 # This code is similar to code in chmm.pyx, but with several small differences. 302 303 # Compare transition matrices 304 for i from 0 <= i < self.m.N: 305 for j from 0 <= j < self.m.s[i].out_states: 306 if self.m.s[i].out_a[j] < o.m.s[i].out_a[j]: 307 return -1 308 elif self.m.s[i].out_a[j] > o.m.s[i].out_a[j]: 309 return 1 310 311 # Compare emission matrices 312 for i from 0 <= i < self.m.N: 313 for j from 0 <= j < self.B._ncols: 314 if self.m.s[i].b[j] < o.m.s[i].b[j]: 315 return -1 316 elif self.m.s[i].b[j] > o.m.s[i].b[j]: 317 return 1 318 319 # Compare initial state probabilities 320 for 0 <= i < self.m.N: 321 if self.m.s[i].pi < o.m.s[i].pi: 322 return -1 323 elif self.m.s[i].pi > o.m.s[i].pi: 324 return 1 325 326 # Compare emission symbols 327 return cmp(self._emission_symbols, o._emission_symbols) 328 329 def __reduce__(self): 330 """ 331 Used in pickling. 332 333 EXAMPLES: 334 sage: m = hmm.DiscreteHiddenMarkovModel([[0.4,0.6],[0.1,0.9]], [[0.0,1.0],[1,1]], [0,1], ['a','b'], name='test model') 335 sage: f,g = m.__reduce__() 336 sage: f(*g) == m 337 True 338 """ 339 return unpickle_discrete_hmm_v0, (self.transition_matrix(), self.emission_matrix(), 340 self.initial_probabilities(), self._emission_symbols, self.name()) 341 342 def __dealloc__(self): 343 """ 344 Deallocate memory allocated by the HMM. 345 346 EXAMPLES: 347 sage: m = hmm.DiscreteHiddenMarkovModel([[0.2,0.8],[0.5,0.5]], [[1,0],[0,1]], [0,1]) # indirect doctest 348 sage: del m 349 """ 350 if self.initialized: 351 ghmm_dmodel_free(&self.m) 352 353 def fix_emissions(self, Py_ssize_t i, bint fixed=True): 354 """ 355 Sets the i-th emission parameters to be either fixed or not 356 fixed. If it is fixed, then running the Baum-Welch algorithm 357 will not change the emission parameters for the i-th state. 358 359 INPUT: 360 i -- nonnegative integer < self.m.N 361 fixed -- bool 362 363 EXAMPLES: 364 First without calling fix_emissions: 365 sage: m = hmm.DiscreteHiddenMarkovModel([[0.4,0.6],[0.1,0.9]],[[.5,.5],[.5,.5]]) 366 sage: m.baum_welch([0,0,0,1,1,1]) 367 sage: m.emission_matrix() 368 [ 1.0 0.0] 369 [3.92881039079e-05 0.999960711896] 370 371 We call fix_emissions on the first state and notice that the first 372 row of the emission matrix does not change: 373 sage: m = hmm.DiscreteHiddenMarkovModel([[0.4,0.6],[0.1,0.9]],[[.5,.5],[.5,.5]]) 374 sage: m.fix_emissions(0) 375 sage: m.baum_welch([0,0,0,1,1,1]) 376 sage: m.emission_matrix() 377 [ 0.5 0.5] 378 [0.000542712675606 0.999457287324] 379 380 We call fix_emissions on the second state and notice that the second 381 row of the emission matrix does not change: 382 sage: m = hmm.DiscreteHiddenMarkovModel([[0.4,0.6],[0.1,0.9]],[[.5,.5],[.5,.5]]) 383 sage: m.fix_emissions(1) 384 sage: m.baum_welch([0,0,0,1,1,1]) 385 sage: m.emission_matrix() 386 [ 0.999999904763 9.52366620142e-08] 387 [ 0.5 0.5] 388 389 TESTS: 390 Make sure that out of range indices are handled correctly with an IndexError. 391 sage: m = hmm.DiscreteHiddenMarkovModel([[0.4,0.6],[0.1,0.9]],[[.5,.5],[.5,.5]]) 392 sage: m.fix_emissions(2) 393 Traceback (most recent call last): 394 ... 395 IndexError: index out of range 396 sage: m.fix_emissions(-1) 397 Traceback (most recent call last): 398 ... 399 IndexError: index out of range 400 """ 401 if i < 0 or i >= self.m.N: 402 raise IndexError, "index out of range" 403 self.m.s[i].fix = fixed 404 405 def __repr__(self): 406 """ 407 Return string representation of this HMM. 125 - eps -- nonnegative real number 408 126 409 127 OUTPUT: 410 string411 128 412 EXAMPLES: 413 sage: a = hmm.DiscreteHiddenMarkovModel([[0.1,0.9],[0.1,0.9]], [[0.9,0.1],[0.1,0.9]], [0.5,0.5], [3/4, 'abc']) 414 sage: a.__repr__() 415 "Discrete Hidden Markov Model with 2 States and 2 Emissions\nTransition matrix:\n[0.1 0.9]\n[0.1 0.9]\nEmission matrix:\n[0.9 0.1]\n[0.1 0.9]\nInitial probabilities: [0.5, 0.5]\nEmission symbols: [3/4, 'abc']" 129 - a digraph 130 131 EXAMPLES:: 132 133 sage: m = hmm.DiscreteHiddenMarkovModel([[.3,0,.7],[0,0,1],[.5,.5,0]], [[.5,.5,.2]]*3, [1/3]*3) 134 sage: G = m.graph(); G 135 Looped multi-digraph on 3 vertices 136 sage: G.edges() 137 [(0, 0, 0.3), (0, 2, 0.7), (1, 2, 1.0), (2, 0, 0.5), (2, 1, 0.5)] 138 sage: G.plot() 416 139 """ 417 s = "Discrete Hidden Markov Model%s with %s States and %s Emissions"%( 418 ' ' + self.m.name if self.m.name else '', 419 self.m.N, self.m.M) 420 s += '\nTransition matrix:\n%s'%self.transition_matrix() 421 s += '\nEmission matrix:\n%s'%self.emission_matrix() 422 s += '\nInitial probabilities: %s'%self.initial_probabilities() 423 if self._emission_symbols_dict: 424 s += '\nEmission symbols: %s'%self._emission_symbols 425 return s 140 cdef int i, j 141 m = self.transition_matrix() 142 for i in range(self.N): 143 for j in range(self.N): 144 if m[i,j] < eps: 145 m[i,j] = 0 146 from sage.graphs.all import DiGraph 147 return DiGraph(m, weighted=True) 426 148 427 def name(self): 428 """ 429 Return the name of this model. 430 431 OUTPUT: 432 string or None 433 434 EXAMPLES: 435 sage: m = hmm.DiscreteHiddenMarkovModel([[0.4,0.6],[0.1,0.9]], [[0.0,1.0],[1,1]], [1,2], name='test model') 436 sage: m.name() 437 'test model' 438 439 If the model is not explicitly named then this function returns None: 440 sage: m = hmm.DiscreteHiddenMarkovModel([[0.4,0.6],[0.1,0.9]], [[0.0,1.0],[1,1]], [1,2]) 441 sage: m.name() is None 442 True 443 """ 444 if self.m.name: 445 s = str(self.m.name) 446 return s 447 else: 448 return None 449 450 def initial_probabilities(self): 451 """ 452 Return the list of initial state probabilities. 453 454 OUTPUT: 455 list of floats 456 457 EXAMPLES: 458 sage: a = hmm.DiscreteHiddenMarkovModel([[0.9,0.1],[0.9,0.1]], [[0.5,0.5,0,0],[0,0,.5,.5]], [0.5,0.5], [1,-1,3,-3]) 459 sage: a.initial_probabilities() 460 [0.5, 0.5] 461 """ 462 cdef Py_ssize_t i 463 return [self.m.s[i].pi for i in range(self.m.N)] 464 465 def transition_matrix(self, list_only=True): 466 """ 467 Return the hidden state transition matrix. 468 469 EXAMPLES: 470 sage: a = hmm.DiscreteHiddenMarkovModel([[0.9,0.1],[0.9,0.1]], [[0.5,0.5,0,0],[0,0,.5,.5]], [0.5,0.5], [1,-1,3,-3]) 471 sage: a.transition_matrix() 472 [0.9 0.1] 473 [0.9 0.1] 474 """ 475 cdef Py_ssize_t i, j 476 for i from 0 <= i < self.m.N: 477 for j from 0 <= j < self.m.s[i].out_states: 478 self.A.set_unsafe_double(i,j,self.m.s[i].out_a[j]) 479 return self.A 480 481 def emission_matrix(self, list_only=True): 482 """ 483 Return the emission probability matrix. 484 485 EXAMPLES: 486 sage: a = hmm.DiscreteHiddenMarkovModel([[0.9,0.1],[0.9,0.1]], [[0.5,0.5,0,0],[0,0,.5,.5]], [0.5,0.5], [1,-1,3,-3]) 487 sage: a.emission_matrix() 488 [0.5 0.5 0.0 0.0] 489 [0.0 0.0 0.5 0.5] 490 """ 491 cdef Py_ssize_t i, j 492 for i from 0 <= i < self.m.N: 493 for j from 0 <= j < self.B._ncols: 494 self.B.set_unsafe_double(i,j,self.m.s[i].b[j]) 495 return self.B 496 497 def normalize(self): 498 """ 499 Normalize the transition and emission probabilities, if applicable. 500 501 EXAMPLES: 502 sage: a = hmm.DiscreteHiddenMarkovModel([[0.5,1],[1.2,0.9]], [[1,0.5],[0.5,1]], [0.1,1.2]) 503 sage: a.normalize() 504 sage: a 505 Discrete Hidden Markov Model with 2 States and 2 Emissions 506 Transition matrix: 507 [0.333333333333 0.666666666667] 508 [0.571428571429 0.428571428571] 509 Emission matrix: 510 [0.666666666667 0.333333333333] 511 [0.333333333333 0.666666666667] 512 Initial probabilities: [0.076923076923076927, 0.92307692307692302] 513 """ 514 ghmm_dmodel_normalize(self.m) 515 516 def sample(self, long length, number=None): 149 def sample(self, Py_ssize_t length, number=None): 517 150 """ 518 151 Return number samples from this HMM of given length. 519 152 … … 529 162 sage: set_random_seed(0) 530 163 sage: a = hmm.DiscreteHiddenMarkovModel([[0.1,0.9],[0.1,0.9]], [[1,0],[0,1]], [0,1]) 531 164 sage: print a.sample(10, 3) 532 [[1, 0, 1, 1, 0, 1, 1, 0, 1, 0], [1, 1, 1, 1, 1, 1, 1, 1, 1, 0], [1, 1, 1, 1, 1, 1, 1, 1, 1, 1]]165 [[1, 0, 1, 1, 1, 1, 0, 1, 1, 1], [1, 1, 0, 0, 1, 1, 1, 1, 1, 1], [1, 1, 1, 1, 0, 1, 0, 1, 1, 1]] 533 166 sage: a.sample(15) 534 [1, 1, 1, 1, 1, 1, 0, 1, 1, 0, 1, 0, 1, 1, 1] 167 [1, 1, 1, 1, 0 ... 1, 1, 1, 1, 1] 168 sage: a.sample(3, 1) 169 [[1, 1, 1]] 535 170 sage: list(a.sample(1000)).count(0) 536 95171 88 537 172 538 173 If the emission symbols are set 539 174 sage: set_random_seed(0) 540 175 sage: a = hmm.DiscreteHiddenMarkovModel([[0.5,0.5],[0.1,0.9]], [[1,0],[0,1]], [0,1], ['up', 'down']) 541 176 sage: a.sample(10) 542 ['down', 'up', 'down', 'down', ' up', 'down', 'down', 'up', 'down', 'up']177 ['down', 'up', 'down', 'down', 'down', 'down', 'up', 'up', 'up', 'up'] 543 178 """ 544 seed = random()545 179 if number is None: 546 number = 1 547 single = True 548 else: 549 single = False 550 cdef ghmm_dseq *d = ghmm_dmodel_generate_sequences(self.m, seed, length, number, length) 551 cdef Py_ssize_t i, j 552 v = [[d.seq[j][i] for i in range(length)] for j in range(number)] 553 if self._emission_symbols_dict: 554 w = self._emission_symbols 555 v = [[w[i] for i in z] for z in v] 556 if single: 557 return v[0] 558 return v 180 return self.generate_sequence(length)[0] 559 181 560 def emission_symbols(self): 561 """ 562 Return a copy of the list of emission symbols of self. 182 cdef Py_ssize_t i 183 return [self.generate_sequence(length)[0] for i in range(number)] 563 184 564 Use set_emission_symbols to set the list of emission symbols. 185 cdef class DiscreteHiddenMarkovModel(HiddenMarkovModel): 186 """ 187 A discrete Hidden Markov model implemented using double precision 188 floating point arithmetic. 565 189 566 EXAMPLES: 567 sage: a = hmm.DiscreteHiddenMarkovModel([[0.5,0.5],[0.1,0.9]], [[1,0],[0,1]], [0,1], ['up', -3/179]) 568 sage: a.emission_symbols() 569 ['up', -3/179] 570 """ 571 return list(self._emission_symbols) 190 EXAMPLES:: 572 191 573 def set_emission_symbols(self, emission_symbols): 574 """ 575 Set the list of emission symbols. 576 577 EXAMPLES: 578 sage: set_random_seed(0) 579 sage: a = hmm.DiscreteHiddenMarkovModel([[0.5,0.5],[0.1,0.9]], [[1,0],[0,1]], [0,1], ['up', 'down']) 580 sage: a.set_emission_symbols([3,5]) 581 sage: a.emission_symbols() 582 [3, 5] 583 sage: a.sample(10) 584 [5, 3, 5, 5, 3, 5, 5, 3, 5, 3] 585 sage: a.set_emission_symbols([pi,5/9+e]) 586 sage: a.sample(10) 587 [e + 5/9, e + 5/9, e + 5/9, e + 5/9, e + 5/9, e + 5/9, pi, pi, e + 5/9, pi] 588 """ 589 if emission_symbols is None: 590 self._emission_symbols = range(self.B.ncols()) 591 self._emission_symbols_dict = None 592 else: 593 self._emission_symbols = list(emission_symbols) 594 if self._emission_symbols != range(self.B.ncols()): 595 self._emission_symbols_dict = dict([(x,i) for i, x in enumerate(emission_symbols)]) 596 597 598 #################################################################### 599 # HMM Problem 1 -- Computing likelihood: Given the parameter set 600 # lambda of an HMM model and an observation sequence O, determine 601 # the likelihood P(O | lambda). 602 #################################################################### 603 def log_likelihood(self, seq): 604 r""" 605 HMM Problem 1: Likelihood. Return $\log ( P[seq | model] )$, 606 the log of the probability of seeing the given sequence given 607 this model, using the forward algorithm and assuming 608 independence of the sequence seq. 609 610 INPUT: 611 seq -- a list; sequence of observed emissions of the HMM 612 613 OUTPUT: 614 float -- the log of the probability of seeing this sequence 615 given the model 616 617 WARNING: By convention we return -inf for 0 probability events. 618 619 EXAMPLES: 620 sage: a = hmm.DiscreteHiddenMarkovModel([[0.1,0.9],[0.1,0.9]], [[1,0],[0,1]], [0,1]) 621 sage: a.log_likelihood([1,1]) 622 -0.10536051565782635 623 sage: a.log_likelihood([1,0]) 624 -2.3025850929940459 625 626 Notice that any sequence starting with 0 can't occur, since 627 the above model always starts in a state that produces 1 with 628 probability 1. By convention log(probability 0) is -inf. 629 sage: a.log_likelihood([0,0]) 630 -inf 631 632 Here's a special case where each sequence is equally probable. 633 sage: a = hmm.DiscreteHiddenMarkovModel([[0.5,0.5],[0.5,0.5]], [[1,0],[0,1]], [0.5,0.5]) 634 sage: a.log_likelihood([0,0]) 635 -1.3862943611198906 636 sage: log(0.25) 637 -1.38629436111989 638 sage: a.log_likelihood([0,1]) 639 -1.3862943611198906 640 sage: a.log_likelihood([1,0]) 641 -1.3862943611198906 642 sage: a.log_likelihood([1,1]) 643 -1.3862943611198906 644 """ 645 if self._emission_symbols_dict: 646 seq = [self._emission_symbols_dict[z] for z in seq] 647 cdef double log_p 648 cdef int* O = to_int_array(seq) 649 cdef int ret = ghmm_dmodel_logp(self.m, O, len(seq), &log_p) 650 sage_free(O) 651 if ret == -1: 652 # forward returned -1: sequence can't be built 653 return -float('Inf') 654 return log_p 655 656 #################################################################### 657 # HMM Problem 2 -- Decoding: Given the complete parameter set that 658 # defines the model and an observation sequence seq, determine the 659 # best hidden sequence Q. Use the Viterbi algorithm. 660 #################################################################### 661 def viterbi(self, seq): 662 """ 663 HMM Problem 2: Decoding. Determine a hidden sequence of 664 states that is most likely to produce the given sequence seq 665 of observations. 666 667 INPUT: 668 seq -- sequence of emitted symbols 669 670 OUTPUT: 671 list -- a most probable sequence of hidden states, i.e., the 672 Viterbi path. 673 float -- log of the probability that the sequence of hidden 674 states actually produced the given sequence seq. 675 676 EXAMPLES: 677 sage: a = hmm.DiscreteHiddenMarkovModel([[0.1,0.9],[0.1,0.9]], [[0.9,0.1],[0.1,0.9]], [0.5,0.5]) 678 sage: a.viterbi([1,0,0,1,0,0,1,1]) 679 ([1, 0, 0, 1, ..., 0, 1, 1], -11.06245322477221...) 680 681 We predict the state sequence when the emissions are 3/4 and 'abc'. 682 sage: a = hmm.DiscreteHiddenMarkovModel([[0.1,0.9],[0.1,0.9]], [[0.9,0.1],[0.1,0.9]], [0.5,0.5], [3/4, 'abc']) 683 684 Note that state 0 is common below, despite the model trying hard to 685 switch to state 1: 686 sage: a.viterbi([3/4, 'abc', 'abc'] + [3/4]*10) 687 ([0, 1, 1, 0, 0, ..., 0, 0, 0, 0, 0, 0, 0], -25.299405845367794) 688 """ 689 if len(seq) == 0: 690 return [], 0.0 691 if self._emission_symbols_dict: 692 seq = [self._emission_symbols_dict[z] for z in seq] 693 cdef int* path 694 cdef int* O = to_int_array(seq) 695 cdef int pathlen 696 cdef double log_p 697 698 path = ghmm_dmodel_viterbi(self.m, O, len(seq), &pathlen, &log_p) 699 sage_free(O) 700 if not path: 701 raise RuntimeError, "error computing viterbi path" 702 p = [path[i] for i in range(pathlen)] 703 sage_free(path) 704 705 return p, log_p 706 707 #################################################################### 708 # HMM Problem 3 -- Learning: Given an observation sequence O and 709 # the set of states in the HMM, improve the HMM to increase the 710 # probability of observing O. 711 #################################################################### 712 def baum_welch(self, training_seqs, nsteps=None, log_likelihood_cutoff=None): 713 """ 714 HMM Problem 3: Learning. Given an observation sequence O and 715 the set of states in the HMM, improve the HMM using the 716 Baum-Welch algorithm to increase the probability of observing O. 717 718 INPUT: 719 training_seqs -- a list of lists of emission symbols (or a single list) 720 nsteps -- integer or None (default: None) maximum number 721 of Baum-Welch steps to take 722 log_likehood_cutoff -- positive float or None (default: 723 None); the minimal improvement in likelihood 724 with respect to the last iteration required to 725 continue. Relative value to log likelihood 726 727 OUTPUT: 728 changes the model in places, or raises a RuntimError 729 exception on error 730 731 EXAMPLES: 732 We make a model that has two states and is equally likely to output 733 0 or 1 in either state and transitions back and forth with equal 734 probability. 735 sage: a = hmm.DiscreteHiddenMarkovModel([[0.5,0.5],[0.5,0.5]], [[0.5,0.5],[0.5,0.5]], [0.5,0.5]) 736 737 We give the model some training data this much more likely to 738 be 1 than 0. 739 sage: a.baum_welch([[1,1,1,1,0,1], [1,0,1,1,1,1]]) 740 741 Now the model's emission matrix changes since it is much 742 more likely to emit 1 than 0. 743 sage: a 744 Discrete Hidden Markov Model with 2 States and 2 Emissions 745 Transition matrix: 746 [0.5 0.5] 747 [0.5 0.5] 748 Emission matrix: 749 [0.166666666667 0.833333333333] 750 [0.166666666667 0.833333333333] 751 Initial probabilities: [0.5, 0.5] 752 753 Note that 1/6 = 1.666...: 754 sage: 1.0/6 755 0.166666666666667 756 757 We run Baum-Welch on a single sequence: 758 sage: a = hmm.DiscreteHiddenMarkovModel([[0.5,0.5],[0.5,0.5]], [[0.5,0.5],[0.5,0.5]], [0.5,0.5]) 759 sage: a.baum_welch([1,0,1]*10) 760 sage: a 761 Discrete Hidden Markov Model with 2 States and 2 Emissions 762 Transition matrix: 763 [0.5 0.5] 764 [0.5 0.5] 765 Emission matrix: 766 [0.333333333333 0.666666666667] 767 [0.333333333333 0.666666666667] 768 Initial probabilities: [0.5, 0.5] 769 770 We compare using a non-default number of steps and non-default log likelihood cutoff: 771 sage: h = hmm.DiscreteHiddenMarkovModel([[.1,.9],[.4,.6]], [[1/2]*2]*2, [1/2]*2) 772 sage: h.baum_welch([1,1,1,1,1,0,1,1,1,1,1,0,0]) 773 sage: h 774 Discrete Hidden Markov Model with 2 States and 2 Emissions 775 Transition matrix: 776 [ 0.643888431046 0.356111568954] 777 [0.00232031442167 0.997679685578] 778 Emission matrix: 779 [ 0.0 1.0] 780 [0.296080269308 0.703919730692] 781 Initial probabilities: [1.0, 3.662034729231...e-20] 782 sage: h = hmm.DiscreteHiddenMarkovModel([[.1,.9],[.4,.6]], [[1/2]*2]*2, [1/2]*2) 783 sage: h.baum_welch([1,1,1,1,1,0,1,1,1,1,1,0,0], nsteps=1) 784 sage: h 785 Discrete Hidden Markov Model with 2 States and 2 Emissions 786 Transition matrix: 787 [0.1 0.9] 788 [0.4 0.6] 789 Emission matrix: 790 [0.222426510432 0.777573489568] 791 [0.234678486789 0.765321513211] 792 Initial probabilities: [0.5, 0.5] 793 sage: h = hmm.DiscreteHiddenMarkovModel([[.1,.9],[.4,.6]], [[1/2]*2]*2, [1/2]*2) 794 sage: h.baum_welch([1,1,1,1,1,0,1,1,1,1,1,0,0], log_likelihood_cutoff=0.01) 795 sage: h 796 Discrete Hidden Markov Model with 2 States and 2 Emissions 797 Transition matrix: 798 [0.0999471960754 0.900052803925] 799 [ 0.399248483887 0.600751516113] 800 Emission matrix: 801 [0.209328925617 0.790671074383] 802 [ 0.24081056723 0.75918943277] 803 Initial probabilities: [0.508779821929586..., 0.491220178070413...] 804 805 TESTS: 806 We test training with non-default string symbols: 807 sage: a = hmm.DiscreteHiddenMarkovModel([[0.5,0.5],[0.5,0.5]], [[0.5,0.5],[0.5,0.5]], [0.5,0.5], ['up','down']) 808 sage: a.baum_welch([['up','up'], ['down','up']]) 809 sage: a 810 Discrete Hidden Markov Model with 2 States and 2 Emissions 811 Transition matrix: 812 [0.5 0.5] 813 [0.5 0.5] 814 Emission matrix: 815 [0.75 0.25] 816 [0.75 0.25] 817 Initial probabilities: [0.5, 0.5] 818 Emission symbols: ['up', 'down'] 819 820 NOTE: Training for models including silent states is not yet supported. 821 822 REFERENCES: 823 Rabiner, L.R.: "`A Tutorial on Hidden Markov Models and Selected 824 Applications in Speech Recognition"', Proceedings of the IEEE, 825 77, no 2, 1989, pp 257--285. 826 """ 827 if len(training_seqs) > 0 and not isinstance(training_seqs[0], (list, tuple)): 828 training_seqs = [training_seqs] 829 830 if self._emission_symbols_dict: 831 seqs = [[self._emission_symbols_dict[z] for z in x] for x in training_seqs] 832 else: 833 seqs = training_seqs 834 835 cdef ghmm_dseq* d = malloc_ghmm_dseq(seqs) 836 837 if nsteps or log_likelihood_cutoff: 838 if nsteps is None: 839 nsteps = GHMM_MAX_ITER_BW 840 if log_likelihood_cutoff is None: 841 log_likelihood_cutoff = GHMM_EPS_ITER_BW 842 if ghmm_dmodel_baum_welch_nstep(self.m, d, nsteps, log_likelihood_cutoff): 843 raise RuntimeError, "error running Baum-Welch algorithm" 844 else: 845 if ghmm_dmodel_baum_welch(self.m, d): 846 raise RuntimeError, "error running Baum-Welch algorithm" 847 848 ghmm_dseq_free(&d) 849 850 851 ################################################################################## 852 # Helper Functions 853 ################################################################################## 854 855 cdef ghmm_dseq* malloc_ghmm_dseq(seqs) except NULL: 856 """ 857 Allocate a discrete sequence of samples. 858 859 INPUT: 860 seqs -- a list of sequences 861 862 OUTPUT: 863 C pointer to ghmm_dseq 864 """ 865 cdef ghmm_dseq* d = ghmm_dseq_calloc(len(seqs)) 866 if d == NULL: 867 raise MemoryError 868 cdef int i, j, m, n 869 m = len(seqs) 870 d.seq_number = m 871 d.capacity = m 872 d.total_w = m 873 for i from 0 <= i < m: 874 v = seqs[i] 875 n = len(v) 876 d.seq[i] = <int*> safe_malloc(sizeof(int) * n) 877 for j from 0 <= j < n: 878 d.seq[i][j] = v[j] 879 d.seq_len[i] = n 880 d.seq_id[i] = i 881 d.seq_w[i] = 1 882 d.flags = 0 883 return d 884 885 886 def unpickle_discrete_hmm_v0(A, B, pi, emission_symbols,name): 887 """ 888 TESTS: 889 sage: m = hmm.DiscreteHiddenMarkovModel([[0.4,0.6],[0.1,0.9]], [[0.0,1.0],[0.5,0.5]], [1,0], name='test model') 890 sage: loads(dumps(m)) == m 891 True 892 sage: loads(dumps(m)).name() 893 'test model' 894 sage: sage.stats.hmm.hmm.unpickle_discrete_hmm_v0(m.transition_matrix(), m.emission_matrix(), m.initial_probabilities(), ['a','b'], m.name()) 895 Discrete Hidden Markov Model test model with 2 States and 2 Emissions 192 sage: m = hmm.DiscreteHiddenMarkovModel([[0.4,0.6],[0.1,0.9]], [[0.1,0.9],[0.5,0.5]], [.5,.5]); m 193 Discrete Hidden Markov Model with 2 States and 2 Emissions 896 194 Transition matrix: 897 195 [0.4 0.6] 898 196 [0.1 0.9] 899 197 Emission matrix: 900 [0. 0 1.0]198 [0.1 0.9] 901 199 [0.5 0.5] 902 Initial probabilities: [1.0, 0.0] 903 Emission symbols: ['a', 'b'] 200 Initial probabilities: [0.5000, 0.5000] 201 sage: m.log_likelihood([0,1,0,1,0,1]) 202 -4.66693474691329... 203 sage: m.viterbi([0,1,0,1,0,1]) 204 ([1, 1, 1, 1, 1, 1], -5.3788328422087481) 205 sage: m.baum_welch([0,1,0,1,0,1]) 206 (0.0, 22) 207 sage: m 208 Discrete Hidden Markov Model with 2 States and 2 Emissions 209 Transition matrix: 210 [1.0134345614...e-70 1.0] 211 [ 1.0 3.99743527136e-19] 212 Emission matrix: 213 [7.3802215662...e-54 1.0] 214 [ 1.0 3.9974352626e-19] 215 Initial probabilities: [0.0000, 1.0000] 216 sage: m.sample(10) 217 [0, 1, 0, 1, 0, 1, 0, 1, 0, 1] 218 sage: m.graph().plot() 904 219 """ 905 return DiscreteHiddenMarkovModel(A,B,pi,emission_symbols,name) 220 cdef TimeSeries B 221 cdef int n_out 222 cdef object _emission_symbols, _emission_symbols_dict 223 224 def __init__(self, A, B, pi, emission_symbols=None, bint normalize=True): 225 """ 226 Create a discrete emissions HMM with transition probability 227 matrix A, emission probabilities given by B, initial state 228 probabilities pi, and given emission symbols (which default 229 to the first few nonnegative integers). 230 231 INPUT:: 232 233 - A -- a list of lists or a square N x N matrix, whose 234 (i,j) entry gives the probability of transitioning from 235 state i to state j. 236 237 - B -- a list of N lists or a matrix with N rows, such that 238 B[i,k] gives the probability of emitting symbol k while 239 in state i. 240 241 - pi -- the probabilities of starting in each initial 242 state, i.e,. pi[i] is the probability of starting in 243 state i. 244 245 - emission_symbols -- None or list (default: None); if 246 None, the emission_symbols are the ints [0..N-1], where N 247 is the number of states. Otherwise, they are the entries 248 of the list emissions_symbols, which must all be hashable. 249 250 - normalize --bool (default: True); if given, input is 251 normalized to define valid probability distributions, 252 e.g., the entries of A are made nonnegative and the rows 253 sum to 1. 254 255 EXAMPLES:: 256 257 sage: hmm.DiscreteHiddenMarkovModel([.5,0,-1,.5], [[1],[1]],[.5,.5]).transition_matrix() 258 [1.0 0.0] 259 [0.0 1.0] 260 sage: hmm.DiscreteHiddenMarkovModel([0,0,.1,.9], [[1],[1]],[.5,.5]).transition_matrix() 261 [0.5 0.5] 262 [0.1 0.9] 263 sage: hmm.DiscreteHiddenMarkovModel([-1,-2,.1,.9], [[1],[1]],[.5,.5]).transition_matrix() 264 [0.5 0.5] 265 [0.1 0.9] 266 sage: hmm.DiscreteHiddenMarkovModel([1,2,.1,1.2], [[1],[1]],[.5,.5]).transition_matrix() 267 [ 0.333333333333 0.666666666667] 268 [0.0769230769231 0.923076923077] 269 """ 270 self.pi = util.initial_probs_to_TimeSeries(pi, normalize) 271 self.N = len(self.pi) 272 self.A = util.state_matrix_to_TimeSeries(A, self.N, normalize) 273 self._emission_symbols = emission_symbols 274 if self._emission_symbols is not None: 275 self._emission_symbols_dict = dict([(y,x) for x,y in enumerate(emission_symbols)]) 276 277 # TODO : normalization 278 if not is_Matrix(B): 279 B = matrix(B) 280 if B.nrows() != self.N: 281 raise ValueError, "number of rows of B must equal number of states" 282 self.B = TimeSeries(B.list()) 283 self.n_out = B.ncols() 284 if emission_symbols is not None and len(emission_symbols) != self.n_out: 285 raise ValueError, "number of emission symbols must equal number of output states" 286 cdef Py_ssize_t i 287 if normalize: 288 for i in range(self.N): 289 util.normalize_probability_TimeSeries(self.B, i*self.n_out, (i+1)*self.n_out) 290 291 def __reduce__(self): 292 """ 293 Used in pickling. 294 295 EXAMPLES: 296 sage: m = hmm.DiscreteHiddenMarkovModel([[0.4,0.6],[0.1,0.9]], [[0.0,1.0],[1,1]], [0,1], ['a','b']) 297 sage: loads(dumps(m)) == m 298 True 299 """ 300 return unpickle_discrete_hmm_v1, \ 301 (self.A, self.B, self.pi, self.n_out, self._emission_symbols, self._emission_symbols_dict) 302 303 def __cmp__(self, other): 304 """ 305 EXAMPLES:: 306 307 sage: m = hmm.DiscreteHiddenMarkovModel([[0.4,0.6],[0.1,0.9]], [[0.0,1.0],[0.5,0.5]], [.5,.5]) 308 sage: n = hmm.DiscreteHiddenMarkovModel([[0.4,0.6],[0.1,0.9]], [[0.0,1.0],[0.5,0.5]], [1,0]) 309 sage: m == n 310 False 311 sage: m == m 312 True 313 sage: m < n 314 True 315 sage: n < m 316 False 317 """ 318 if not isinstance(other, DiscreteHiddenMarkovModel): 319 raise ValueError 320 return cmp(self.__reduce__()[1], other.__reduce__()[1]) 321 322 323 def emission_matrix(self): 324 """ 325 Return the matrix whose i-th row specifies the emission 326 probability distribution for the i-th state. More precisely, 327 the i,j entry of the matrix is the probability of the Markov 328 model outputing the j-th symbol when it is in the i-th state. 329 330 OUTPUT: 331 332 - a Sage matrix with real double precision (RDF) entries. 333 334 EXAMPLES:: 335 336 sage: m = hmm.DiscreteHiddenMarkovModel([[0.4,0.6],[0.1,0.9]], [[0.1,0.9],[0.5,0.5]], [.5,.5]) 337 sage: E = m.emission_matrix(); E 338 [0.1 0.9] 339 [0.5 0.5] 340 341 The returned matrix is mutable, but changing it does not 342 change the transition matrix for the model:: 343 344 sage: E[0,0] = 0; E[0,1] = 1 345 sage: m.emission_matrix() 346 [0.1 0.9] 347 [0.5 0.5] 348 """ 349 from sage.matrix.constructor import matrix 350 from sage.rings.all import RDF 351 return matrix(RDF, self.N, self.n_out, self.B.list()) 352 353 354 def __repr__(self): 355 """ 356 Return string representation of this discrete hidden Markov model. 357 358 EXAMPLES:: 359 360 sage: m = hmm.DiscreteHiddenMarkovModel([[0.4,0.6],[0.1,0.9]], [[0.1,0.9],[0.5,0.5]], [.2,.8]) 361 sage: m.__repr__() 362 'Discrete Hidden Markov Model with 2 States and 2 Emissions\nTransition matrix:\n[0.4 0.6]\n[0.1 0.9]\nEmission matrix:\n[0.1 0.9]\n[0.5 0.5]\nInitial probabilities: [0.2000, 0.8000]' 363 """ 364 s = "Discrete Hidden Markov Model with %s States and %s Emissions"%( 365 self.N, self.n_out) 366 s += '\nTransition matrix:\n%s'%self.transition_matrix() 367 s += '\nEmission matrix:\n%s'%self.emission_matrix() 368 s += '\nInitial probabilities: %s'%self.initial_probabilities() 369 if self._emission_symbols is not None: 370 s += '\nEmission symbols: %s'%self._emission_symbols 371 return s 372 373 def _emission_symbols_to_IntList(self, obs): 374 """ 375 Internal function used to convert a list of emission symbols to an IntList. 376 377 INPUT: 378 379 - obs -- a list of objects 380 381 OUTPUT: 382 383 - an IntList 384 385 EXAMPLES:: 386 387 sage: m = hmm.DiscreteHiddenMarkovModel([[0.4,0.6],[0.1,0.9]], [[0.1,0.9],[0.5,0.5]], [.2,.8], ['smile', 'frown']) 388 sage: m._emission_symbols_to_IntList(['frown','smile']) 389 [1, 0] 390 """ 391 d = self._emission_symbols_dict 392 return IntList([d[x] for x in obs]) 393 394 def _IntList_to_emission_symbols(self, obs): 395 """ 396 Internal function used to convert a list of emission symbols to an IntList. 397 398 INPUT: 399 400 - obs -- a list of objects 401 402 OUTPUT: 403 404 - an IntList 405 406 EXAMPLES:: 407 408 sage: m = hmm.DiscreteHiddenMarkovModel([[0.4,0.6],[0.1,0.9]], [[0.1,0.9],[0.5,0.5]], [.2,.8], ['smile', 'frown']) 409 sage: m._IntList_to_emission_symbols(stats.IntList([0,0,1,0])) 410 ['smile', 'smile', 'frown', 'smile'] 411 """ 412 d = self._emission_symbols 413 return [d[x] for x in obs] 414 415 def log_likelihood(self, obs, bint scale=True): 416 """ 417 Return the logarithm of the probability that this model produced the given 418 observation sequence. Thus the output is a non-positive number. 419 420 INPUT: 421 422 - obs -- sequence of observations 423 424 - scale -- boolean (default: True); if True, use rescaling 425 to overoid loss of precision due to the very limited 426 dynamic range of floats. You should leave this as True 427 unless the obs sequence is very small. 428 429 EXAMPLES:: 430 431 sage: m = hmm.DiscreteHiddenMarkovModel([[0.4,0.6],[0.1,0.9]], [[0.1,0.9],[0.5,0.5]], [.2,.8]) 432 sage: m.log_likelihood([0, 1, 0, 1, 1, 0, 1, 0, 0, 0]) 433 -7.3301308009370825 434 sage: m.log_likelihood([0, 1, 0, 1, 1, 0, 1, 0, 0, 0], scale=False) 435 -7.3301308009370816 436 sage: m.log_likelihood([]) 437 0.0 438 439 sage: m = hmm.DiscreteHiddenMarkovModel([[0.4,0.6],[0.1,0.9]], [[0.1,0.9],[0.5,0.5]], [.2,.8], ['happy','sad']) 440 sage: m.log_likelihood(['happy','happy']) 441 -1.6565295199679506 442 sage: m.log_likelihood(['happy','sad']) 443 -1.4731602941415523 444 445 Overflow from not using the scale option:: 446 447 sage: m = hmm.DiscreteHiddenMarkovModel([[0.4,0.6],[0.1,0.9]], [[0.1,0.9],[0.5,0.5]], [.2,.8]) 448 sage: m.log_likelihood([0,1]*1000, scale=True) 449 -1433.8206666527281 450 sage: m.log_likelihood([0,1]*1000, scale=False) 451 -inf 452 """ 453 if len(obs) == 0: 454 return 0.0 455 if self._emission_symbols is not None: 456 obs = self._emission_symbols_to_IntList(obs) 457 if not isinstance(obs, IntList): 458 obs = IntList(obs) 459 if scale: 460 return self._forward_scale(obs) 461 else: 462 return self._forward(obs) 463 464 def _forward(self, IntList obs): 465 """ 466 Memory-efficient implementation of the forward algorithm, without scaling. 467 468 INPUT: 469 470 - obs -- an integer list of observation states. 471 472 OUTPUT: 473 474 - float -- the log of the probability that the model produced this sequence 475 476 EXAMPLES:: 477 478 sage: m = hmm.DiscreteHiddenMarkovModel([[0.4,0.6],[0.1,0.9]], [[0.1,0.9],[0.5,0.5]], [.2,.8]) 479 sage: m._forward(stats.IntList([0,1]*10)) 480 -14.378663512219456 481 482 The forward algorithm computes the log likelihood:: 483 484 sage: m.log_likelihood(stats.IntList([0,1]*10), scale=False) 485 -14.378663512219456 486 487 But numerical overflow will happen (without scaling) for long sequences:: 488 489 sage: m._forward(stats.IntList([0,1]*1000)) 490 -inf 491 """ 492 if obs.max() > self.N or obs.min() < 0: 493 raise ValueError, "invalid observation sequence, since it includes unknown states" 494 495 cdef Py_ssize_t i, j, t, T = len(obs) 496 497 cdef TimeSeries alpha = TimeSeries(self.N), \ 498 alpha2 = TimeSeries(self.N) 499 500 # Initialization 501 for i in range(self.N): 502 alpha[i] = self.pi[i] * self.B[i*self.n_out + obs._values[0]] 503 alpha[i] = self.pi[i] * self.B[i*self.n_out + obs._values[0]] 504 505 # Induction 506 cdef double s 507 for t in range(1, T): 508 for j in range(self.N): 509 s = 0 510 for i in range(self.N): 511 s += alpha._values[i] * self.A._values[i*self.N + j] 512 alpha2._values[j] = s * self.B._values[j*self.n_out+obs._values[t]] 513 for j in range(self.N): 514 alpha._values[j] = alpha2._values[j] 515 516 # Termination 517 return log(alpha.sum()) 518 519 def _forward_scale(self, IntList obs): 520 """ 521 Memory-efficient implementation of the forward algorithm, with scaling. 522 523 INPUT: 524 525 - obs -- an integer list of observation states. 526 527 OUTPUT: 528 529 - float -- the log of the probability that the model produced this sequence 530 531 EXAMPLES:: 532 533 sage: m = hmm.DiscreteHiddenMarkovModel([[0.4,0.6],[0.1,0.9]], [[0.1,0.9],[0.5,0.5]], [.2,.8]) 534 sage: m._forward_scale(stats.IntList([0,1]*10)) 535 -14.378663512219454 536 537 The forward algorithm computes the log likelihood:: 538 539 sage: m.log_likelihood(stats.IntList([0,1]*10)) 540 -14.378663512219454 541 542 Note that the scale algorithm ensures that numerical overflow 543 won't happen for long sequences like it does with the forward 544 non-scaled algorithm:: 545 546 sage: m._forward_scale(stats.IntList([0,1]*1000)) 547 -1433.8206666527281 548 sage: m._forward(stats.IntList([0,1]*1000)) 549 -inf 550 551 A random sequence produced by the model is more likely:: 552 553 sage: set_random_seed(0); v = m.sample(1000) 554 sage: m._forward_scale(v) 555 -686.87531893650555 556 """ 557 # This is just like self._forward(obs) above, except at every step of the 558 # algorithm, we rescale the vector alpha so that the sum of 559 # the entries in alpha is 1. Then, at the end of the 560 # algorithm, the sum of probabilities (alpha.sum()) is of 561 # course just 1. However, the true probability that we want 562 # is the product of the numbers that we divided by when 563 # rescaling, since at each step of the iteration the next term 564 # depends linearly on the previous one. Instead of returning 565 # the product, we return the sum of the logs, which avoid 566 # numerical error. 567 cdef Py_ssize_t i, j, t, T = len(obs) 568 569 # The running some of the log probabilities 570 cdef double log_probability = 0, sum, a 571 572 cdef TimeSeries alpha = TimeSeries(self.N), \ 573 alpha2 = TimeSeries(self.N) 574 575 # Initialization 576 sum = 0 577 for i in range(self.N): 578 a = self.pi[i] * self.B[i*self.n_out + obs._values[0]] 579 alpha[i] = a 580 sum += a 581 582 log_probability = log(sum) 583 alpha.rescale(1/sum) 584 585 # Induction 586 # The code below is just an optimized version of: 587 # alpha = (alpha * A).pairwise_product(B[O[t+1]]) 588 # alpha = alpha.scale(1/alpha.sum()) 589 # along with keeping track of the log of the scaling factor. 590 cdef double s 591 for t in range(1, T): 592 sum = 0 593 for j in range(self.N): 594 s = 0 595 for i in range(self.N): 596 s += alpha._values[i] * self.A._values[i*self.N + j] 597 a = s * self.B._values[j*self.n_out+obs._values[t]] 598 alpha2._values[j] = a 599 sum += a 600 601 log_probability += log(sum) 602 for j in range(self.N): 603 alpha._values[j] = alpha2._values[j] / sum 604 605 # Termination 606 return log_probability 607 608 def generate_sequence(self, Py_ssize_t length): 609 """ 610 Return a sample of the given length from this HMM. 611 612 INPUT: 613 614 - length -- positive integer 615 616 OUTPUT: 617 618 - an IntList or list of emission symbols 619 - IntList of the actual states the model was in when 620 emitting the corresponding symbols 621 622 EXAMPLES:: 623 624 In this example, the emission symbols are not set:: 625 626 sage: set_random_seed(0) 627 sage: a = hmm.DiscreteHiddenMarkovModel([[0.1,0.9],[0.1,0.9]], [[1,0],[0,1]], [0,1]) 628 sage: a.generate_sequence(5) 629 ([1, 0, 1, 1, 1], [1, 0, 1, 1, 1]) 630 sage: list(a.generate_sequence(1000)[0]).count(0) 631 90 632 633 Here the emission symbols are set:: 634 635 sage: set_random_seed(0) 636 sage: a = hmm.DiscreteHiddenMarkovModel([[0.5,0.5],[0.1,0.9]], [[1,0],[0,1]], [0,1], ['up', 'down']) 637 sage: a.generate_sequence(5) 638 (['down', 'up', 'down', 'down', 'down'], [1, 0, 1, 1, 1]) 639 """ 640 if length < 0: 641 raise ValueError, "length must be nonnegative" 642 643 # Create Integer lists for states and observations 644 cdef IntList states = IntList(length) 645 cdef IntList obs = IntList(length) 646 if length == 0: 647 # A special case 648 if self._emission_symbols is None: 649 return states, obs 650 else: 651 return states, [] 652 653 # Setup variables, including random state. 654 cdef Py_ssize_t i, j 655 cdef randstate rstate = current_randstate() 656 cdef int q = 0 657 cdef double r, accum 658 r = rstate.c_rand_double() 659 660 # Now choose random variables from our discrete distribution. 661 662 # This standard naive algorithm has complexity that is linear 663 # in the number of states. It might be possible to replace it 664 # by something that is more efficient. However, make sure to 665 # refactor this code into distribution.pyx before doing so. 666 # Note that state switching involves the same algorithm 667 # (below). Use GSL as described here to make this O(1): 668 # http://www.gnu.org/software/gsl/manual/html_node/General-Discrete-Distributions.html 669 670 # Choose initial state: 671 accum = 0 672 for i in range(self.N): 673 if r < self.pi._values[i] + accum: 674 q = i 675 else: 676 accum += self.pi._values[i] 677 678 states._values[0] = q 679 # Generate a symbol from state q 680 obs._values[0] = self._gen_symbol(q, rstate.c_rand_double()) 681 682 cdef double* row 683 cdef int O 684 _sig_on 685 for i in range(1, length): 686 # Choose next state 687 accum = 0 688 row = self.A._values + q*self.N 689 r = rstate.c_rand_double() 690 for j in range(self.N): 691 if r < row[j] + accum: 692 q = j 693 break 694 else: 695 accum += row[j] 696 states._values[i] = q 697 # Generate symbol from this new state q 698 obs._values[i] = self._gen_symbol(q, rstate.c_rand_double()) 699 _sig_off 700 701 if self._emission_symbols is None: 702 # No emission symbol mapping 703 return obs, states 704 else: 705 # Emission symbol mapping, so change our intlist into a list of symbols 706 return self._IntList_to_emission_symbols(obs), states 707 708 cdef int _gen_symbol(self, int q, double r): 709 """ 710 Generate a symbol in state q using the randomly chosen 711 floating point number r, which should be between 0 and 1. 712 713 INPUT: 714 715 - q -- a nonnegative integer, which specifies a state 716 717 - r -- a real number between 0 and 1 718 719 OUTPUT: 720 721 - a nonnegative int 722 723 EXAMPLES: 724 725 sage: m = hmm.DiscreteHiddenMarkovModel([[0.4,0.6],[0.1,0.9]], [[0.1,0.9],[0.9,0.1]], [.2,.8]) 726 sage: set_random_seed(0) 727 sage: m.generate_sequence(10) # indirect test 728 ([0, 1, 0, 0, 0, 0, 1, 0, 1, 1], [1, 0, 1, 1, 1, 1, 0, 0, 0, 0]) 729 730 """ 731 cdef Py_ssize_t j 732 cdef double a, accum = 0 733 # See the comments above about switching to GSL for this; also note 734 # that this should get factored out into a call to something 735 # defined in distributions.pyx. 736 for j in range(self.n_out): 737 a = self.B._values[q*self.N+j] 738 if r < a + accum: 739 return j 740 else: 741 accum += a 742 # None of the values was selected: shouldn't this only happen if the 743 # distribution is invalid? Anyway, this will get factored out. 744 return self.n_out - 1 745 746 def viterbi(self, obs, log_scale=True): 747 """ 748 Determine "the" hidden sequence of states that is most likely 749 to produce the given sequence seq of observations, along with 750 the probability that this hidden sequence actually produced 751 the observation. 752 753 INPUT: 754 755 - seq -- sequence of emitted ints or symbols 756 757 - log_scale -- bool (default: True) whether to scale the 758 sequence in order to avoid numerical overflow. 759 760 OUTPUT: 761 762 - list -- "the" most probable sequence of hidden states, i.e., 763 the Viterbi path. 764 765 - float -- log of probability that the observed sequence 766 was produced by the Viterbi sequence of states. 767 768 EXAMPLES:: 769 770 sage: a = hmm.DiscreteHiddenMarkovModel([[0.1,0.9],[0.1,0.9]], [[0.9,0.1],[0.1,0.9]], [0.5,0.5]) 771 sage: a.viterbi([1,0,0,1,0,0,1,1]) 772 ([1, 0, 0, 1, ..., 0, 1, 1], -11.06245322477221...) 773 774 We predict the state sequence when the emissions are 3/4 and 'abc'.:: 775 776 sage: a = hmm.DiscreteHiddenMarkovModel([[0.1,0.9],[0.1,0.9]], [[0.9,0.1],[0.1,0.9]], [0.5,0.5], [3/4, 'abc']) 777 778 Note that state 0 is common below, despite the model trying hard to 779 switch to state 1:: 780 781 sage: a.viterbi([3/4, 'abc', 'abc'] + [3/4]*10) 782 ([0, 1, 1, 0, 0 ... 0, 0, 0, 0, 0], -25.299405845367794) 783 """ 784 if self._emission_symbols is not None: 785 obs = self._emission_symbols_to_IntList(obs) 786 elif not isinstance(obs, IntList): 787 obs = IntList(obs) 788 if log_scale: 789 return self._viterbi_scale(obs) 790 else: 791 return self._viterbi(obs) 792 793 cpdef _viterbi(self, IntList obs): 794 """ 795 Used internally to compute the viterbi path, without 796 rescaling. This can be useful for short sequences. 797 798 INPUT: 799 800 - obs -- IntList 801 802 OUTPUT: 803 804 - IntList (most likely state sequence) 805 806 - log of probability that the observed sequence was 807 produced by the Viterbi sequence of states. 808 809 EXAMPLES:: 810 811 sage: m = hmm.DiscreteHiddenMarkovModel([[0.1,0.9],[0.9,0.1]], [[1,0],[0,1]], [.2,.8]) 812 sage: m._viterbi(stats.IntList([1]*5)) 813 ([1, 1, 1, 1, 1], -9.4334839232903924) 814 sage: m._viterbi(stats.IntList([0]*5)) 815 ([0, 0, 0, 0, 0], -10.819778284410283) 816 817 Long sequences will overflow:: 818 819 sage: m._viterbi(stats.IntList([0]*1000)) 820 ([0, 0, 0, 0, 0 ... 0, 0, 0, 0, 0], -inf) 821 """ 822 cdef Py_ssize_t t, T = obs._length 823 cdef IntList state_sequence = IntList(T) 824 if T == 0: 825 return state_sequence, 0.0 826 827 cdef int i, j, N = self.N 828 829 # delta[i] is the maximum of the probabilities over all 830 # paths ending in state i. 831 cdef TimeSeries delta = TimeSeries(N), delta_prev = TimeSeries(N) 832 833 # We view psi as an N x T matrix of ints. The quantity 834 # psi[N*t + j] 835 # is a most probable hidden state at time t-1, given 836 # that j is a most probable state at time j. 837 cdef IntList psi = IntList(N * T) # initialized to 0 by default 838 839 # Initialization: 840 for i in range(N): 841 delta._values[i] = self.pi._values[i] * self.B._values[self.n_out*i + obs._values[0]] 842 843 # Recursion: 844 cdef double mx, tmp 845 cdef int index 846 for t in range(1, T): 847 delta_prev, delta = delta, delta_prev 848 for j in range(N): 849 # delta_t[j] = max_i(delta_{t-1}(i) a_{i,j}) * b_j(obs[t]) 850 mx = -1 851 index = -1 852 for i in range(N): 853 tmp = delta_prev._values[i]*self.A._values[i*N+j] 854 if tmp > mx: 855 mx = tmp 856 index = i 857 delta._values[j] = mx * self.B._values[self.n_out*j + obs._values[t]] 858 psi._values[N*t + j] = index 859 860 # Termination: 861 mx, index = delta.max(index=True) 862 863 # Path (state sequence) backtracking: 864 state_sequence._values[T-1] = index 865 t = T-2 866 while t >= 0: 867 state_sequence._values[t] = psi._values[N*(t+1) + state_sequence._values[t+1]] 868 t -= 1 869 870 return state_sequence, log(mx) 871 872 873 cpdef _viterbi_scale(self, IntList obs): 874 """ 875 Used internally to compute the viterbi path with rescaling. 876 877 INPUT: 878 879 - obs -- IntList 880 881 OUTPUT: 882 883 - IntList (most likely state sequence) 884 885 - log of probability that the observed sequence was 886 produced by the Viterbi sequence of states. 887 888 EXAMPLES:: 889 890 sage: m = hmm.DiscreteHiddenMarkovModel([[0.1,0.9],[0.9,0.1]], [[.5,.5],[0,1]], [.2,.8]) 891 sage: m._viterbi_scale(stats.IntList([1]*10)) 892 ([1, 0, 1, 0, 1, 0, 1, 0, 1, 0], -4.6371240950343733) 893 894 Long sequences should not overflow:: 895 896 sage: m._viterbi_scale(stats.IntList([1]*1000)) 897 ([1, 0, 1, 0, 1 ... 0, 1, 0, 1, 0], -452.05188897345...) 898 """ 899 # The algorithm is the same as _viterbi above, except 900 # we take the logarithms of everything first, and add 901 # instead of multiply. 902 cdef Py_ssize_t t, T = obs._length 903 cdef IntList state_sequence = IntList(T) 904 if T == 0: 905 return state_sequence, 0.0 906 cdef int i, j, N = self.N 907 908 # delta[i] is the maximum of the probabilities over all 909 # paths ending in state i. 910 cdef TimeSeries delta = TimeSeries(N), delta_prev = TimeSeries(N) 911 912 # We view psi as an N x T matrix of ints. The quantity 913 # psi[N*t + j] 914 # is a most probable hidden state at time t-1, given 915 # that j is a most probable state at time j. 916 cdef IntList psi = IntList(N * T) # initialized to 0 by default 917 918 # Log Preprocessing 919 cdef TimeSeries A = self.A.log() 920 cdef TimeSeries B = self.B.log() 921 cdef TimeSeries pi = self.pi.log() 922 923 # Initialization: 924 for i in range(N): 925 delta._values[i] = pi._values[i] + B._values[self.n_out*i + obs._values[0]] 926 927 # Recursion: 928 cdef double mx, tmp, minus_inf = float('-inf') 929 cdef int index 930 931 for t in range(1, T): 932 delta_prev, delta = delta, delta_prev 933 for j in range(N): 934 # delta_t[j] = max_i(delta_{t-1}(i) a_{i,j}) * b_j(obs[t]) 935 mx = minus_inf 936 index = -1 937 for i in range(N): 938 tmp = delta_prev._values[i] + A._values[i*N+j] 939 if tmp > mx: 940 mx = tmp 941 index = i 942 delta._values[j] = mx + B._values[self.n_out*j + obs._values[t]] 943 psi._values[N*t + j] = index 944 945 # Termination: 946 mx, index = delta.max(index=True) 947 948 # Path (state sequence) backtracking: 949 state_sequence._values[T-1] = index 950 t = T-2 951 while t >= 0: 952 state_sequence._values[t] = psi._values[N*(t+1) + state_sequence._values[t+1]] 953 t -= 1 954 955 return state_sequence, mx 956 957 cdef TimeSeries _backward_scale_all(self, IntList obs, TimeSeries scale): 958 r""" 959 Return the scaled matrix of values `\beta_t(i)` that appear in 960 the backtracking algorithm. This function is used internally 961 by the Baum-Welch algorithm. 962 963 The matrix is stored as a TimeSeries T, such that 964 `\beta_t(i) = T[t*N + i]` where N is the number of states of 965 the Hidden Markov Model. 966 967 The quantity beta_t(i) is the probability of observing the 968 sequence obs[t+1:] assuming that the model is in state i at 969 time t. 970 971 INPUT: 972 973 - obs -- IntList 974 - scale -- series that is *changed* in place, so that 975 after calling this function, scale[t] is value that is 976 used to scale each of the `\beta_t(i)`. 977 978 OUTPUT: 979 980 - a TimeSeries of values beta_t(i). 981 - the input object scale is modified 982 """ 983 cdef Py_ssize_t t, T = obs._length 984 cdef int N = self.N, i, j 985 cdef double s 986 cdef TimeSeries beta = TimeSeries(N*T, initialize=False) 987 988 # 1. Initialization 989 for i in range(N): 990 beta._values[(T-1)*N + i] = 1/scale._values[T-1] 991 992 # 2. Induction 993 t = T-2 994 while t >= 0: 995 for i in range(N): 996 s = 0 997 for j in range(N): 998 s += self.A._values[i*N+j] * \ 999 self.B._values[j*self.n_out+obs._values[t+1]] * beta._values[(t+1)*N+j] 1000 beta._values[t*N + i] = s/scale._values[t] 1001 t -= 1 1002 return beta 1003 1004 cdef _forward_scale_all(self, IntList obs): 1005 """ 1006 Return scaled values alpha_t(i), the sequence of scalings, and 1007 the log probability. 1008 1009 INPUT: 1010 1011 - obs -- IntList 1012 1013 OUTPUT: 1014 1015 - TimeSeries alpha with alpha_t(i) = alpha[t*N + i] 1016 - TimeSeries scale with scale[t] the scaling at step t 1017 - float -- log_probability of the observation sequence 1018 being produced by the model. 1019 """ 1020 cdef Py_ssize_t i, j, t, T = len(obs) 1021 cdef int N = self.N 1022 1023 # The running some of the log probabilities 1024 cdef double log_probability = 0, sum, a 1025 1026 cdef TimeSeries alpha = TimeSeries(N*T, initialize=False) 1027 cdef TimeSeries scale = TimeSeries(T, initialize=False) 1028 1029 # Initialization 1030 sum = 0 1031 for i in range(self.N): 1032 a = self.pi._values[i] * self.B._values[i*self.n_out + obs._values[0]] 1033 alpha._values[0*N + i] = a 1034 sum += a 1035 1036 scale._values[0] = sum 1037 log_probability = log(sum) 1038 for i in range(self.N): 1039 alpha._values[0*N + i] /= sum 1040 1041 # Induction 1042 # The code below is just an optimized version of: 1043 # alpha = (alpha * A).pairwise_product(B[O[t+1]]) 1044 # alpha = alpha.scale(1/alpha.sum()) 1045 # along with keeping track of the log of the scaling factor. 1046 cdef double s 1047 for t in range(1,T): 1048 sum = 0 1049 for j in range(N): 1050 s = 0 1051 for i in range(N): 1052 s += alpha._values[(t-1)*N + i] * self.A._values[i*N + j] 1053 a = s * self.B._values[j*self.n_out + obs._values[t]] 1054 alpha._values[t*N + j] = a 1055 sum += a 1056 1057 log_probability += log(sum) 1058 scale._values[t] = sum 1059 for j in range(self.N): 1060 alpha._values[t*N + j] /= sum 1061 1062 # Termination 1063 return alpha, scale, log_probability 1064 1065 cdef TimeSeries _baum_welch_gamma(self, TimeSeries alpha, TimeSeries beta): 1066 """ 1067 Used internally to compute the scaled quantity gamma_t(j) 1068 appearing in the Baum-Welch reestimation algorithm. 1069 1070 The quantity gamma_t(j) is the (scaled!) probability of being 1071 in state j at time t, given the observation sequence. 1072 1073 INPUT: 1074 1075 - alpha -- TimeSeries as output by the scaled forward algorithm 1076 - beta -- TimeSeries as output by the scaled backward algorithm 1077 1078 OUTPUT: 1079 1080 - TimeSeries gamma such that gamma[t*N+j] is gamma_t(j). 1081 """ 1082 cdef int j, N = self.N 1083 cdef Py_ssize_t t, T = alpha._length//N 1084 cdef TimeSeries gamma = TimeSeries(alpha._length, initialize=False) 1085 cdef double denominator 1086 _sig_on 1087 for t in range(T): 1088 denominator = 0 1089 for j in range(N): 1090 gamma._values[t*N + j] = alpha._values[t*N + j] * beta._values[t*N + j] 1091 denominator += gamma._values[t*N + j] 1092 for j in range(N): 1093 gamma._values[t*N + j] /= denominator 1094 _sig_off 1095 return gamma 1096 1097 cdef TimeSeries _baum_welch_xi(self, TimeSeries alpha, TimeSeries beta, IntList obs): 1098 """ 1099 Used internally to compute the scaled quantity xi_t(i,j) 1100 appearing in the Baum-Welch reestimation algorithm. 1101 1102 INPUT: 1103 1104 - alpha -- TimeSeries as output by the scaled forward algorithm 1105 - beta -- TimeSeries as output by the scaled backward algorithm 1106 - obs -- IntList of observations 1107 1108 OUTPUT: 1109 1110 - TimeSeries xi such that xi[t*N*N + i*N + j] = xi_t(i,j). 1111 """ 1112 cdef int i, j, N = self.N 1113 cdef double sum 1114 cdef Py_ssize_t t, T = alpha._length//N 1115 cdef TimeSeries xi = TimeSeries(T*N*N, initialize=False) 1116 _sig_on 1117 for t in range(T-1): 1118 sum = 0.0 1119 for i in range(N): 1120 for j in range(N): 1121 xi._values[t*N*N + i*N + j] = alpha._values[t*N + i]*beta._values[(t+1)*N + j]*\ 1122 self.A._values[i*N + j] * self.B._values[j*self.n_out + obs._values[t+1]] 1123 sum += xi._values[t*N*N + i*N + j] 1124 for i in range(N): 1125 for j in range(N): 1126 xi._values[t*N*N + i*N + j] /= sum 1127 _sig_off 1128 return xi 1129 1130 def baum_welch(self, obs, int max_iter=100, double log_likelihood_cutoff=1e-4, bint fix_emissions=False): 1131 """ 1132 Given an observation sequence obs, improve this HMM using the 1133 Baum-Welch algorithm to increase the probability of observing obs. 1134 1135 INPUT: 1136 1137 - obs -- list of emissions 1138 1139 - max_iter -- integer (default: 100) maximum number 1140 of Baum-Welch steps to take 1141 1142 - log_likehood_cutoff -- positive float (default: 1e-4); 1143 the minimal improvement in likelihood with respect to 1144 the last iteration required to continue. Relative value 1145 to log likelihood. 1146 1147 - fix_emissions -- bool (default: False); if True, do not 1148 change emissions when updating 1149 1150 OUTPUT: 1151 1152 - changes the model in places, and returns the log 1153 likelihood and number of iterations. 1154 1155 EXAMPLES:: 1156 1157 sage: m = hmm.DiscreteHiddenMarkovModel([[0.1,0.9],[0.9,0.1]], [[.5,.5],[0,1]], [.2,.8]) 1158 sage: m.baum_welch([1,0]*20, log_likelihood_cutoff=0) 1159 (0.0, 4) 1160 sage: m 1161 Discrete Hidden Markov Model with 2 States and 2 Emissions 1162 Transition matrix: 1163 [1.35152697077e-51 1.0] 1164 [ 1.0 0.0] 1165 Emission matrix: 1166 [ 1.0 6.46253713885e-52] 1167 [ 0.0 1.0] 1168 Initial probabilities: [0.0000, 1.0000] 1169 1170 The following illustrates how Baum-Welch is only a local 1171 optimizer, i.e., the above model is far more likely to produce 1172 the sequence [1,0]*20 than the one we get below:: 1173 1174 sage: m = hmm.DiscreteHiddenMarkovModel([[0.5,0.5],[0.5,0.5]], [[.5,.5],[.5,.5]], [.5,.5]) 1175 sage: m.baum_welch([1,0]*20, log_likelihood_cutoff=0) 1176 (-27.725887222397784, 1) 1177 sage: m 1178 Discrete Hidden Markov Model with 2 States and 2 Emissions 1179 Transition matrix: 1180 [0.5 0.5] 1181 [0.5 0.5] 1182 Emission matrix: 1183 [0.5 0.5] 1184 [0.5 0.5] 1185 Initial probabilities: [0.5000, 0.5000] 1186 1187 We illustrate fixing emissions:: 1188 1189 sage: m = hmm.DiscreteHiddenMarkovModel([[0.1,0.9],[0.9,0.1]], [[.5,.5],[.2,.8]], [.2,.8]) 1190 sage: set_random_seed(0); v = m.sample(100) 1191 sage: m.baum_welch(v,fix_emissions=True) 1192 (-66.986308569187742, 100) 1193 sage: m.emission_matrix() 1194 [0.5 0.5] 1195 [0.2 0.8] 1196 sage: m = hmm.DiscreteHiddenMarkovModel([[0.1,0.9],[0.9,0.1]], [[.5,.5],[.2,.8]], [.2,.8]) 1197 sage: m.baum_welch(v) 1198 (-66.7823606592935..., 100) 1199 sage: m.emission_matrix() 1200 [0.530308574863 0.469691425137] 1201 [0.290977555017 0.709022444983] 1202 """ 1203 if self._emission_symbols is not None: 1204 obs = self._emission_symbols_to_IntList(obs) 1205 elif not isinstance(obs, IntList): 1206 obs = IntList(obs) 1207 cdef IntList _obs = obs 1208 cdef TimeSeries alpha, beta, scale, gamma, xi 1209 cdef double log_probability, log_probability_prev, delta 1210 cdef int i, j, k, N, n_iterations 1211 cdef Py_ssize_t t, T 1212 cdef double denominator_A, numerator_A, denominator_B, numerator_B 1213 1214 # Initialization 1215 alpha, scale, log_probability = self._forward_scale_all(_obs) 1216 beta = self._backward_scale_all(_obs, scale) 1217 gamma = self._baum_welch_gamma(alpha, beta) 1218 xi = self._baum_welch_xi(alpha, beta, _obs) 1219 log_probability_prev = log_probability 1220 N = self.N 1221 n_iterations = 0 1222 T = len(_obs) 1223 1224 # Re-estimation 1225 while True: 1226 1227 # Reestimate frequency of state i in time t=0 1228 for i in range(N): 1229 self.pi._values[i] = gamma._values[0*N+i] 1230 1231 # Reestimate transition matrix and emissions probability in 1232 # each state. 1233 for i in range(N): 1234 denominator_A = 0.0 1235 for t in range(T-1): 1236 denominator_A += gamma._values[t*N+i] 1237 for j in range(N): 1238 numerator_A = 0.0 1239 for t in range(T-1): 1240 numerator_A += xi._values[t*N*N+i*N+j] 1241 self.A._values[i*N+j] = numerator_A / denominator_A 1242 1243 if not fix_emissions: 1244 denominator_B = denominator_A + gamma._values[(T-1)*N + i] 1245 for k in range(self.n_out): 1246 numerator_B = 0.0 1247 for t in range(T): 1248 if _obs._values[t] == k: 1249 numerator_B += gamma._values[t*N + i] 1250 self.B._values[i*self.n_out + k] = numerator_B / denominator_B 1251 1252 # Initialization 1253 alpha, scale, log_probability = self._forward_scale_all(_obs) 1254 beta = self._backward_scale_all(_obs, scale) 1255 gamma = self._baum_welch_gamma(alpha, beta) 1256 xi = self._baum_welch_xi(alpha, beta, _obs) 1257 1258 # Compute the difference between the log probability of 1259 # two iterations. 1260 delta = log_probability - log_probability_prev 1261 log_probability_prev = log_probability 1262 n_iterations += 1 1263 1264 # If the log probability does not change by more than 1265 # delta, then terminate 1266 if delta <= log_likelihood_cutoff or n_iterations >= max_iter: 1267 break 1268 1269 return log_probability, n_iterations 1270 1271 1272 # Keep this -- it's for backwards compatibility with the GHMM based implementation 1273 def unpickle_discrete_hmm_v0(A, B, pi, emission_symbols, name): 1274 """ 1275 TESTS:: 1276 1277 sage: m = hmm.DiscreteHiddenMarkovModel([[0.4,0.6],[0.1,0.9]], [[0.0,1.0],[0.5,0.5]], [1,0]) 1278 sage: sage.stats.hmm.hmm.unpickle_discrete_hmm_v0(m.transition_matrix(), m.emission_matrix(), m.initial_probabilities(), ['a','b'], 'test model') 1279 Discrete Hidden Markov Model with 2 States and 2 Emissions... 1280 """ 1281 return DiscreteHiddenMarkovModel(A,B,pi,emission_symbols,normalize=False) 1282 1283 def unpickle_discrete_hmm_v1(A, B, pi, n_out, emission_symbols, emission_symbols_dict): 1284 """ 1285 TESTS:: 1286 1287 sage: m = hmm.DiscreteHiddenMarkovModel([[0.4,0.6],[0.1,0.9]], [[0.0,1.0],[0.5,0.5]], [1,0],['a','b']) 1288 sage: loads(dumps(m)) == m # indirect test 1289 True 1290 """ 1291 cdef DiscreteHiddenMarkovModel m = PY_NEW(DiscreteHiddenMarkovModel) 1292 m.A = A 1293 m.B = B 1294 m.pi = pi 1295 m.n_out = n_out 1296 m._emission_symbols = emission_symbols 1297 m._emission_symbols_dict = emission_symbols_dict 1298 return m 1299 1300 -
deleted file sage/stats/hmm/misc.pxi
diff --git a/sage/stats/hmm/misc.pxi b/sage/stats/hmm/misc.pxi deleted file mode 100644
+ - 1 #############################################################################2 # Copyright (C) 2008 William Stein <wstein@gmail.com>3 # Distributed under the terms of the GNU General Public License (GPL),4 # version 2 or any later version at your option.5 # The full text of the GPL is available at:6 # http://www.gnu.org/licenses/7 #############################################################################8 9 cdef extern from "string.h":10 char *strcpy(char *s1, char *s2)11 void* memcpy(void* dst, void* src, size_t len)12 13 cdef double* to_double_array(v) except NULL:14 """15 Transform a Python list of floats to a C array of doubles. The caller is16 responsible for deallocating the resulting memory.17 18 INPUT:19 v -- a list of objects coercible to floats20 OUTPUT:21 a newly allocated C array of doubles22 """23 cdef double x24 cdef double* w = <double*> safe_malloc(sizeof(double)*len(v))25 cdef Py_ssize_t i = 026 for x in v:27 w[i] = x28 i += 129 return w30 31 cdef int* to_int_array(v) except NULL:32 """33 Transform a Python list of ints to a C array of ints. The caller is34 responsible for deallocating the resulting memory.35 36 INPUT:37 v -- a list of objects coercible to ints38 OUTPUT:39 a newly allocated C array of ints40 """41 cdef int x42 cdef int* w = <int*> safe_malloc(sizeof(int)*len(v))43 cdef Py_ssize_t i = 044 for x in v:45 w[i] = x46 i += 147 return w48 49 cdef void* safe_malloc(int bytes) except NULL:50 """51 malloc the given bytes of memory and check that the malloc52 succeeds -- if not raise a MemoryError.53 54 INPUT:55 bytes -- a nonnegatie integer56 57 OUTPUT:58 void pointer or raise a MemoryError.59 """60 cdef void* t = sage_malloc(bytes)61 if not t:62 raise MemoryError, "error allocating memory for Hidden Markov Model"63 return t64 -
new file sage/stats/hmm/util.pxd
diff --git a/sage/stats/hmm/util.pxd b/sage/stats/hmm/util.pxd new file mode 100644
- + 1 from sage.finance.time_series cimport TimeSeries 2 3 cdef class HMM_Util: 4 cpdef normalize_probability_TimeSeries(self, TimeSeries T, Py_ssize_t i, Py_ssize_t j) 5 cpdef TimeSeries initial_probs_to_TimeSeries(self, pi, bint normalize) 6 cpdef TimeSeries state_matrix_to_TimeSeries(self, A, int N, bint normalize) 7 -
new file sage/stats/hmm/util.pyx
diff --git a/sage/stats/hmm/util.pyx b/sage/stats/hmm/util.pyx new file mode 100644
- + 1 """ 2 Hidden Markov Models -- Utility functions 3 4 AUTHOR: 5 6 - William Stein, 2010-03 7 """ 8 9 ############################################################################# 10 # Copyright (C) 2010 William Stein <wstein@gmail.com> 11 # Distributed under the terms of the GNU General Public License (GPL) v2+. 12 # The full text of the GPL is available at: 13 # http://www.gnu.org/licenses/ 14 ############################################################################# 15 16 include "../../ext/stdsage.pxi" 17 18 from sage.matrix.all import is_Matrix 19 from sage.misc.flatten import flatten 20 21 cdef class HMM_Util: 22 """ 23 A class used in order to share cdef's methods between different files. 24 """ 25 cpdef normalize_probability_TimeSeries(self, TimeSeries T, Py_ssize_t i, Py_ssize_t j): 26 """ 27 This function is used internally by the Hidden Markov Models code. 28 29 Replace entries of T[i:j] in place so that they are all 30 nonnegative and sum to 1. Negative entries are replaced by 0 and 31 T[i:j] is then rescaled to ensure that the sum of the entries in 32 each row is equal to 1. If all entries are 0, replace them 33 by 1/(j-i). 34 35 INPUT: 36 37 - T -- a TimeSeries 38