IEICE TRANS. INF. & SYST., VOL.E85–D, NO.3 MARCH 2002

455

INVITED PAPER

Special Issue on the 2000 IEICE Excellent Paper Award

Multi-Space Probability Distribution HMM∗∗
Keiichi TOKUDA† , Takashi MASUKO†† , Noboru MIYAZAKI††∗ ,and Takao KOBAYASHI†† , Regular Members

SUMMARY
This paper proposes a new kind of hidden
Markov model (HMM) based on multi-space probability distribution, and derives a parameter estimation algorithm for theextended HMM. HMMs are widely used statistical models forcharacterizing sequences of speech spectra, and have been successfully applied to speech recognition systems. HMMs are categorized into discrete HMMs and continuous HMMs, which canmodel sequences of discrete symbols and continuous vectors, respectively. However, we cannot apply both the conventional discrete and continuous HMMs to observation sequences which consist of continuous values and discrete symbols: F0 pattern modeling of speech is a good illustration. The proposed HMM includes discrete HMM and continuous HMM as special cases, andfurthermore, can model sequences which consist of observation
vectors with variable dimensionality and discrete symbols.
key words: hidden Markov model, text-to-speech, F0, multispace probability distribution

1.

Introduction

The hidden Markov models (HMMs) are widely usedstatistical models, and have been successfully applied
to modeling sequences of speech spectra in speech
recognition systems.
The performance of HMMbased speech recognition systems has been improvedby techniques which utilize the ﬂexibility of HMMs:
context-dependent modeling [1], dynamic feature parameters [2], mixtures of Gaussian densities [3], tying
techniques (e.g., [4]), and speaker/environment adaptation techniques (e.g., [5]).
HMMs are categorized into discrete and continuous
HMMs, which can model sequences of discrete symbols
and continuous vectors, respectively. However, we cannot apply both the conventional discrete and continuous HMMs to observations which consist of continuous
values and discrete symbols. Modeling the fundumental
frequency (F0) pattern of speech is a good illustration:
we cannot directly apply both the conventional discrete
Manuscript received November 6, 2001.
The author is with the Department of Computer Science, Nagoya Institute of Technology, Nagoya-shi, 466–8555
Japan.
††
The authors are with the Interdisciplinary Graduate
School of Science and Engineering, Tokyo Institute of Technology, Yokohama-shi, 226–8502 Japan.
∗
Presently, with the NTT Communication Science Laboratories.
∗∗
This article was originally published in IEICE Trans.
Inf. & Syst. (Japanese Edition) vol.J83-D-II, no.7, pp.1579–
1589, July 2000.
†

and continuous HMMs to F0 pattern modeling since
F0 values are not deﬁned in the unvoiced region, i.e.,the observation sequence of an F0 pattern is composedof one-dimensional continuous values and discrete symbols which represent “unvoiced.” Several methods havebeen investigated for handling the unvoiced region: i)
replacing each “unvoiced” symbol by a random vector generated from a probability density function (pdf)with a large variance and then modeling the random
vectors explicitly in the continuous HMMs [6], ii) modeling the “unvoiced” symbols explicitly in the continuous HMMs by replacing each “unvoiced” symbol with 0
and adding an extra pdf for the “unvoiced” symbol toeach mixture [7], iii) assuming that F0 values always exist but they cannot be observed in the unvoiced region
and applying the EM algorithm [8]. Although approachiii) is appropriate from the viewpoint of statistical modeling, it intends to estimate F0 values which are not existent contradictorily. Approaches i) and ii) are basedon heuristic assumptions. As a result, we cannot derive
statistical techniques, e.g., context-dependent modeling, speaker/environment adaptation techniques, in astatistically correct manner.
This paper describes a new kind of HMM in which
the state output probabilities are deﬁned by multispace probability distributions, and derives its reestimation formulas. Each space in the multi-space probability distribution has its weight and a continuous pdf
whose dimensionality depends on the space. An observation consists of an n-dimensional continuous vector
and a set of space indices which specify n-dimensional
spaces. As a result, the extended HMM includes both
the discrete and continuous mixture HMMs as special
cases, and furthermore, can model the sequences of observation vectors with variable dimensionality including
zero-dimensional observations, i.e., discrete symbols.
This paper is organized as follows. Multi-spaceprobability distribution and multi-space probability
distribution HMM (MSD-HMM) are deﬁned in Sects. 2
and 3, respectively. A reestimation algorithm for MSDHMMs is derived in Sect. 4. The relation between theconventional and the proposed HMMs, and the application of MSD-HMM to F0 pattern modeling, are discussed in Sect. 5. Concluding remarks and our plans
for future work are presented in the ﬁnal section.

IEICE TRANS. INF. & SYST., VOL.E85–D, NO.3 MARCH 2002

456

2.

Multi-Space Probability Distribution

We consider a sample space Ω shown in Fig. 1, which
consists of G spaces:
Ω=

G


Ωg ,

(1)

g=1

