This paper discusses a Bayesian approach to regularizing hidden Markov models and demonstrates an application of this scheme to Bioinformatics

Husmeier, D.

English

Enlighten: Publications

Bayesian regularization of hidden Markov models with an application to bioinformatics

Enlighten

       Husmeier, D. (2000) Bayesian regularization of hidden Markov models with an application to bioinformatics. Neural Network World, 10 (4). pp. 589-595. ISSN 1210-0552  Copyright © 2000 Akademie Ved Ceske Republiky  A copy can be downloaded for personal non-commercial research or study, without prior permission or charge   The content must not be changed in any way or reproduced in any format or medium without the formal permission of the copyright holder(s)   When referring to this work, full bibliographic details must be given      http://eprints.gla.ac.uk/85684/       Deposited on:  13 September 2013          Enlighten – Research publications by members of the University of Glasgow http://eprints.gla.ac.uk Bayesian Regularization of Hidden MarkovModels with an Application to BioinformaticsDirk HusmeierBiomathematics and Statistics Scotland (BioSS)SCRI, Dundee DD2 5DA, United KingdomEmail: dirk@bioss.sari.ac.ukThis paper discusses a Bayesian approach to regularizing hidden Markovmodels and demonstrates an application of this scheme to Bioinformatics.1 IntroductionHidden Markov models (HMMs) are close relatives of neural networks (NNs) in manyrespects. Both models are composed of a network of visible and hidden nodes connectedby edges, which can be described in the more general framework of probabilistic graphicalmodels [5]. The forward-backward algorithm for training HMMs [10] shows strikingresemblance to the backpropagation algorithm for NNs. Finally, many applications,especially in signal processing [2] and bioinformatics [1], employ hybrid schemes whichcombine HMMs and NNs. This article focuses on a problem common to both models.For sparse data, the classical training algorithms (mentioned above), which derive froma maximum likelihood (ML) approach, are sub-optimal due to overfitting. Following [6],much research on NNs in the last few years has explored Bayesian methods as a possibleremedy in this respect. The objective of this article is to take a similar route and testif the generalization performance of HMMs can be improved by Bayesian free energyminimization.2 Classical HMMsConsider a vector of observations, yt, where t represents some ordered label, e.g., timein signal processing, or the site of a DNA strand in bioinformatics (see below). Corre-1s_Ns_2...s_1y_2y_1 y_Nwνs_Ns_2...s_1y_2y_1 y_NwβανFigure 1: Left: HMM drawn as a probabilistic graphical model, where nodes repre-sent random variables and arcs indicate conditional dependencies between these nodes.Empty circles show hidden states, filled circles represent vectors of observations. The pa-rameters of the model are shown as squares. Right: A Bayesian approach to parameteradaptation. The parameters of the model, ν and w, are themselves treated as randomvariables, whose prior distribution depends on further so-called hyperparameters (shownas squares).lations between these observations are assumed to be indirect via unobserved, discretehidden states st ∈ {1, . . . , K}. The model is characterized by two sets of conditionalindependence relations [10]: (1) The observation at time or site t, yt, is independentof all other observations and states given st. (2) The state at t, st, is independent ofthe previous sequence of states given st−1. Using these independence relations, the jointprobability for the sequence of states and observations can be factored asP (Y, s) = P (y1, . . . ,yN , s1, . . . , sN) = P (s1)N∏t=2P (st|st−1)N∏t=1P (yt|st) (1)The corresponding graphical model is shown in Figure 1. In general, the transitionprobabilities P (st|st−1) and the emission probabilities P (yt|st) depend on some modelparameters ν andw. A standard training scheme is to optimize these parameters so as tomaximize the likelihood L = P (Y|ν,w) =∑sP (Y, s|ν,w). This can be accomplishedwith the forward-backward and the Baum-Welch algorithm, as described in [10].23 A Bayesian approach3.1 General outline of the methodologyA possible shortcoming of the maximum likelihood (ML) method is its susceptibility tooverfitting when the data are sparse. We therefore adopt a Bayesian approach, whichhas proven to improve the generalization performance of neural networks [6], [9], andtreat the parameters ν and w as random variables. Starting from some prior distribu-tion P (ν,w), the ultimate objective of learning is to determine the posterior distributionP (s,w, ν|Y). As this is analytically intractable, we need to make certain approxima-tions. We can sample from the posterior distribution numerically with Markov chainMonte Carlo (MCMC). A good overview of this method can be found in [4]. Alterna-tively, we can assume a simpler, analytically tractable distribution Q(s,w, ν) and adaptits parameters so as to minimize the Kullback-Leibler divergence between Q and thetrue posterior. This so-called method of variational free energy minimization has be-come very popular in the context of learning in graphical models [5]. A first applicationto HMMs can be found in [7].We will constrain our approximating distribution to be separable, such that1Q(s,w, ν) = Q(s)Q(w)Q(ν) (2)The closeness of this distribution to the true posterior is measured by the Kullback-Leibler divergenceF (Q) =∑s∫dw∫dνQ(s,w, ν) ln[Q(s,w, ν)P (Y, s,w, ν)](3)whereP (Y, s,w, ν) =N∏t=2P (st|st−1, ν)N∏t=1P (yt|st,w)P (s1)P (ν)P (w) (4)and we have introduced P (ν) and P (w) as priors on ν and w. Our objective is to opti-mize Q(s,w, ν) so as to minimize F . This can be accomplished in an iterative scheme,where each of the marginal distributions Q(s), Q(w), Q(ν) is optimized separately whilekeeping the remaining two distributions fixed.1We use the convention that different arguments indicate different probability distributions.33.2 Explicit form of the update schemeWe here assume a specific form for the transition probabilities, first suggested in [8]:P (st|st−1) = νδ(st, st−1) +1− νK − 1[1− δ(st, st−1)] (5)This choice is motivated by the phylogenetic problem of detecting recombinations inDNA sequence alignments, where 1−ν can be interpreted as a recombination probability(see below). Since ν is a binomial random variable, the natural conjugate prior is a betadistribution:P (ν) = D(ν|α, β) :=1Zνα−1(1− ν)β−1 (6)where α, β > 0 are hyperparameters that determine the mean and the dispersion of thedistribution, while Z is a normalization factor. The resulting model is depicted on theright of Figure 1.First, we consider F (Q) as a functional of Q(ν) with Q(s) and Q(w) fixed. Introducingthe definitionsΨ =∑sN∑t=2Q(s)δ(st, st−1) α˜ = α+Ψ β˜ = β +N − 1−Ψ (7)and dropping terms that are independent of ν, this givesF [Q(ν)] =∫dνQ(ν)[lnQ(ν)−∑sQ(s)N∑t=2lnP (st|st−1, ν) − lnD(ν|α, β)]=∫dνQ(ν) [lnQ(ν)−Ψ lnν − (N − 1−Ψ) ln(1− ν)− lnD(ν|α, β)]=∫dνQ(ν) ln[Q(ν)D(ν|α˜, β˜)]+ C (8)This expression is minimized forQ(ν) = D(ν|α˜, β˜) (9)where D(.) denotes the beta distribution, defined in (6). Second, we have to optimizeQ(w) with Q(s) and Q(ν) fixed. At the current stage this has not been made explicit,though, and we rather assume a simple delta distribution located at the ML estimatewˆ: Q(w) = δ(w − wˆ). Finally, the optimal distribution over the hidden states st for4fixed Q(w) and Q(ν) is given, after some algebra2, by:Q(s) =1ZsP (s1)∏ta(st+1, st)∏tP (yt|st, wˆ) (10)in which Zs is a normalization constant, and where we have defined:ln a(st, st−1) = δ(st, st−1)Υ(α˜)−Υ(α˜+ β˜) + [1− δ(st, st−1)](Υ(β˜)− ln(K − 1))(11)Υ(x) =ddxln Γ(x) =ddxln∫∞0ux−1 exp(−u)du (12)Note that, after a sufficient number of iterations,Q(s) converges to the optimal posteriordistribution P (s|Y), and the most likely state sequence is given by maximizing eqn.(10).4 Simulation experiments4.1 A synthetic toy problemConsider a rogue casino in which a fair die is occasionally replaced by a loaded one.The latter has probability 0.5 of a six and probability 0.1 for the numbers one to five.The visitor to the casino observes a sequence of rolls, Y = (y1, . . . ,yN), but he doesnot know which rolls used a loaded (L) die and which used a fair (F) one. The natureof the die in the tth roll is represented by the hidden state st ∈ {L, F}, and theobjective is to infer the true sequence of states s = (s1, . . . , sN) from the sequence ofobservations. We parameterized the transition probabilities according to eqn.(5) witha value of ν = 0.9, and generated 20 training sets of length Ntrain = 50 from the truemodel. The emission parameters were kept at their correct values, while three differenttraining schemes were applied for adapting the transition parameter (always startingfrom a neutral value of ν = 0.5): (1) ML, (2) the Bayesian scheme with a fairly vagueprior on ν (α = 2.5, β = 1.5, corresponding to a mean of µ = 0.625 and a standarddeviation of σ = 0.22), (3) the Bayesian scheme using a more narrow prior with a higherprobability mass in the region around the correct value (α = 17, β = 3, corresponding toµ = 0.85, σ = 0.08). The prediction performance of the 20 ensuing models was assessedon an independent test set of length Ngen = 1000. The results are shown in Figure 2and suggest that both Bayesian training schemes lead to a significant improvement overresults obtained with ML.2See http://www.bioss.sari.ac.uk/∼dirk/ for a complete derivation.50 50 100020406080100Viterbi0 50 100020406080100Single Mode0 50 100020406080100Viterbi0 50 100020406080100Single Mode−1.79 −1.68−1.79−1.68Egen−1.79 −1.68−1.79−1.68EgenFigure 2: Comparison between training simulations with ML (horizontal axis) and thevariational Bayesian approach (vertical axis). Two different prior distributions werechosen. Top row: Informative prior, α = 17, β = 3. Bottom row: Vague prior,α = 2.5, β = 1.5. The methods were compared with three performance criteria. Leftcolumn: Percent correct classification of hidden states, as predicted by the mode of thejoint distribution P (s|Y) = P (s1, . . . , sN |Y) (Viterbi path). Middle column: Percentcorrect classification based on the single-site mode, P (st|yt). Right column: Generaliza-tion performance in terms of the normalized test-set log likelihood, Egen =1NlnP (Y).The crosses compare the performance of the ML (horizontal axis) and the Bayesianapproaches (vertical axis) for the various training sets. The diagonal line indicates anequal performance of the two methods. Crosses above this line represent simulations forwhich the Bayesian approach outperforms the ML method.4.2 Detection of recombination in DNA sequence alignmentThe second study is related to the bioinformatics problem of estimating phylogenetictrees from DNA sequence alignments. A phylogenetic tree is a directed acyclic graph,whose topology conveys information about evolutionary relatedness between differentspecies, whereas the edge lengths indicate evolutionary time in terms of the averagenumber of mutations per site. Standard training methods aim to optimize the topologyand the edge lengths by maximizing the likelihood of a given DNA sequence alignment[1], [3]. This approach, however, cannot deal with recombination, which is the exchangeof DNA subsequences between different strains or species and, as discussed in [8], cor-responds to a change of the tree topology in the affected region. We here follow [8] andmodel the process with an HMM, where each state st represents a different topology,670 75 80 85 90 95 100707580859095100SensitivityMLBayes90 92 94 96 98 10090919293949596979899100SpecificityMLBayesFigure 3: Comparison of the sensitivities (left) and specificities (right) obtained with theML (horizontal axis) and the Bayesian (vertical axis) training schemes. Each cross orcircle represents a different combination of training and test sets, where the difference inthe symbols is related to the different lengths of the recombinant regions in the respectivetest set (crosses: 200, 200; circles: 50, 250). The diagonal dashed line indicates an equalperformance of the ML and the Bayesian scheme, symbols above this line indicate aperformance improvement as a result of applying the Bayesian scheme (and vice versa).while the emission parameters are given by the branch lengths of the tree.In our experiment, we simulated the evolution and recombination processes in a 4-species tree according to the way described in [3] and [8]. All synthetic DNA sequencealignments contained 1000 sites and two recombinant zones, but varied in the edgelengths of the true tree and the locations and lengths of the recombinant regions. Eachof the data sets in turn was used for training an HMMwith either ML or according to theBayesian scheme described above [using two different settings of the hyperparameters:(α, β) = (10, 10), (150, 50)]. The prediction performance was subsequently assessed onthe remaining data sets not used for training, where we chose two performance measures.From the Viterbi path (that is, the mode of P (s|Y) = P (s1, . . . , sN |Y)) we determinedthe true prediction rate for detecting recombinant sites (the sensitivity) as well as thetrue prediction rate for non-recombinant sites (the specificity). A comparison of theperformance between training with ML and the Bayesian scheme is shown in Figure 3.33This comparison is based on the more informative prior, (α, β) = (150, 50). For the vaguer prior,(α, β) = (10, 10), the pattern in the plots is similar, but less pronounced.75 DiscussionWhile the Bayesian approach clearly outperformed ML on the synthetic toy problem,the results on the phylogenetic data sets are less clear-cut. The simulations suggest thatthe Bayesian approach improves the sensitivity for detecting recombinations at the costof a slight degradation in the specificity. Since the values for the latter are usually quitehigh (over 95%), the improvement in the former seems to be more relevant for practicalapplications. Further simulation studies are required to clarify this issue.References[1] P. Baldi and P. Brunak. Bioinformatics - The Machine Learning Approach. MITPress, 1988.[2] Y. Bengio and P. Frasconi. Input-Output HMMs for Sequence Processing. IEEETransactions on Neural Networks, 7(5), 1996.[3] R. Durbin, S. R. Eddy, A. Krogh, and Mitchison G. Biological sequence analysis.Probabilistic models of proteins and nucleic acids. Cambridge University Press,Combridge, UK, 1998.[4] W.R. Gilks, S. Richardson, and D.J. Spiegelhalter. Markov Chain Monte Carlo inPractice. Chapman and Hall, Suffolk, 1996. ISBN: 0-412-05551-1.[5] M. I. Jordan, Z. Ghahramani, T. S. Jaakola, and L. K. Saul. An Introduction toVariational Methods for Graphical Models. In M. I. Jordan, editor, Learning inGraphical Models, pages 105–161, Cambridge, MA, 1999. MIT Press.[6] D. J. C. MacKay. A practical Bayesian framework for backpropagation networks.Neural Computation, 4:448–472, 1992.[7] D.J.C. MacKay. Ensemble Learning for Hidden Markov Models. Technical report,Cavendish Laboratory, Cambridge CB3 0HE, UK, 1998.[8] G. McGuire. Statistical Methods for DNA Sequences: Detection of Recombinationand Distance Estimation. PhD thesis, University of Edinburgh, 1998.[9] R. M. Neal. Bayesian Learning for Neural Networks, volume 118 of Lecture Notesin Statistics. Springer, New York, 1996. ISBN 0-387-94724-8.[10] L.R. Rabiner. A tutorial on hidden Markov models and selected applications inspeech recognition. Proceedings of the IEEE, 77(2):257–286, 1989.8

http://eprints.gla.ac.uk/85684/1/85684.pdf

Bayesian regularization of hidden Markov models with an application to bioinformatics

Abstract

Similar works

Full text

Available Versions

Enlighten: Publications

Enlighten