where Ωg is an ng -dimensional real space Rng , speciﬁed by space index g. While each space has its own
dimensionality, some of them may have the same dimensionality.
Each space Ωg has its probability wg , i.e., P (Ωg ) =Gwg , where g=1 wg = 1. If ng > 0, each space has a
pdf function Ng (x), x ∈ Rng , where Ng (x)dx = 1.
We assume that Ωg contains only one sample point if
ng = 0. Accordingly, we have P (Ω) = 1.
Each event E, which will be considered in this paper, is represented by a random vector o which consistsof a set of space indices X and a continuous randomvariable x ∈ Rn , that is,o = (X, x),(2)where all spaces speciﬁed by X are n-dimensional. On
the other hand, X does not necessarily include all indices which specify n-dimensional spaces (see o1 ando2 in Fig. 1). It is noted that not only the observation
vector x but also the space index set X is a randomvariable, which is determined by an observation device(or feature extractor) at each observation. The observation probability of o is deﬁned bywg Ng (V (o)),(3)b(o) =g∈S(o)

where
S(o) = X,

V (o) = x.

(4)

It is noted that, although Ng (x) does not exist for ng =0 since Ωg contains only one sample point, for simplicity
of notation, we deﬁne Ng (x) ≡ 1 for ng = 0.
Some examples of observations are shown in Fig. 1.
An observation o1 consists of a set of space indices
X1 = {1, 2, G} and a three-dimensional vector x1 ∈R3 . Thus the random variable x is drawn from one ofthree spaces Ω1 , Ω2 , ΩG ∈ R3 , and its pdf is given by
w1 N1 (x) + w2 N2 (x) + wG NG (x).
The probability distribution deﬁned above, which
will be referred to as multi-space probability distribution (MSD) in this paper, is the same as the discrete
distribution when ng ≡ 0. Furthermore, if ng ≡ m > 0
and S(o) ≡ {1, 2, . . . , G}, the multi-space probability
distribution is represented by a G-mixture pdf. Thus
the multi-space probability distribution is more generalthan either discrete or continuous distributions.
The following example shows that the multi-space

Fig. 1

Multi-space probability distribution and observations.

probability distribution conforms to statistical phenomena in the real world:
A man is ﬁshing in a pond. There are red ﬁshes,
blue ﬁshes, and tortoises in the pond. In addition, some junk articles are in the pond. Whenhe catches a ﬁsh, he is interested in the kind of
the ﬁsh and its size, for example, the length andheight. When he catches a tortoise, it is suﬃcient
to measure the diameter if we assume that the tortoise has a circular shape. Furthermore, when hecatches a junk article, he takes no interest in its
size.
In this case, the sample space consists of four spaces:Ω1 : two-dimensional space correponding to lengths andheights of red ﬁshes.
Ω2 : two-dimensional space correponding to lengths andheights of blue ﬁshes.
Ω3 : one-dimensional space correponding to diametersof tortoises.
Ω4 : zero-dimensional space correponding to junk articles.
The weights w1 , w2 , w3 , w4 are determined by the ratio of red ﬁshes, blue ﬁshes, tortoises, and junk articles in the pond. Functions N1 (·) and N2 (·) are twodimensional pdfs for sizes (lengths and heights) of red
ﬁshes and blue ﬁshes, respectively. The function N3 (·)
is the one-dimensional pdf for diameters of tortoises.
For example, when the man catches a red ﬁsh, theobservation is given by o = ({1}, x), where x is a
two-dimensional vector which represents the length andheight of the red ﬁsh. Suppose that he is ﬁshing day
and night, and during the night, he cannot distinguishbetween the colors of ﬁshes, while he can measure theirlengths and heights. In this case, the observation of aﬁsh at night is given by o = ({1, 2}, x).

TOKUDA et al.: MULTI-SPACE PROBABILITY DISTRIBUTION HMM

457




···





aq1 q2 wq1 g Nq2 g (V (o2 ))

g∈S(o2 )aqT −1 qT wqT g NqT g (V (oT ))

g∈S(oT )=T 

aqt−1 qt wqt lt Nqt lt (V (ot )),(6)

all q,l t=1

where q = {q1 , q2 , . . . , qT } is a possible state sequence,l = {l1 , l2 , . . . , lT } ∈ {S(o1 ) × S(o2 ) × . . . × S(oT )}is a sequence of space indices which is possible for theobservation sequence O, and aq0 j denotes πj .
Equation (6) can be calculated eﬃciently throughthe forward and backward variables:
Fig. 23.

αt (i) = P (o1 , o2 , . . . , ot , qt = i|λ)βt (i) = P (ot+1 , ot+2 , . . . , oT |qt = i, λ),An HMM based on multi-space probability distribution.

HMMs Based on Multi-Space Probability
Distribution

By using the multi-space distribution, we deﬁne a newkind of HMM. In this paper, we call it multi-spaceprobability distribution HMM (MSD-HMM). The output probability in each state of MSD-HMM is givenby the multi-space probability distribution deﬁned in
the previous section. An N -state MSD-HMM λ is
speciﬁed by the initial state probability distribution
π = {πj }N
j=1 , the state transition probability distribution A = {aij }N
i, j=1 , and the state output probability
distribution B = {bi (·)}N
i=1 , wherebi (o) =wig Nig (V (o)).

P (O|λ)=aqt−1 qt bqt (ot )Taqt−1 qt

all q t=1

=


all q





g∈S(ot )g∈S(o1 )

α1 (i) = πi bi (o1 ), 1 ≤ i ≤ N
1≤i≤N
βT (i) = 1,

(9)

2. Recursion:Nαt+1 (i) = αt (j)aji  bi (ot+1 ),j=1

1 ≤ i ≤ N,
βt (i) =N

t = 1, 2, . . . , T − 1(10)

aij bj (ot+1 )βt+1 (j),j=1

1 ≤ i ≤ N,

P (O|λ) =wqt g Nqt g (V (ot ))

aq0 q1 wq1 g Nq1 g (V (o1 ))

N

i=1

t = T − 1, 2, . . . , 1.

(11)

αT (i) =N

aq0 i bi (o1 )β1 (i).

(12)

i=1

The forward and backward variables are also used for
calculating the reestimation formulas derived in thenext section (i.e., calculation of Eqs. (15) and (16)).
4.

all q t=1=

1. Initialization:

According to the deﬁnitions, (6) can be calculated as

As shown in Fig. 2, each state i has G pdfs Ni1 (·),
Ni2 (·), . . . , NiG (·), and their weights wi1 , wi2 , . . . , wiG ,where G
g=1 wig = 1. The observation probability ofO = {o1 , o2 , . . . , oT } can be written asT

which can be calculated with the forward-backward inductive procedure in a manner similar to conventional
HMMs:

(5)

g∈S(o)(7)(8)

Reestimation Algorithm

For a given observation sequence O and a particular
choice of MSD-HMM, the objective in maximum likelihood estimation is to maximize the observation likelihood P (O|λ) given by Eq. (6), over all parameters in
λ. In a manner similar to those reported in [9] and[3], we derive reestimation formulas for the maximum
likelihood estimation of MSD-HMM.

IEICE TRANS. INF. & SYST., VOL.E85–D, NO.3 MARCH 2002

458

P (qt = i, O|λ)P (lt = h|qt = i, O, λ)P (O|λ)αt (i)βt (i)wih Nih (V (ot ))
= N
· 
.

wig Nig (V (ot ))αt (j)βt (j) g∈S(o )=4.1 Q-Function
An auxiliary function Q(λ , λ) of current parameters λand new parameters λ is deﬁned as follows:Q(λ , λ) =P (O, q, l|λ ) log P (O, q, l|λ). (13)

t

j=1

(15)

all q,lIn the following, we assume Nig (·) to be the Gaussian
density with mean vector µig and covariance matrix
Σig . However, extension to elliptically symmetric densities which satisfy the consistency conditions of Kolmogorov is straightforward. We present the following
three theorems:

Similarly, the posterior probability of transitions fromstate i to state j at time t + 1 is given by
ξt (i, j) = P (qt = i, qt+1 = j|O, λ)P (qt = i, qt+1 = j, O|λ)=P (O|λ)αt (i)aij bj (ot+1 )βt+1 (j)
.
= N N

αt (m)amk bk (ot+1 )βt+1 (k)Theorem 1:Q(λ , λ) ≥ Q(λ , λ ) → P (O, λ) ≥ P (O, λ )(14)Theorem 2: If, for each space Ωg , there are among
V (o1 ), V (o2 ), . . . , V (oT ), ng + 1 observations g ∈S(ot ), any ng of which are linearly independent,Q(λ , λ) has a unique global maximum as a functionof λ, and this maximum is the one and only criticalpoint.
Theorem 3: A parameter set λ is a critical point of
the likelihood P (O|λ) if and only if it is a critical pointof the Q-function.
Theorems 1 and 3 can be proved in a similar manner tothe conventional HMM. We have to newly prove Theorem 2, which conﬁrms that the Q-function has a uniqueglobal maximum as a function of λ because the proposed HMM has a diﬀerent state output probability
distribution from the conventional descrete or continuous HMMs. The proof of Theorem 2 is given in Appendix.
We deﬁne the parameter reestimates to be thosewhich maximize Q(λ , λ) as a function of λ, λ being thelatest estimates. Because of the above theorems, thesequence of reestimates obtained in this way produces amonotonic increase in the likelihood unless λ is a criticalpoint of the likelihood.
4.2 Maximization of Q-Function
For given observation sequence O and model λ , wederive parameters of λ which maximize Q(λ , λ).
The posterior probability of being in state i at timet, given the observation sequence O and model λ, isgiven byγt (i, h)= P (qt = i, lt = h|O, λ)= P (qt = i|O, λ)P (lt = h|qt = i, O, λ)m=1 k=1(16)
We deﬁne a function T (O, g) which returns a set oftime t at which the space index set S(ot ) includes spaceindex g:T (O, g) = {t|g ∈ S(ot )}.

(17)

By introducing this fucntion, the following manipulations of the equations can be carried out in a similarmanner to the conventional continuous mixture HMMs.
From Eq. (6), log P (O, q, l|λ) can be written aslog P (O, q, l|λ)=Tlog aqt−1 qt + log wqt lt + log Nqt lt (V (ot )) . (18)t=1Hence, Q-function Eq. (13) can be written as
Q(λ , λ)P (O, q, l|λ ) log aq0 q1=
all q,l+

P (O, q, l|λ )+

P (O, q, l|λ )all q,lT

log wqt lt

t=1

all q,l+log aqt qt+1

t=1

all q,lT−1P (O, q, l|λ )Tlog Nqt lt (V (ot )).

t=1

(19)
The ﬁrst term of Eq. (19), which is related to aq0 q1 , i.e.,
πq1 , is given byTOKUDA et al.: MULTI-SPACE PROBABILITY DISTRIBUTION HMM

459



N
under the constraint
i=1 yi = 1, yi ≥ 0. Therefore, the parameters πi , aij , and wig which maximize Eq. (20), subject to the stochastic constraints
NG
N
i=1 πi = 1,
j=1 aij = 1, and
g=1 wg = 1, respectively, can be derived as

P (O, q, l|λ ) log aq0 q1

all q,l=N P (O, q1 = i, l|λ ) log aq0 i

i=1 all l=N

P (O, q1 = i|λ ) log πi .

(20)

i=1

The second term of Eq. (19), which is related to aij , isgiven by

P (O, q, l|λ )log aqt qt+1

N 
N T
−1 



P (O, q1 = i|λ )P (O, q1 = i|λ )=N
P (O|λ )
P (O, q1 = j|λ )j=1= P (q1 = i|O, λ )
=γ1 (i, g)

P (O, qt = i, qt+1 = j, l|λ ) log aij

T−1i=1 j=1 t=1 all l

=

N T
−1



P (O, qt = i, qt+1 = j|λ ) log aij .

aij =i,j=1 t=1

The third term of Eq. (19), which is related to wig , isgiven by

P (O, q, l|λ )

=

N 
T



N 
G


= t=1 T −1

P (O, qt = i, lt = g|λ ) log wig .

The fourth term of Eq. (19), which is related to Nig (·),
is given by

N 
T





wig =

· log Nig (V (ot ))=P (O, qt = i, lt = h|λ )h=1 t∈T (O,h)



i=1 g=1 t∈T (O,g)

(23)

N
Equations (20)–(22) have the form of
i=1 ui log yi ,which attains its unique maximum pointui(24)
yi = Nuj
j=1

P (O, qt = i, lt = g|λ )t∈T (O,g)G
P (O, qt = i, lt = g|λ )· log Nig (V (ot )).

γt (i, g)

t=1 g∈S(ot )P (O, qt = i, lt = g|λ )



(26)

log Nqt lt (V (ot ))i=1 t=1 g∈S(ot )

N 
G


ξt (i, j)= T −1 t=1
 

t=1

all q,l=T

P (qt = i|O, λ )

t=1
T
−1


(22)

P (O, q, l|λ )P (qt = i, qt+1 = j|O, λ )= t=1 T −1

i=1 g=1 t∈T (O,g)

P (O, qt = i|λ )

t=1
T
−1


P (O, qt = i, lt = g|λ ) log wig



P (O, qt = i, qt+1 = k|λ )P (O, qt = i, qt+1 = j|λ )log wqt lt

i=1 t=1 g∈S(ot )

=

t=1
N T
−1



T
−1


t=1

all q,lP (O, qt = i, qt+1 = j|λ )k=1 t=1

(21)

T


(25)

g∈S(o1 )t=1

all q,l

=

T
−1


πi =

=

γt (i, g)t∈T (O,g)
G




.

(27)

γt (i, h)h=1 t∈T (O,h)

When Nig (·), ng > 0 is the ng -dimensional Gaussian
density function with mean vector µig and covariancematrix Σig , Eq. (23) is maximized by setting the partial
derivatives with respect to µig and Σig −1 :

IEICE TRANS. INF. & SYST., VOL.E85–D, NO.3 MARCH 2002

460

∂
∂µig





P (O, qt = i, lt = g|λ )t∈T (O,g)· log Nig (V (ot )) = 0
∂P (O, qt = i, lt = g|λ )−1∂Σig t∈T (O,g)· log Nig (V (ot )) = 0.

t∈T (O,g)(28)(V (ot ) − µig )(V (ot ) − µig )T


Σig =
(29)

P (O, qt = i, lt = g|λ )t∈T (O,g)Fromγt (i, g)t∈T (O,g)

∂
∂µig
=



(V (ot ) − µig )(V (ot ) − µig )TP (O, qt = i, lt = g|λ ) log Nig (V (ot ))=t∈T (O,g)

P (O, qt = i, lt = g|λ )t∈T (O,g)

(30)

and

∂
P (O, qt = i, lt = g|λ ) log Nig (V (ot ))−1∂Σig t∈T (O,g)
P (O, qt = i, lt = g|λ )=t∈T (O,g)∂log Nig (V (ot ))
∂Σ−1
ig
1 
=
P (O, qt = i, lt = g|λ )2·t∈T (O,g)(31)µig and Σig are given byµig =P (O, qt = i, lt = g|λ )V (ot )t∈T (O,g)

P (O, qt = i, lt = g|λ )t∈T (O,g)=γt (i, g)V (ot )t∈T (O,g)t∈T (O,g)γt (i, g),

ng > 0

,

(32)

(33)

respectively. From the condition mentioned in Theorem
2, it can be shown that each Σig is positive deﬁnite.
From Sect. 4.1, by iterating the following procedure: 1) calculating λ which maximizes Q(λ , λ) by
Eqs. (25)–(27), (32), and (33), and 2) substituting theobtained λ for λ , we can obtain a critical point of
P (O|λ).
5.

· Σig − (V (ot ) − µig )(V (ot ) − µig )T= 0,

γt (i, g)

ng > 0,∂log Nig (V (ot ))
·
∂µig

=
P (O, qt = i, lt = g|λ )·Σ−1ig (V (ot ) − µig )=0t∈T (O,g)t∈T (O,g)and

P (O, qt = i, lt = g|λ )Discussion

5.1 Relation to Discrete Distribution HMM and Continuous Distribution HMM
The MSD-HMM includes the discrete HMM and thecontinuous mixture HMM as special cases since themulti-space probability distribution includes the discrete distribution and the continuous distribution. If
ng ≡ 0, the MSD-HMM is the same as the discrete
HMM. In the case where S(ot ) speciﬁes one space, i.e.,|S(ot )| ≡ 1, the MSD-HMM is exactly the same as
the conventional discrete HMM. If |S(ot )| ≥ 1, theMSD-HMM is the same as the discrete HMM basedon the multi-labeling VQ [10]. If ng ≡ m > 0 andS(o) ≡ {1, 2, . . . , G}, the MSD-HMM is the same as
the continuous G-mixture HMM. These can also be conﬁrmed by the fact that if ng ≡ 0 and |S(ot )| ≡ 1,
the reestimation formulas Eqs. (25)–(27) are the sameas those for discrete HMM of codebook size G, andif ng ≡ m and S(ot ) ≡ {1, 2, . . . , G}, the reestimation formulas Eqs. (25)–(33) are the same as thosefor continuous HMM with m-dimensional G-mixture
densities. Accordingly, MSD-HMM includes the discrete and continuous mixture HMMs as special cases,
and furthermore, can model the sequence of observation vectors with variable dimensionality including
zero-dimensional observations, i.e., discrete symbols.
In addition, multi-channel HMMs [11] are also related to MSD-HMMs. Multi-channel HMMs have a
special structure similar to MSD-HMMs. However,TOKUDA et al.: MULTI-SPACE PROBABILITY DISTRIBUTION HMM

461

they assume that each channel always observes a discrete symbol, and they cannot be applied to the observation sequence composed of continuous vectors withvariable dimensionality including zero-dimensional observations, i.e., discrete symbols. On the other hand,
MSD-HMM includes the multi-channel HMM which
was ﬁnally derived in [11] as a special case under thefollowing conditions:• The sample space consists of zero-dimensional
spaces, each of which has a one-to-one correspondence with each symbol used in the multi-channel
HMM.
• The observation consists of M space indices, eachof which has a one-to-one correspondence with achannel and is drawn from symbols used in thechannel.

(a) natural F0(b) generated F0
Fig. 3 Generated F0 pattern for the sentence “heikiNbairitsuwo sageta keisekiga aru.”

5.2 Application to F0 Pattern Modeling
While the observation of F0 has a continuous value inthe voiced region, there exist no values for the unvoiced
region. We can model this kind of observation sequence
assuming that the observed F0 value occurs from onedimensional spaces and the “unvoiced” symbol occurs
from the zero-dimensional space deﬁned in Sect. 2, thatis, by setting ng = 1 (g = 1, 2, . . . , G − 1), nG = 0and
{1, 2, . . . , G − 1},(voiced)(34)
S(ot ) ={G},(unvoiced)
the MSD-HMM can cope with F0 patterns including
the unvoiced region without heuristic assumptions. In
this case, the observed F0 value is assumed to be drawn
from a continuous (G − 1)-mixture pdf.
Experiments reported in [12] have shown that thelikelihood of the MSD-HMM for the training data increases monotonically by calculating the reestimation
formulas iteratively. From the trained MSD-HMMs,
we can generate F0 patterns which approximate thoseof natural speech by using an algorithm (described in[13] as the “case 1” algorithm) for speech parameter
generation from HMMs with dynamic features . An
example is shown in Fig. 3, without an explanation of
the experimental conditions because of limitations ofspace† .
Real world phenomena of time sequences are notnecessarily observed as a sequence of discrete symbols or continuous vectors. Accordingly, the proposed
HMMs can be applied to not only F0 pattern modeling but also the modeling of various kinds of time sequences which consist of continuous vectors with variable dimensionality including zero-dimensional vectors,i.e., discrete symbols. As a result, MSD-HMM expected
to be useful in such research areas as the prediction ofhuman actions and economic forecasting.
Algorithms based on a statistical framework, whichhave been developed for the conventional HMMs, canbe applied or extended to MSD-HMMs since the proposed scheme does not utilize any heuristic assumption
or approximation. Actually, in the example shown in
Fig. 3, we derived a model clustering scheme based on
the minimum description length (MDL) principle, andused it in the model training.
6.

Conclusion

A multi-space probability distribution HMM has been
proposed and its reestimation formulas are derived.
The MSD-HMM includes the discrete HMM and thecontinuous mixture HMM as special cases, and furthermore, can cope with the sequence of observation
vectors with variable dimensionality including zerodimensional observations, i.e., discrete symbols. As
a result, MSD-HMMs can model F0 patterns without
heuristic assumptions.
In the near future, we will present a speech synthesis system in which sequences of speech spectra [15], F0
patterns [12] and state durations [16] are modeled by
MSD-HMM in a uniﬁed framework [14]. Fundamental
frequency (F0) pattern modeling based on MSD-HMM
may also be useful for enhancing the performance ofspeech recognition systems.
Acknowledgment
This work was partially supported by the Ministry of
Education, Culture, Sports, Science and Technology,
Grant-in-Aid for Scientiﬁc Research (B), 10555125, andthe Research Foundation for the Electrotechnology of
Chubu, Grant for Research. The authors are gratefulto Mr. Takayoshi Yoshimura, PhD student, Nagoya Institute of Technology, for carrying out experiments.
†

Detailed experimental conditions can be found in [14].

IEICE TRANS. INF. & SYST., VOL.E85–D, NO.3 MARCH 2002

462

References[1] S. Schwartz, Y-L. Chow, O. Kimball, S. Roucos,
M. Krasner, and J. Makhoul, “Context-dependent modeling for acoustic-phonetic of continuous speech,” Proc.
ICASSP85, pp.1205–1208, March 1985.
[2] S. Furui, “Speaker independent isolated word rcognition
using dynamic features of speech spectrum,” IEEE Trans.
Acoust., Speech & Signal Process., vol.34, no.1, pp.52–59,
Feb. 1986.
[3] B.-H. Juang, “Maximum-likelihood estimation for mixture
multivariate stochastic observations of Markov chains,”
AT&T Technical Journal, vol.64, no.6, pp.1235–1249, 1985.
[4] J.J. Odell, “The use of context in large vocabulary speech
recognition,” PhD Thesis, Cambridge University, March
1995.
[5] C.-H. Lee, C.H. Lin, and B.H. Juang, “A study on speaker
adaptation of the parameters of continuous density hidden
Markov models,” IEEE Trans. Acoust., Speech & Signal
Process., vol.39, no.4, pp.806–814, April 1992.
[6] G.J. Freij and F. Fallside, “Lexical stress recognition using hidden Markov models,” Proc. ICASSP88, pp.135–138,
April 1988.
[7] U. Jensen, R.K. Moore, P. Dalsgaard, and B. Lindberg,“Modeling intonation contours at the phrase level using continuous density hidden Markov models,” Computer
Speech and Language, vol.8, no.3, pp.247–260, July 1994.
[8] K. Ross and M. Ostendorf, “A dynamical system model for
generating F0 for synthesis,” Proc. ESCA/IEEE Workshop
on Speech Synthesis, pp.131–134, Sept. 1994.
[9] L.A. Liporace, “Maximum likelihood estimation for multivariate observations of Markov sources,” IEEE Trans. Inf.
Theory, vol.28, no.5, pp.729–734, Sept. 1982.
[10] M. Nishimura and K. Toshioka, “HMM-based speech
recognition using multi-dimensional multi-labeling,” Proc.
ICASSP87, pp.1163–1166, May 1987.
[11] D. Xu, C. Fancourt, and C.Wang, “Multi-channel HMM,”
Proc. ICASSP96, pp.841–844, May 1996.
[12] N. Miyazaki, K. Tokuda, T. Masuko, and T. Kobayashi,“A study on pitch pattern generation using HMMs basedon multi-space probability distributions,” IEICE Technical
Report, SP98-12, 1998.
[13] K. Tokuda, T. Yoshimura, T. Masuko, T. Kobayashi, andT. Kitamura, “Speech parameter generation algorithms
for HMM-based speech synthesis,” ICASSP2000, pp.1315–
1318, June 2000.
[14] T. Yoshimura, K. Tokuda, T. Masuko, T. Kobayashi, andT. Kitamura, “Simultaneous modeling of spectrum, pitchand duration in HMM-based speech synthesis,” Proc. EUROSPEECH99, pp.2347–2350, Sept. 1999.
[15] T. Masuko, K. Tokuda, T. Kobayashi, and S. Imai, “Speechsynthesis from HMMs using dynamic features,” Proc.
ICASSP96, pp.389–392, May 1996.
[16] T. Yoshimura, K. Tokuda, T. Masuko, T. Kobayashi, andT. Kitamura, “Duration modeling for HMM-based speechsynthesis,” Proc. ICSLP98, pp.29–32, Nov.–Dec. 1998.

Appendix: Proof of Theorem 2
The proof consists of the following three parts:a) The second derivative of the Q-function along any
direction in the parameter space is strictly negativeat a critical point. This implies that any criticalpoint is a relative maximum and that if there are
more than one they are isolated.
b) Q(λ , λ) → −∞ as λ approaches the ﬁnite boundary of the parameter space or the point at ∞. This
property implies that the global maximum is a critical point.
c) The critical point is unique.
Proof (a)From Sect. 4.2, the Q-function can be written as
Q(λ , λ)=NP (O, q1 = i|λ ) log πi

i=1

+

N T
−1



P (O, qt = i, qt+1 = j|λ ) log aij

i,j=1 t=1

+

N 
G




P (O, qt = i, lt = g|λ )i=1 g=1 t∈T (O,g)

1nglog(2π) + log |C ig |· log wig −22
1− (V (ot ) − µig )T C ig (V (ot ) − µig ) ,2(A· 1)
where C ig = Σ−1ig . From the condition on observations
ot , described in Theorem 2, Σig and Σ−1
ig are positivedeﬁnite if Σig is calculated by Eq. (33).
Let us express λ as a linear combination of two
arbitrary points: λ = θλ(1) + (1 − θ)λ(2) , where 0 ≤ θ ≤1. That is,(1)

πi = θπi(2)+ (1 − θ)πi

(A· 2)(1)(2)
aij = θaij + (1 − θ)aij(A· 3)(1)(2)
wig = θwig + (1 − θ)wig(A· 4)(1)(2)
C ig = θC ig + (1 − θ)C ig(A· 5)(1)(2)µig = θµig + (1 − θ)µig .

(A· 6)

Substituting these equations for Eq. (A· 1) and taking
the second derivative with respect to θ, we obtain
∂2Q
∂θ2
=

N


P (O, q1 = i|λ )

i=1

+

N

i,j=1

(1)

−(πi
(1)

(θπi

(2)

− πi )2(2)+ (1 − θ)πi )2

P (O, qt = i, qt+1 = j|λ )TOKUDA et al.: MULTI-SPACE PROBABILITY DISTRIBUTION HMM

463
(1)

·

(2)

−(aij − aij )2

+

(1)
(2)
(θaij + (1 − θ)aij )2

+

N 
G




P (O, qt = i, lt = g|λ )+· log wig −
+

(1)

(2)

+2(µig − µig )T (C ig − C ig )
(1)

(2)

·[V (ot ) − (θµig + (1 − θ)µig )] ,
(1)

(2)

(1)

(A· 7)(2)

where xigk and xigk satisfy xigk = θxigk + (1 − θ)xigk
for xigk which are the diagonal entries of U ig C ig U −1
ig ,and the orthogonal matrix U ig diagonalizes C ig .
At a critical point, from the relation

∂Q 

∂µig 
(1)(2)µig = θµig + (1 − θ)µig
(1)

(2)

= (θC ig + (1 − θ)C ig )·P (O, qt = i, lt = g|λ )

ng
ng
1
1
2
log yigs −yigs ztigs
,
2 s=1
2 s=1

(A· 10)1) πi → 02) aij → 023) ztigs→∞4) wig → 05) yigs → 06) yigs → ∞When one of the conditions 1)–5) holds, it is obvious that Q(λ , λ) → −∞ because one of the terms
in Eq. (A· 10) approaches −∞. In the case where
yigs → ∞, from the condition on observations ot , de2
has a nonzero positive value
scribed in Theorem 2, ztigs
at some t. Thus,(A· 11)

As a result, Q(λ , λ) → −∞, as λ approaches the ﬁnite
boundary of the parameter space or the point at ∞.



(1)
(2)
· V (ot ) − (θµig + (1 − θ)µig )= 0,ng
log(2π)22log yigs − ztigs
yigs → −∞.

P (O, qt = i, lt = g|λ )t∈T (O,g)

Proof (c)(A· 8)

the sum involving the term bracketed by [ ] in Eq. (A· 7)
vanishes. All of the remaining terms have negative values. Therefore, independent of λ(1) and λ(2) ,From Proof (a), if there are multiple critical points,
they are isolated. Assume that C ig = τ Tig τ ig , whereτ ig ’s are triangular and positive deﬁnite. We canrewrite Eq. (A· 1) as
Q(λ , λ)

∂2Q
≤0
∂θ2

(A· 9)=

along any direction.

N

P (O, q1 = i|λ ) log πii=1+Proof (b)

N T
−1



P (O, qt = i, qt+1 = j|λ ) log aij

i,j=1 t=1

The Q-function Q(λ , λ) can be rewritten as

N

i=1

+

N 
G




P (O, qt = i, lt = g|λ )i=1 g=1 t∈T (O,g)

Q(λ , λ)=where yigs , s = 1, 2, . . . , ng are are eigenvalues of C ig ,eigs , s = 1, 2, . . . , ng are an orthonormal set of eigenvectors of C ig , and ztigs = (V (ot ) − µig )T eigs .
When λ approaches ∞ or the boundary of the parameter space, one of the following conditions holds.

(2)−(µig − µig )T


(1)
(2)
(1)
(2)
· θC ig + (1 − θ)C ig (µig − µig )
(2)

N 
G


i=1 g=1 t∈T (O,g)(1)(2)−(wig − wig )2
·
(1)
(2)
(θwig + (1 − θ)wig )2
ng
(1)
(2)
−(xigk − xigk )2
1
+
(1)
(2) 2
2
k=1 (θxigk + (1 − θ)xigk )(1)

P (O, qt = i, qt+1 = j|λ ) log aij

i,j=1 t=1

i=1 g=1 t∈T (O,g)

(1)

N T
−1



P (O, q1 = i|λ ) log πinglog(2π) + log |τ ig |· log wig −2
1− ||τ ig (V (ot ) − µig )||2 .
2(A· 12)

IEICE TRANS. INF. & SYST., VOL.E85–D, NO.3 MARCH 2002

464

The change of variables {πi , aij , wig , µig , C ig } → {πi ,aij , wig , µig , τ ig }, which is a diﬀeomorphism, maps
critical points onto critical points. Therefore, the global
maximum is the unique ciritical point since Eq. (A· 12)
is convex with respect to πi , aij , wig , µig , τ ig .

Keiichi Tokuda
received the B.E. degree in electrical and electronic engineering from the Nagoya Institute of Technology, Nagoya, Japan, in 1984, and the M.E.
and Dr.Eng. degrees in information processing from the Tokyo Institute of Technology, Tokyo, Japan, in 1986 and 1989,
respectively. From 1989 to 1996 he was a
Research Associate at the Department ofElectronic and Electric Engineering, Tokyo Institute of Technology. Since 1996 hehas been with the Department of Computer Science, Nagoya Institute of Technology as Associate Professor. He is a co-recipientof both the Best Paper Award and the Inose Award from the
IEICE in 2001, and the TELECOM System Technology Prize
from the Telecommunications Advancement Foundation Award,
Japan, in 2001. His research interests include speech coding,speech synthesis and recognition and multimodal signal processing. He is a member of IEEE, ISCA, IPSJ, ASJ and JSAI.

Takashi Masuko
received the B.E.
degree in computer science, and the M.E.
degree in intelligence science from Tokyo
Institute of Technology, Tokyo, Japan, in
1993 and 1995, respectively. In 1995, hejoined the Precision and Intelligence Laboratory, Tokyo Institute of Technology, asa Research Associate. He is currently
a Research Associate at the Interdisciplinary Graduate School of Science andEngineering, Tokyo Institute of Technology, Yokohama, Japan. He is a co-recipient of both the Best
Paper Award and the Inose Award from the IEICE in 2001. His
research interests include speech synthesis, speech recognition,
speech coding, and multimodal interface. He is a member ofIEEE, ISCA and ASJ.

Noboru Miyazaki
was born inEhime, Japan, in 1972. He received the
B.E. degree in information engineeringand the M.E. degree in intelligence science from Tokyo Institute of Technology, Tokyo, Japan, in 1995 and 1997, respectively. In 1997, he joined the Basic Research Laboratories, Nippon Telegraph and Telephone Corporation (NTT),
Japan as a researcher. He is currently
a researcher at the NTT Communication
Science Laboratories. He is a co-recipient of both the Paper
Award and the Inose Award from the IEICE in 2001. His research interests include speech recognition and spoken dialogue
processing. He is a member of ASJ.

Takao Kobayashi
received the B.E.
degree in electrical engineering, the M.E.
and Dr.Eng. degrees in information processing from Tokyo Institute of Technology, Tokyo, Japan, in 1977, 1979, and1982, respectively. In 1982, he joined
the Research Laboratory of Precision Machinery and Electronics, Tokyo Instituteof Technology as a Research Associate.
He became an Associate Professor at thesame Laboratory in 1989. He is currently
a Professor at the Interdisciplinary Graduate School of Science and Engineering, Tokyo Institute of Technology, Yokohama,Japan. He is a co-recipient of both the Best Paper Award andthe Inose Award from the IEICE in 2001, and the TELECOM
System Technology Prize from the Telecommunications Advancement Foundation Award, Japan, in 2001. His research interests include speech analysis and synthesis, speech coding, speech
recognition, and multimodal interface. He is a member of IEEE,
ISCA, IPSJ and ASJ.