Generative Statistical Models with Self-Emergent Grammar of Chord Sequences
Hiroaki Tsushima1 , Eita Nakamura1 , Katsutoshi Itoyama1 , and Kazuyoshi Yoshii1,2
1

Graduate School of Informatics, Kyoto University, Kyoto 606-8501, Japan
2
RIKEN AIP, Tokyo 103-0027, Japan

A weakness of commonly used Markov models is their
limitation in learning syntactic similarities among the output symbols. For example, in the first-order model, there is
no easy way to learn tied transition probabilities P ( · |x) =
P ( · |x0 ) for a pair (and in general a group) of symbols x and
x0 without knowing the pair/group a priori1 . Empirically,
on the other hand, we know that such similarities exist for
chord symbols: e.g. composers often substitute, in C major key, G with G7 or F with Dm and vice versa. This is
one motivation for introducing categories of chords called
harmonic functions (e.g. tonic, dominant, and subdominant) in harmony theories [2–6]. If we can appropriately
capture those similarities and effectively reduce model parameters, the overfitting problem would be alleviated and
models with better performances could be obtained.
In generative models for natural language processing
(NLP), syntactic categories like parts of speech (POS)
are often represented as latent variables, e.g. latent states
of hidden Markov models (HMMs) [35, 36] and nonterminals of probabilistic context-free grammar (PCFG) models [37–39]. As the grammatical structure is described in
terms of the latent variables, with a linear-chain (Markov)
process in HMMs and a tree generation process in PCFG
models, the syntactic similarity among words in each category is naturally captured. Furthermore, unsupervised
learning techniques can be applied to refining categories according to data [37,38] or totally inducing the latent grammar with unknown categories directly from data [35,36,39].
We use the terms self-emergent HMMs and self-emergent
PCFG models to indicate those models whose latent variables and grammar are completely inferred through unsupervised learning using data without annotations (e.g.
POS tags) as in the last case. The use of these models
opens the possibility of constructing generative chord sequence models incorporating syntactic categories of chords
that can be adapted to the size and the style of given data.
Analogies between grammatical structures of music and
language have been repeatedly discussed in the literature
[18, 40]. Especially, the phrase structure grammar [41]
has been adapted for chord sequences and harmonic functions have been represented as nonterminals in contextfree grammar (CFG) or similar tree-structured models
[18, 19, 23], and more recently, in their probabilistic extensions [24, 27]. In these models, the nonterminals and
production rules are manually determined according to

arXiv:1708.02255v3 [] 2 Mar 2018

Abstract
Generative statistical models of chord sequences play crucial roles in music processing. To capture syntactic similarities among certain chords (e.g. in C major key, between G and G7 and between F and Dm), we study hidden Markov models and probabilistic context-free grammar
models with latent variables describing syntactic categories
of chord symbols and their unsupervised learning techniques for inducing the latent grammar from data. Surprisingly, we find that these models often outperform conventional Markov models in predictive power, and the selfemergent categories often correspond to traditional harmonic functions. This implies the need for chord categories
in harmony models from the informatics perspective.
Keywords: statistical music language model; harmonic
(tonal) function; unsupervised grammar induction;
Bayesian inference; hidden Markov model; probabilistic
context-free grammar.
1

Introduction

It is widely accepted that harmony, or the progression of
chords, plays an influential role in tonal music [1–6]. In the
field of music informatics, computational models of chord
sequences are regarded as a necessary ingredient for applications including chord recognition [7–11], harmonisation [12–17], music analysis [18–24], and music arrangement [25–27]. In particular, generative statistical models
of chord sequences with higher predictive power are not
only useful for chord sequence generation [28] and harmony analysis [20, 24] but also important for improving
the accuracies of chord recognition [8, 9] and melody harmonisation [12, 14]. This situation parallels the case for
language processing [29–31].
An advantage of using statistical models is that the
model parameters can be learned from data. Note that
models with more precise structure do not always perform
better in practice due to the overfitting problem [31, 32].
Thus, the performance of statistical models like predictive
power should be discussed with the generalisation ability
and dependence on the size of training data. This is especially important for music modelling, for which largesize data are often not available. For example, the Beatles
dataset [33], which is currently one of the largest datasets
used for chord recognition, contains about 200 pieces. For
attempts to improve Markov models by applying smoothing techniques and/or incorporating musical knowledge,
see Refs. [7–9, 14, 34].

1 The number of possible ways to tie transition probabilities is
given by the partition number. For example, with 10 symbols, it
is 42 for a first-order model and about 1.9 × 108 for a second-order
model (with bigram contexts).

1

music theories or intuitions. Recent studies have shown
that chord categories similar to traditional harmonic functions can be obtained from data by unsupervised learning [42, 43], which implies that the generative grammar
could also be inferred from data.

ourselves to information-scientific analyses in this study,
we hope that the approach will be expanded to the area of
computational musicology as discussed in the Conclusion.
2

In this paper, we formulate self-emergent HMMs and
PCFG models of chord sequences whose latent variables
correspond to unknown chord categories and describe
linear-chain grammar in HMMs and tree-structured grammar in PCFG models, and study unsupervised learning
techniques for them. To evaluate the predictive power and
generalisation ability, we compare these models and conventionally used Markov models in terms of test-data perplexity and symbol-wise predictive accuracy for different
sizes of training data. We find that these models outperform Markov models in certain cases, typically with small
training-data sizes and small numbers of possible chord
symbols. Dependence on initialisation and learning scheme
is also investigated, which is known to be significant for
models’ performance in NLP [35, 36, 39]. We also carry
out information-theoretical analyses on the structure of obtained HMMs and confirm that they indeed capture syntactic similarities among chord symbols. Interestingly, the
learned categories often correspond to traditional harmonic
functions when the model size is small, which confirms the
results of unsupervised chord categorisation [42, 43].

Related Work

In this section, we review related studies on unsupervised
learning of models for music and natural language.
2.1

Unsupervised Learning of Music Models

Unsupervised learning of chord categories has been studied
in the context of music analysis. In one study [42], chords
(pitch class sets) taken from a corpus of Bach’s chorales
were clustered based on the hierarchical cluster analysis
according to the similarity in antecedent and consequent
transition probabilities of chords. In another study [43],
probabilistic categorisation of chord symbols was formulated based on combined optimisation of the perplexity of
categories and the predictive accuracy of chords. In both
studies, chord categories similar to traditional harmonic
functions were automatically obtained without referring to
the pitch contents, which indicates that harmonic functions
are closely related to chords’ similarities in the sequential
context. In the latter study, a supplementary analysis using a self-emergent HMM with three states has been carried
out and it was reported that the learned categories were
not consonant with traditional harmonic functions for some
data. These studies focused on determining the optimal
chord categorisation from data and did not argue whether
chord categories are necessary on the ground of predictive
modelling or music processing. For example, in the formulation of the latter study [43], any categorisation scheme
cannot exceed Markov models in the predictive power.
In Ref. [27], unsupervised learning was applied to a variant of PCFG models for chord symbols in which nonterminals describe harmonic functions. In this study, the nonterminals and the production rules were manually determined in accordance with existing data of harmony analysis and only the probabilities were learned unsupervisedly.
Our PCFG models can be considered as extensions of this
model as we allow larger numbers of nonterminals and they
have no prior labels or restricted production rules.
Unsupervised learning of self-emergent HMMs have been
applied for analyses of melodies [44, 45], which demonstrated that learned HMM states can capture the syntactic roles of musical symbols such as pitch and note value.
In the first study [44], pitch sequences of Greek church
chants were modelled with HMMs and it was found that
the learned HMM structures can be related to melodic cadential formulas. In the second study [45], sequences of
note values extracted from Palestrina’s masses were analysed similarly and it was found that the learned HMM
states corresponded to metric positions of note onsets. Although the EM algorithm was used for learning, whose results generally depend on initial parameters, the influence

Our main contribution is to quantitatively show the
potential of self-emergent HMMs and PCFG models as
generative statistical models of chord sequences, which
can be applied for multiple problems in music processing. As explained in Sec. 2, these models have been
explored in limited ways for music, and our study covers their Bayesian extensions and comparison between the
expectation-maximisation (EM) algorithm and Gibbs sampling (GS) algorithm for learning. Systematic investigations of the influence of initialisation and the relation between HMMs and PCFG models, which have practical importance, have not been discussed, to our knowledge, in
the context of music modelling. The theoretical part is
intended to include a concise review of commonly used
Markov models, HMMs, and PCFG models and their relations, which could be useful generally for music and language processing communities.
Our approach can be contrasted with music modelling
based on expert knowledge, which has been the mainstream of music analysis. The use of unsupervised learning removes the reliance on knowledge about functional
harmony theory etc. and the models can be automatically adapted to the size and style of given training data.
Whereas we use chord-sequence data in the popular/jazz
notation, which is of interest for applications in music processing, the method is independent of particular music
styles or systems of chord notation and is applicable to
other chord notations such as the Roman numeral notation and the figured bass notation. Although we confine
2

of initialisation was not discussed. The optimal number
of latent states was estimated by the minimum description
length principle, but the dependence on the training-data
size was not studied and there was no comparison with
Markov models or other models in terms of perplexity etc.
In the context of melody modelling, comparison of
Markov models, HMMs, and PCFG models have been carried out [46]. Zeroth- and first-order Markov models of
pitches and of pitch intervals, self-emergent HMMs with
up to 12 latent states, and PCFG models with pre-assigned
rules were formulated as Bayesian models and compared in
terms of variational free energy. The results showed that
the relative performance of the models depended on the
data: for Bach’s chorales a PCFG model was the best and
for a collection of folk songs the first-order Markov model of
pitch intervals was the best. Because higher-order Markov
models or larger-size HMMs were not tested and models
were learned and tested using the same data, the overfitting effect could not be observed. The influence of initialisation was also not discussed.
2.2

3

Models

Here we describe Markov models and self-emergent HMMs
and PCFG models for chord symbols and their learning methods. The problem of initialisation for unsupervised learning of self-emergent HMMs and PCFG models
is discussed and two learning methods, the expectationmaximisation (EM) algorithm and Gibbs sampling (GS)
method, are explained.
A sequence of chord symbols is written as
x = x1:N = (x1 , . . . , xN ) ,

(1)

where each xn takes values in the set of chord symbols
(output symbol space) Ω = {x}. The number of possible
symbols (symbol-space size) is denoted by NΩ = #Ω. In
this study, a generative statistical model of chord sequences
is defined as a model that yields the probability P (x1:N )
for any sequence x1:N .
3.1

Markov Model

The first-order Markov model is defined with an initial
probability (χini ) and transition probabilities (χ) as

Unsupervised Learning of Language Models

Self-emergent HMMs have been applied for unsupervised
POS tagging [35, 36]. In Ref. [36], Bayesian HMMs were
learned using GS and they were shown to yield significantly
higher tagging accuracies than HMMs trained by the EM
algorithm. In Ref. [35], it was shown that the performance
of HMMs learned by the EM algorithm is largely influenced
by the random initialisation and when the size of the state
space is small (e.g. 25), they yielded a similar tagging accuracy as the Bayesian HMMs. Since the tagging accuracy
was the focus of these studies, the perplexities of the output symbols (words) were not measured.
Unsupervised learning of PCFG models was studied in
the context of grammar induction and parsing [39]. The
maximum likelihood estimation using the EM algorithm
and Bayesian estimation using GS were compared and the
latter method was shown to give higher accuracies in some
cases. It was argued that even though for both learning
schemes fully unsupervised grammar induction using simple PCFG models is difficult, Bayesian learning seems to
induce more meaningful grammar owing to its ability to
prefer sparse grammars.
Unsupervised learning has also been applied to adapt
statistical language models for particular data. For NLP,
there exist corpora of syntactic trees that are based on
widely accepted linguistic theories (e.g. [41]), from which
production probabilities can be obtained by supervised
learning. On the other hand, studies have shown that extending the standard annotation symbols (e.g. NP) to incorporate context dependence and subcategories of words,
often called symbol refinement, improves the accuracy of
parsing [37]. Unsupervised learning of variants of PCFG
models has been successfully applied to find optimal refinement of nonterminal symbols from data [38].

P (x1 ) = χini
x1 ,

P (xn |xn−1 ) = χxn−1 xn .

(2)

The probability of x1:N is then given as
P (x1:N ) = χini
x1

N
Y

χxn xn−1 .

(3)

n=2

As a generalisation, a k th-order Markov model is defined
with the following initial and transition probabilities:
P (x1 ) = χini
x1 ,
P (x2 |x1 ) =
..
.

(4)

χini(2)
x1 x2 ,

(5)

ini(k)

(6)

P (xk | x1 , . . . , xk−1 ) = χx1 ···xk−1 xk ,

P (xn | xn−k , . . . , xn−1 ) = χxn−k ···xn−1 xn .

(7)

The probability P (x1:N ) is given similarly as in Eq. (3).
Higher-order Markov models are more precise than lowerorder Markov models. This is because for any k th-order
Markov model we can choose parameters for a (k+1)thorder Markov model so that these models yield the same
probability P (x1:N ) for any sequence x1:N .
Since each xn can take NΩ values, for the first-order
model, the initial probability has NΩ parameters and the
transition probability has NΩ2 parameters. These parameters are not independent, however, since the probabilities
must add up to unity. After normalisation, we have NΩ −1
free parameters for the initial probability and NΩ (NΩ − 1)
for the transition probability, and hence (1 + NΩ )(NΩ − 1)
free parameters in total. Similar calculations yield the
3

Model

We can also define higher-order HMMs by considering
higher-order Markov models for latent states. For selfemergent HMMs (i.e. if the state space is not fixed), however, first-order HMMs are general enough that higherorder HMMs can be described by using them. To see this,
consider a second-order HMM with a state space Γ = {z}
and transition probability P (zn |zn−1 , zn−2 ). If we define
a state space Γ̂ as the product space Γ × Γ and identify
ẑn ∈ Γ̂ as (zn , zn−1 ) ∈ Γ × Γ, the transition probability
P (ẑn |ẑn−1 ) is equivalent to P (zn |zn−1 , zn−2 ) since

Number of free parameters

Markov
(1 + NΩ + · · · + NΩk )(NΩ − 1)
HMM (1 + NΓ )(NΓ − 1) + NΓ (NΩ − 1)
2
PCFG
(1 + N∆ )(N∆
− 1) + N∆ NΩ
Table 1: Number of model parameters after normalisation
for a k th order Markov model, an HMM with NΓ states,
and a PCFG model with N∆ nonterminals and the constraint ψS→x ≡ 0. NΩ is the symbol-space size.

P (zn , zn−1 | zn−1 , zn−2 ) = P (zn | zn−1 , zn−2 ).

number of free parameters for a k th-order Markov model
as in Table 1.
Given a training dataset, the parameters of Markov
models can be learned based on the maximum likelihood
principle. For example, the transition probability P (x0 |x)
for the first-order Markov model is estimated as
P (x0 |x) ∝ C(xx0 ),

A similar argument is valid for the initial and output probabilities and for higher-order HMMs.
In the unsupervised case, i.e. when the training data
only has the information of output symbols but not corresponding latent states, the parameters of HMMs cannot be
learned as simply as the case for Markov models. Nevertheless, the EM algorithm can be applied for the maximum
likelihood estimation. The EM algorithm (also called the
Baum-Welch algorithm) is based on iterative optimisation.
Given training data X̄ = x1:N̄ (for notational simplicity, we
here consider that the training data has only one sequence;
the extension for data with multiple sequences is straightforward) and a set of parameters Π = (πzini , πzz0 , φzx ), the
updated parameters Π̃ = (π̃zini , π̃zz0 , φ̃zx ) are given by

(8)

where C(xx0 ) denotes the number of times that a bigram
xx0 appears in the training data. Other parameters can be
obtained similarly. Especially for small training data and
for high-order Markov models, however, smoothing techniques are necessary to avoid the data sparseness problem. For the additive smoothing, we simply add a positive constant number  to the right-hand side (RHS) of
Eq. (8) and similarly for other cases. Based on a previous study that compared various smoothing methods for
language data, the most effective smoothing method was
the (interpolated) Knesner-Ney (KN) smoothing and the
modified KN (MKN) smoothing, which incrementally use
lower-order transition probabilities to obtain higher-order
transition probabilities [31].

π̃zini = P (z1 = z | X̄, Π),
N̄P
−1

π̃zz0 =

n=1

n=1

3.2.1

HMM

φ̃zx =

Basic Model and the EM Algorithm

P (zn |zn−1 ) = πzn−1 zn ,

P (xn |zn ) = φzn xn .

N̄
P

,

(13)

P (zn = z | X̄, Π)

δxn x P (zn = z | X̄, Π)

n=1

In HMMs, we introduce a latent variable (state), denoted
by zn , corresponding to each symbol xn and consider a
first-order Markov model for the sequence of latent states.
Let Γ = {z} denote the set of possible latent states called
the state space and NΓ = #Γ denote its size. The initial
and transition probabilities will be denoted by π ini and
π. In addition, each xn is considered to be drawn from a
distribution depending on zn called the output probability
(φ). Thus, the model parameters are as follows:
P (z1 ) = πzini
,
1

n=1

(12)

P (zn = z, zn+1 = z 0 | X̄, Π)
N̄P
−1

N̄
P

3.2

(11)

.

(14)

P (zn = z | X̄, Π)

Here, the probabilities in the RHS are defined with the
parameters before the update and can be computed by
the forward-backward algorithm (see e.g. Ref. [30] for
details). These updates always increase the likelihood,
P (X̄|Π̃) ≥ P (X̄|Π), until it converges to a local maximum [32]. According to this monotonic increase of the
likelihood, the EM algorithm can also get stuck in a local
maximum and miss the global maximum.
There is no guarantee that the EM algorithm finds
the global maximum likelihood solution. Thus we have
the problem of initialisation for unsupervised learning of
HMMs. Possible approaches are (i) to start with several
random initial values and obtain the result with the maximum likelihood and (ii) to start with a potentially ‘good’
choice of parameters, e.g. parameters estimated by supervised learning with small-size data; parameters set by existing theories, intuition, etc.

(9)
(10)

The number of parameters is as shown in Table 1. We
interpret the latent states as chord categories. The transition probabilities describe the linear-chain grammar for
chord categories and the output probabilities relate the
categories to chord symbols.
4

3.2.2

output probabilities (ψ) relate the categories to chord symbols. The probabilities are normalised as
X
X
1=
θS→zL zR +
ψS→x ,
(18)

Bayesian Model and the GS Method

Bayesian learning can be applied for learning HMM parameters to alleviate the potential problems of the EM algorithm [35,36]. Firstly, with the introduction of the prior
distribution of the parameters, one can incorporate prior
knowledge or preferences into parameter learning. Secondly and more importantly for us, there are samplingbased algorithms that can search parameters without getting stuck in local optima.
To utilise Bayesian learning techniques, we extend the
model with prior distributions on the parameters. Since
the initial, transition, and output probabilities are all categorical distributions, we put the conjugate Dirichlet disini
tributions on their parameters π ini = (πzini )z , πz = (πzw
)w ,
and φz = (φzx )x as
π ini ∼ Dir(λini ),

πz ∼ Dir(λz ),

1=

φz ∼ Dir(ρz ), (15)

θ̃z→zL zR ∝

P (z → zL zR ) = θz→zL zR ,

P (z → x) = ψz→x .

ψz→x .

(19)

x∈Ω

X
T ∈T(X̄)

c(z → zL zR ; T )P (T | X̄, Θ)

(20)

x

zL ,zR

·

Y

c(z→zL zR ;T )
θz→z
L zR

z,zL ,zR

Y

c(z→x;T )
ψz→x

z,x

and the RHS of Eq. (20) can be computed by the insideoutside algorithm [29]. The discussion about initialisation
for HMMs is also valid for PCFG models.
3.3.2

Bayesian Model and the GS Method

The Bayesian extension for PCFG models is also similar as
the case for HMMs. We put conjugate priors for the model
parameters as

PCFG models are defined with a set of output symbols
(also called terminals) Ω = {x}, a set of latent variables
called nonterminals ∆ = {z} (its size is denoted by N∆ ), a
start symbol S, and the following set of production probabilities and output probabilities2 :
P (S → x) = ψS→x ,

θz→zL zR +

and similarly for θ̃S→zL zR , ψ̃S→x , and ψ̃z→x . Here, T denotes a derivation tree, T(X̄) is the set of possible derivation trees that yield X̄, and c(R; T ) is the number of times
that a production rule R is used in T . The probability
P (T | X̄, Θ) can be expanded as
Y c(S→z z ;T ) Y c(S→x;T )
P (T | X̄, Θ) =
θS→zL zLR R
ψS→x

Basic Model and the EM Algorithm

P (S → zL zR ) = θS→zL zR ,

X

Notice that the rule S → x is only used for sequences of
length unity. Whereas we retain this probability in this
section, which will be used for theoretical discussions, we
set ψS→x = 0 practically for numerical calculations because there were no chord sequences with length unity in
our datasets. With this constraint, the number of parameters is as shown in Table 1.
We can derive the EM algorithm for PCFG
models similarly as the case for HMMs.
Given
training data X̄ and a set of parameters Θ =
(θS→zL zR , θz→zL zR , ψS→x , ψz→x ), the updated parameters

Θ̃ = θ̃S→zL zR , θ̃z→zL zR , ψ̃S→x , ψ̃z→x are given as

PCFG Model

3.3.1

x∈Ω

X
zL ,zR ∈∆

where λini , λz , and ρz are Dirichlet parameters with dimensions equal to the corresponding distributions.
For Bayesian inference, we estimate the distribution of
the parameters P (Π| X̄, Λ) given the training data X̄ and
hyperparameters Λ = (λini , λz , ρz ). The Gibbs sampling
(GS) is a method to sample Z = z1:N̄ and Π from the
joint distribution P (Z, Π | X̄, Λ). The variables are alternately sampled from the distributions P (Z| Π, X̄, Λ) and
P (Π| Z, X̄, Λ) [36]. In theory, the Gibbs sampler yields
samples from the exact distribution P (Π| X̄, Λ). This
means that the sampled parameters explore globally in
the parameter space with more weight on the regions with
larger likelihood. In practice, we select the maximum likelihood parameters out of samples from the Gibbs sampler.
Because the sampled parameters are not always locally
optimal, we apply the EM algorithm using the selected
sample as initial values, to finally obtain locally optimal
parameters.
3.3

zL ,zR ∈∆

θS ∼ Dir(ξS ),

θz ∼ Dir(ξz ),

(16)

ψS ∼ Dir(ζS ),

ψz ∼ Dir(ζz ),

(21)
(22)

where θS = (θS→zL zR )zL zR , θz = (θz→zL zR )zL ,zR , ψS =
(ψS→x )x , and ψz = (ψz→x )x are the vector representation
of the probabilities and ξS , ξz , ζS , and ζz are Dirichlet
2
2
parameters, with dimension N∆
, N∆
, NΩ , and NΩ .
In the GS for PCFG models, we alternately sample
derivation trees and model parameters. Let us denote
the hyperparameters as Ξ = (ξS , ξz , ζS , ζz ). The variables are sampled from P (T | Θ, X̄, Ξ) and P (Θ| T, X̄, Ξ),

(17)

As the case for HMMs, we interpret the nonterminals as
chord categories. The production probabilities (θ) describe
the tree-structured grammar for chord categories and the
2 To make the relation between PCFG models and HMMs explicit,
we here do not include the start symbol in ∆.

5

S
z2

z3

···

···

construct an HMM that mimics any given Markov model.
For any first-order Markov model, we can indeed construct
such an HMM by setting Γ = Ω, π ini = χini , π = χ,
and φzx = δzx . Since any higher-order Markov model can
be represented as a first-order Markov model with an extended set of output symbols, as we saw in the discussion
for higher-order HMMs, this shows that the set of HMMs
includes Markov models of any order.
The relation between HMMs and PCFG models derives
from a fact that a linear-chain grammar can be realised as
a special case of tree-structured grammar as illustrated in
Fig. 1. Roughly speaking, for a given HMM with NΓ states,
we can define a PCFG model with N∆ = NΓ nonterminals
with the following probabilities

zN −1

z1

z2

z3

···

zN −1

zN

x1

x2

x3

···

xN −1

xN

Figure 1: Realisation of linear-chain grammar in the treestructured grammar.
respectively. The derivation trees can be sampled by the
inside-filtering outside-sampling method and the model parameters can be sampled from the posterior Dirichlet distributions in a standard way [39].
4

π
(1 − κ0 ),
θS→zL zR = πzini
L zL zR
X
ψS→x =
πzini φzx κ0 ,

Model Analysis Techniques

(24)

z

θz→zL zR = δzzL πzzR (1 − κ),

Methods for evaluating and analysing the learned models
are developed here. In Sec. 4.1, we discuss theoretical relations among the three types of models (Markov models,
HMMs, and PCFG models) and derive an initialisation
method for PCFG models using parameters of corresponding HMMs. In Sec. 4.2, we formulate several evaluation
metrics to measure the predictive power of the generative
models. In Sec. 4.3, we develop techniques for examining
the structure of self-emergent HMMs, which will be used in
Sec. 5.4.1 to confirm that those models can indeed capture
syntactic similarities among chord symbols.
4.1

(23)

ψz→x = φzx κ

(25)
(26)

with some constants 0 ≤ κ, κ0 ≤ 1, so that the derivation
tree in Fig. 1 yields the probability (for N > 1)
P (x1:N , z1:N |Θ) = κN (1 − κ)N −2 (1 − κ0 )
· πzini
φz1 x1
1

N
Y

(πzn−1 zn φzn xn ),

(27)

n=2

which reproduces the complete-data probability of the
HMM up to a constant factor involving κ and κ0 . To show
a strict inclusion, however, we need more careful derivation
because of the different normalisation conditions for HMMs
and PCFG models. In fact, the inclusion relation is strictly
true only when HMMs and PCFG models are appropriately reformulated as we show in Appendix A. For here,
we remark that the inclusion is not strictly true with our
definition of HMMs, since the probabilities in Eqs. (23)–
(26) can generate derivation trees with different structures.
It is also remarked that κ0 controls the probability of generating sequences of length unity. Since all chord sequences
used in our study have lengths larger than unity, we exclude this case and consider κ0 = 0.
The inclusion relation between HMMs and PCFG models can be applied to finding initial values for parameter
learning of PCFG models [47]. It has been reported that
the initialisation problem is more severe for PCFG models than HMMs because of the complex model structure
and the existence of many local optima [39]. Thus, if we
were given parameters of HMMs, the initialisation problem could be alleviated if we use them to set initial values
for PCFG models as in Eqs. (23)–(26). Several remarks
should be made here. P
First, we need to specify the constant κ. As we have
x ψz→x = κ from Eq. (26), this
parameter describes the probability that a nonterminal is
converted to an output symbol (rather than being split into

Relations among the Models

We here explain theoretical relations among Markov models, HMMs, and PCFG models, which provide insights
into the model analysis. We use the expression that a
model (say model A) can be mimicked by another model
(say model B) if for all possible parameter configurations
for model A there is always a parameter configuration for
model B that provides the same probability P (X; A) =
P (X; B) for all data X. We also say that a set of models
includes another set of models if any model in the latter
set can be mimicked by a model in the former set. In summary, we show that smaller HMMs/PCFG models with
less number of latent states or nonterminals can be mimicked by larger HMMs/PCFG models. In addition, HMMs
include Markov models and PCFG models include HMMs.
The first assertion is easy to see. The parameters for
an HMM with NΓ states can be extended to those for a
new HMM with NΓ +1 states by assigning πẑini = πzẑ = 0
for the newly added state ẑ. In order to complete the new
HMM, πẑz and φẑx must be specified, but these parameters
do not affect the probability P (x|Π) for any sequence x
since ẑ will be generated with zero probability. The same
argument is valid for PCFG models.
Next, to see that HMMs include Markov models, we can
6

two nonterminals), which affects the expected length hN i
of generated sequences. Indeed, we can derive the following
equation (valid for κ > 1/2 and κ0 = 0):
hN i =

2κ
,
2κ − 1

of the ground-truth symbol. According to the probability in Eq. (30), we can sort the candidate predictions and
calculate the rank r(xn ) of the ground-truth symbol. We
define the reciprocal of the mean reciprocal rank (RMRR;
the harmonic mean of the ranks) as

−1
1 X X
1
,
(32)
R=
|X|
r(xn )
x ∈x

(28)

which can be used to choose κ. Second, since the constraint
z = zL in Eq. (25) cannot be relaxed during the learning
process (for both EM and GS algorithms) once imposed,
we should modify the equation to increase the flexibility of
learning. For this, we can add a small positive constant η
to the RHS of Eq. (25) and renormalise the probability.

x∈X

which is commonly used for evaluating classification methods. Formulas for calculating Eq. (30) for the three types
of models we consider are derived in Appendix C.
4.3

4.2

Evaluation Measures

(29)

where |X| denotes the number of samples in X and P (X|Θ̄)
is called the evidence probability. As noted above, we need
to use the equivalent normalisation condition to compare
PCFG models with other models appropriately. Calculation of the normalised evidence probability for PCFG
models is explained in Appendix B.
The test-data perplexity is (the exponential of) the entropy rate of the test data measured with the learned model
parameters and defined with the evidence probability for
the whole symbol sequences. While it is a measure of
model’s predictive power in a general sense, this measure is
defined without making explicit predictions. We can also
consider quantities that measure how the learned model
can predict each symbol from the rest of the symbols in a
sequence. To define these, we first calculate the probability
P (xn = x |x¬n ) = P (xn = x | x1:n−1 , xn+1:N )

(30)

for each symbol xn in each sequence x = x1:N in X. Here,
x¬n denotes the sequence obtained from x by removing xn
and the dependence on Θ̄ is made implicit for simplicity.
We can choose the symbol that maximises this probability
as a prediction for xn and the symbol-wise predictive error
rate (or simply, error rate) E can be defined as
1 X X
E=
I(xn 6= argmax P (x|x¬n )),
|X|
x
x ∈x
x∈X

Information-Theoretical Analysis of HMMs

As we will see in Sec. 5, HMMs often perform better than
Markov models in terms of predictive power. Here we develop analysis techniques to examine how the structure of
self-emergent HMMs learned from data evolves with increasing sizes of the state space.
In principle, when the model size is increased, HMMs
can evolve in two ways to increase their predictive power
(or decrease the perplexity). One way is to give sparser
output probabilities for each state. Since the output probabilities determine the resolution of assigning symbols to
each state, their sparseness contributes to accurate prediction of symbols. The extreme case is achieved when
output probabilities have all zero probabilities except for
one symbol as in Markov models. The other way is to associate multiple states to one or a group of symbol(s) and
describe finer sequential dependence with distinct transition probabilities. In this way transition probabilities can
incorporate dependence on more precise contexts including
longer-range contexts, similarly as a higher-order HMM is
represented as a first-order HMM. Since these two directions are in a trade-off relation when the state space is
finite, how learned HMMs actually behave depends on (or
characterises) the nature of used data. Thus, although
the following methods are model analysis techniques, they
can also reveal the nature of chord sequence data from the
information-theoretical perspective.
Let us now derive quantities that characterise the structure of HMMs. We first introduce the stationary distribution of latent states P̂ (z), which is defined by the equilibrium equation as
X
P̂ (z) =
P̂ (w)πwz .
(33)

As a measure of evaluating statistical models, the test-data
perplexity (equivalent to the cross entropy) is commonly
used. Denoting the test data by X and the learned parameters by Θ̄, the test-data perplexity P is given as


1
ln P (X|Θ̄) ,
P = exp −
|X|

n

w∈Γ

(Recall the notation for HMMs in Sec. 3.2.1: z and w
denote latent states, Γ the latent state space, and π the
transition probability.) Because of the ergodic property of
HMMs, the stationary distribution can be seen as the unigram probability of latent states generated from the HMM,
and we use the notation P̂ to indicate that the probability
is considered ‘locally’ by ignoring the sequential dependencies. The first quantity we define is the perplexity of the

(31)

n

where I(C) is 1 if condition C is true and 0 otherwise. A
potential problem with the error rate is that it cannot take
into account how the prediction is close to the ground truth
if it is incorrect. An alternative measure can be defined by
ranking the multiple predictions and ‘averaging’ the ranks
7

stationary distribution PΓ given as
 X

PΓ = exp −
P̂ (z) ln P̂ (z) .

Rank

1 to 10 C,F,G,Am,Em,Dm7,G7,Dm,E7,Em7
11 to 20 B[,Am7,FM7,Fm,A[,D,C7,D7,A,A7
21 to 50 E[,E,CM7,Gsus4,Bm7−5,Gm,Csus4,Gm7
D[,Cm,F7,F],F]dim,A[dim,F]m7−5
Bm,Fm7,B,Cadd9 ,F]m,G7sus4,B[m
B7,Fadd9 ,Bm7,Cm7,Bdim,C6,B[7,G6

(34)

z∈Γ

This quantity is interpreted as the effective size of the state
space that is used to describe data.
The sparseness of the output probabilities can be measured by means of their entropy. We define the average
perplexity of output probabilities Pφ as
 X

X
Pφ = exp −
P̂ (z)
φzx ln φzx .
(35)
z∈Γ

P̂ (x)

=

P̂ (z)φzx
P̂ (x)

.

5.1

(37)

z∈Γ

This quantity describes how many states are associated to
each symbol effectively.
The last quantity we define is the average perplexity of
transition probabilities Pπ given as
 X

X
Pπ = exp −
P̂ (z)
πzw ln πzw .
(39)
z∈Γ

w∈Γ

This quantity measures the average sparseness of the transition probabilities, or the effective number of possible
states that can be reached from one state. Low perplexities
of transition probabilities are necessary for accurate prediction of latent states and thence that of symbols. This
quantity should be considered in relation with the stationary perplexity PΓ in Eq. (34), which ignores the sequential
dependence of state transitions.
5

98.20%

Datasets and Setup

Two datasets of chord sequences of popular music pieces
were used for comparative evaluation of the models. The Jpop data consisted of chord sequences of 3500 J-pop songs
that were obtained from a public web page 3 . All songs had
key information in the major mode and were transposed to
C major key. Each sequence in this dataset corresponds to
a whole piece and the average length was 119.5 (chords).
The computation time required for learning and inferring
PCFG models using such long sequences is practically too
large and only Markov models and HMMs were tested using this data. The other dataset was constructed by extracting chord sequences from the Billboard dataset [48],
which is sampled from the Billboard Hot 100 from 1958
to 1991 4 . Using the musical structure annotation, 1702
chord sequences with a length larger than 7 were obtained
by segmenting 468 popular pieces in units of sections such
as verse and chorus. Using the key information assigned to
each section, the tonic of each sequence was transposed to
C. Because the Billboard data has the information about
the tonic but not the mode, the obtained sequences included C minor sequences as well as more frequent C major sequences. The average length of chord sequences was
13.0, and we tested all types of models using this dataset.
For evaluation, we separate each of the datasets into
training data and test data. For the J-pop data 500 sequences and for the Billboard data one tenth of all sequences (i.e. 171 sequences) were randomly chosen as the
test data. The rest of the sequences (3000 sequences for
the J-pop data and 1531 sequences for the Billboard data)
were used as the training data. To examine the generalisation ability of the models, randomly chosen subsets of the
training data were prepared. For the J-pop data we prepared training datasets with 30, 300, and 3000 sequences
and for the Billboard data 30, 300, and 1531 sequences.

The average state association variety V is defined by first
calculating the entropy of P̂ (z|x) with respect to z for each
x and then averaging them with P̂ (x) (and finally exponentiating):
 X

X
V = exp −
P̂ (x)
P̂ (z|x) ln P̂ (z|x) .
(38)
x∈Ω

79.75%
90.92%

discussed. In Sec. 5.4, we examine the structure of selfemergent grammar of learned models. Implications of the
results are discussed in Sec. 5.5.

x∈Ω

z∈Γ

P̂ (z, x)

Proportion

Table 2: Most frequent chord symbols in the J-pop data.

(Recall that x denotes a chord symbol and Ω the set of possible chord symbols.) Intuitively, this quantity tells how
many symbols are assigned to each state on average. Another quantity can be obtained by considering the variety
of association between symbols and states in the opposite
way. The unigram probability of symbols generated by the
HMM P̂ (x) and the probability of states given a symbol
P̂ (z|x) can be expressed as follows:
X
P̂ (x) =
P̂ (z)φzx ,
(36)
P̂ (z|x) =

Chord symbols

Numerical Experiments

3 J-Total Music: http://music.j-total.net (Researchers who
wish to have access to the J-Pop data should contact the authors.)
4 The dataset was downloaded from the McGill Billboard Project
webpage:
http://ddmal.music.mcgill.ca/research/billboard
(Complete Annotation section).

Here we present numerical results. After describing the
data and the setup in Sec. 5.1, the results of model comparisons are presented in Secs. 5.2 and 5.3, where the influence of initialisation for unsupervised learning is also
8

Rank

Chord symbols

1 to 10 C,F,G,Am,B[,Cm,G7,Dm,Dm7,A[
11 to 20 Em7,E[,Am7,Em,F/C,C7,FM7,Fm
Cm7,C5
21 to 50 F7,Gsus4,CM7,D,A,Gm,D[,Fm7,NC
D7,G7sus4,C/E,G7sus4add9 ,C1,A[M7
Gm7,F],C/G,Fadd9 ,G/C,B[/C,A7
B[5,F9,Cadd9 ,C9,B,E7,F6,G/B

0.5416. For both HMMs and PCFG models, the Dirichlet
parameters used for the GS algorithm were all set to 0.1.

Proportion
58.42%

5.2

69.55%

5.2.1
84.89%

Comparison of Predictive Power

Let us first compare the performances of HMMs and
Markov models. The test-data perplexity P and the error rate E for the two datasets are given in Figs. 2, 3, 4,
and 5, where the number of sequences in the training data
denoted by NX was 30, 300, and 3000/1531 sequences for
the J-pop/Billboard data. The number of tested random
number seeds was 30 when the number of latent states
NΓ ≤ 20 and 10 otherwise. The displayed perplexity or
error rate indicates the best one among the tested random
number seeds.
For the J-pop data (Figs. 2 and 3), we see that for all NX
and for all NΩ (the number of possible chord symbols), the
lowest test-data perplexity by the HMMs was less than the
lowest value by the Markov models. For NX = 30 and 300,
one can find that the test-data perplexity for the HMMs
reaches a minimum for some NΓ between 1 and 100 and
it increases for larger NΓ , which is the overfitting effect.
For NX = 3000, both test-data perplexity and error rate
decreased up to NΓ = 100, implying that they could be
further reduced with larger size HMMs. For comparison
between the EM and GS algorithms, there were no significant differences between their test-data perplexities in
the regions without overfitting. Some differences between
the error rates can be found, but we do not observe clear
tendencies.
For the Billboard data (Figs. 4 and 5), we see a similar
tendency of overfitting for NX = 30 and 300 as in the case
for the J-pop data. The lowest test-data perplexity and
error rate obtained for HMMs were respectively less than
those for Markov models. However, for NX = 1531 and
for NΩ = 20+1 and 50+1, the results for the third-order
Markov models were less than or similar to the best results for the HMMs. In this case, both test-data perplexity and error rate decreased up to NΓ = 100 for HMMs
and similarly for Markov models up to the third order, so
both types of models could yield lower values by further increasing the model sizes. If the perplexities or error rates
for models with similar numbers of parameters are compared, HMMs consistently outperformed Markov models.
For this data, we find a tendency that the results of the
GS algorithm were better than those of the EM algorithm
for NX = 1531 and N∆ ≥ 50.

Table 3: Most frequent chord symbols in the Billboard
data.
For music processing, a subset of chord symbols is usually selected and used as inputs/outputs of the models. For
example, the set of 12 major triads and 12 minor triads
with enharmonically inequivalent roots is often used. Because some other chords appear more frequently than some
major/minor triads (e.g. G7 or Dm7 is usually used much
more frequently than F] or B[m), we choose several sets
of most frequent chord symbols for each dataset as listed
in Tables 2 and 3. We used the top 10, 20, and 50 chord
symbols and all other chord symbols were treated as ‘Other’
so that the size of the symbol space NΩ = #Ω was 10+1,
20+1, and 50+1. Note that the Billboard data is more
contaminated with chords in C minor or other keys (due
to temporary modulations) and the proportion of ‘Other’ is
relatively large.
We compared Markov models, HMMs, and PCFG models with different learning schemes described in Sec. 3. For
Markov models, we tested from first-order to third-order
models learned with the additive smoothing (with additive
constant 0.1), KN smoothing, and MKN smoothing. For
HMMs, we tested the following state-space sizes NΓ :
{1, 2, . . . , 9, 10, 15, 20, 25, 30, 40, . . . , 90, 100}.

Comparison between HMMs and Markov
Models

(40)

For PCFG models, the tested numbers of possible nonterminals N∆ were {1, 2, . . . , 9, 10, 15, 20}. We remark that
the upper limit on the order of Markov models and that
on the NΓ or N∆ for HMMs or PCFG models were chosen so that the computational cost is near the practical
limit. Each HMM and PCFG model was learned based
on the EM and GS algorithms with random initial parameters and we tested for each case with several random
number seeds to examine the effect of initialisation. For
HMMs (PCFG models), at most 500 (200) iterations were
carried out in the EM algorithm until the likelihood converges. Whether the relative difference of likelihoods was
less than 10−5 was used as the criterion of convergence.
In the GS algorithm, the maximum likelihood parameter
set was chosen after 500 (200) samplings and then applied
the EM algorithm with a maximum of 50 iterations. For
PCFG models, initial parameters obtained from HMM parameters as in Eqs. (23)–(26), which will be called HMM
initialisations, were also tested. We set η = 0.01/N∆ and
the value of κ determined from data using Eq. (28) was
9

Order of Markov models:
HMM size: 1
3

1st

12

2nd
30

10

3rd
100

HMM (EM)
HMM (GS)
Markov (Add)
Markov (KN)
Markov (MKN)

Test-data perplexity P

11
10

9
8

NX =30

7

NX =300

6
5
4

NX =3000
10

100

1000

10000

# of parameters
(a) NΩ = 10 + 1.
Order of Markov models: 1st
HMM size: 1
3
10
20

Test-data perplexity P

18
16

2nd
100

30

3rd

HMM (EM)
HMM (GS)
Markov (Add)
Markov (KN)
Markov (MKN)

14
12

NX =30

10

NX =300

8

NX =3000

6
10

100

1000

10000

100000

# of parameters
(b) NΩ = 20 + 1.
Order of Markov models:
HMM size: 1
3
10

Test-data perplexity P

30

25

1st
30

2nd

100

3rd

HMM (EM)
HMM (GS)
Markov (Add)
Markov (KN)
Markov (MKN)

20

15

NX =30
NX =300

10

5

NX =3000
100

1000

10000

100000

1x10

6

1x10 7

# of parameters
(c) NΩ = 50 + 1.

Figure 2: Test-data perplexities for the J-pop data. NX is the number of sequences used for learning and NΩ is the
number of possible chord symbols. For example, ‘NΩ = 10 + 1’ means that the most common ten chord symbols were
used, plus a symbol ‘Other’ for all other chords. The horizontal axis indicates the number of parameters given in Table 1
and the corresponding model size (the order of Markov models or the number of states) is indicated on the upper side.

10

Order of Markov models:
HMM size: 1
3

1st

2nd
30

10

0.8

3rd
100

0.75

Error rate E

0.7
0.65

NX =30

0.6

NX =300

0.55
HMM (EM)
HMM (GS)
Markov (Add)
Markov (KN)
Markov (MKN)

0.5
0.45
0.4

10

NX =3000
100

1000

10000

# of parameters
(a) NΩ = 10 + 1.
Order of Markov models: 1st
HMM size: 1
3
10

2nd
100

30

0.8

3rd

0.75

Error rate E

0.7

NX =30

0.65

NX =300

0.6

0.55
0.5

HMM (EM)
HMM (GS)
Markov (Add)
Markov (KN)
Markov (MKN)

0.45
10

NX =3000

100

1000

10000

100000

# of parameters
(b) NΩ = 20 + 1.
Order of Markov models:
HMM size: 1
3
10
0.8

1st
30

2nd

100

3rd

0.75

Error rate E

0.7

NX =30

0.65
0.6

0.55
0.5

NX =300

HMM (EM)
HMM (GS)
Markov (Add)
Markov (KN)
Markov (MKN)
100

NX =3000
1000

10000

100000

1x10 6

1x10 7

# of parameters
(c) NΩ = 50 + 1.

Figure 3: Error rates for the J-pop data. See the caption to Fig. 2 for an explanation of symbols.

11

Order of Markov models:
HMM size: 1
3

1st

2nd
30

10

14

3rd
100

Test-data perplexity P

HMM (EM)
HMM (GS)
Markov (Add)
Markov (KN)
Markov (MKN)

12
10
8

NX =30

6

NX =300

4

NX =1531
2

10

100

1000

10000

# of parameters
(a) NΩ = 10 + 1.
Order of Markov models: 1st
HMM size: 1
3
10

2nd
100

30

20

3rd
HMM (EM)
HMM (GS)
Markov (Add)
Markov (KN)
Markov (MKN)

Test-data perplexity P

18
16
14

NX =30

12
10

NX =300

8
6

NX =1531

4
10

100

1000

10000

100000

# of parameters
(b) NΩ = 20 + 1.
Order of Markov models:
HMM size: 1
3
10

1st
30

2nd

100

Test-data perplexity P

50
40

3rd
HMM (EM)
HMM (GS)
Markov (Add)
Markov (KN)
Markov (MKN)

NX =30

30
20

NX =300

10
0

NX =1531
100

1000

10000

100000

1x10

6

1x10 7

# of parameters
(c) NΩ = 50 + 1.

Figure 4: Test-data perplexities for the Billboard data. See the caption to Fig. 2 for an explanation of symbols.

12

Order of Markov models:
HMM size: 1
3

1st

0.65

2nd
30

10

100

3rd

0.6
0.55

NX =30

Error rate E

0.5

NX =300

0.45
0.4

0.35

0.25
0.2

NX =1531

HMM (EM)
HMM (GS)
Markov (Add)
Markov (KN)
Markov (MKN)

0.3

10

100

1000

10000

# of parameters
(a) NΩ = 10 + 1.
Order of Markov models: 1st
HMM size: 1
3
10

2nd
100

30

3rd

0.7

NX =30

Error rate E

0.6
0.5

NX =300

0.4
0.3

HMM (EM)
HMM (GS)
Markov (Add)
Markov (KN)
Markov (MKN)

0.2
10

NX =1531

100

1000

10000

100000

# of parameters
(b) NΩ = 20 + 1.
Order of Markov models:
HMM size: 1
3
10
0.9

1st
30

2nd

100

3rd

0.8

NX =30

Error rate E

0.7
0.6

NX =300

0.5
0.4
0.3
0.2

HMM (EM)
HMM (GS)
Markov (Add)
Markov (KN)
Markov (MKN)
100

NX =1531

1000

10000

100000

1x10

6

1x10 7

# of parameters
(c) NΩ = 50 + 1.

Figure 5: Error rates for the Billboard data. See the caption to Fig. 2 for an explanation of symbols.

13

14

0.8

13

Test-data perplexity P

0.85

12

Error rate E

0.75

11

0.7

0.6

0.5
5

10

15

20

25

Model with the
lowest train-data
perplexity

9

NX = 30
NX = 300
NX = 3000

0.55

0.45

NX = 30
NX = 3000

10

0.65

8
7
6

30

1

100

(a) NX = 30 and 3000.

(a) Test-data perplexity vs error rate.
14

0.85

13

NX = 300

Test-data perplexity P

0.8

12

0.75

Model with the lowest
train-data perplexity

Error rate E

11

0.7

10

0.65

9

0.6

8

NX = 30
NX = 300
NX = 3000

0.55
0.5
0.45
1.4

10

HMM size NΓ

Test-data perplexity P

1.6

1.8

2

2.2

7
6

2.4

1

10

100

HMM size NΓ

RMRR R

(b) NX = 300.

(b) RMRR vs error rate.

Figure 7: Test-data perplexities for the J-pop data with
different initialisations (NΩ = 20 + 1, GS). See the caption
to Fig. 2 for an explanation of symbols.

Figure 6: Correlations between test-data perplexities, error
rates, and RMRRs (reciprocal of the mean reciprocal rank;
see Eq. (32)) for the J-pop data (NΩ = 20 + 1, GS). See
the caption to Fig. 2 for an explanation of symbols.

fective for application to music processing. In comparison
between the EM and GS algorithms for learning HMMs,
the GS algorithm may perform better especially for a large
number of latent states.

The relation between the test-data perplexity and error rate and that between the error rate and RMRR obtained for differently learned HMMs are shown in Fig. 6,
for the case of NΩ = 20+1 for the J-pop data. We see that
these quantities are highly correlated especially for lower
values (corresponding to higher predictive power), indicating that, despite the different definitions, they consistently
work as measures of predictive power in the limit when
the training-data size and state-space size are large. In
Fig. 6(a), there are two branches of models for NX = 30,
one of which extends to high test-data perplexities with
small changes in the error rate. This branch corresponds
to the overfitting models with NΓ & 10 and is a result of a
general tendency that the overfitting has a relatively little
effect on the error rate, as seen in Figs. (2)–(5). Results
for other cases and for the Billboard data were similar.
In summary, we confirmed that especially for small
training-data sizes and small symbol-space sizes, the predictive power of HMMs exceeds that of Markov models.
For larger training-data size and larger symbol-space size,
HMMs might be able to still outperform Markov models if
we could increase the model size. However, increasing the
HMM size larger than NΓ = 100 is practically difficult due
to the large computational cost and thus may not be ef-

5.2.2

Influence of Initialisation for HMMs

The next problem is how to find the optimal parameter
value among those obtained from different initialisations
using only the training data. This is an important issue
since, in the practical situation of music processing, one
wants to find a model that works best for unseen data. A
natural strategy to estimate the optimal parameter values
is to choose the one that minimises the perplexity or error
rate evaluated on the training data.
Fig. 7 shows the distributions of test-data perplexities for
different initial values for the J-pop data with NΩ = 20+1
and NX = 30, 300, and 3000. These results were obtained
with the GS algorithm and the results with the EM algorithm were similar. In the figure, each point corresponds
to each case of initialisation and the model that had the
lowest perplexity evaluated on the training data for each
model size is indicated with a square. For NX = 30, the
dependencies on the initial values are relatively large and
the parameters that are optimal for the training data does
not necessarily coincide with those for the test data even
14

for small NΓ . For NX = 300, there are still some differences among various initial values but they are smaller than
those in the case of NX = 30. There is a clear tendency
that for each model size the model optimal in terms of the
training-data perplexity is also close to optimal in terms of
the test-data perplexity up to about NΓ = 30, above which
the overfitting occurs. For NX = 3000, the differences
among various initial values are small for larger model sizes
and the correspondence between the optimal values for the
training and test data is evident. Even though it was more
obscure, there was a similar tendency for the error rates.
From these results, we conclude that although the best
way to choose the optimal parameter set is to directly
evaluate HMMs on the ground-truth test data, the optimisation using the training data works reasonably: when
the training-data size is small the optimisation can predict
wrong choices but with increasing data size the accuracy
of optimisation increases and the effect of incorrect choices
becomes smaller. For the data we studied, training data of
size NX ≥ 300 were able to accurately predict the optimal
model for unseen data for model sizes without overfitting.
5.3

Test-data perplexity P

14
12

PCFG (GS)
PCFG (HMM init, GS)
HMM (GS)

NX =30

10
8

NX =300

6
4

NX =1531
1

10

Model size NΓ /N∆

100

(a) Test-data perplexity.
0.7

NX =30

Error rate E

0.6

NX =300

0.5

NX =1531

0.4

PCFG (GS)
PCFG (HMM init, GS)
HMM (GS)

0.3
1

Comparison between PCFG Models and
HMMs

10

Model size NΓ /N∆

100

(b) Error rate.

Let us now compare the results for HMMs and PCFG models tested on the Billboard data. For PCFG models, the
EM and GS algorithms were run with 30 random number
seeds for model sizes N∆ ≤ 6 and between 4 and 6 random
number seeds for N∆ ≥ 7 depending on the computation
time required for learning. To learn PCFG models using the HMM initialisation, the optimal HMM parameter
values by the GS algorithm were chosen for each case according to the test-data perplexity and then the EM and
GS algorithms were applied.
In the following, we show and discuss only the results
for NΩ = 20+1 but the results for other cases were similar.
The test-data perplexities and error rates for the PCFG
models trained with the GS algorithm with random and
HMM initialisations are shown in Fig. 8, where for random
initialisations the best values are shown. The results for
HMMs with the GS algorithm are also shown as references.
For test-data perplexities, results for PCFG models
showed a clear overfitting effect in the case of NX = 30.
For NX = 300 and 1531, the PCFG models worked at
most equivalently and often worse than the HMMs with
NΓ = N∆ . There is a clear tendency that the HMM initialisations yielded better results than the random initialisations in the range 4 ≤ N∆ ≤ 10. These results show the
difficulty of parameter optimisation for PCFG models and
the effectiveness of the HMM initialisation in some cases.
Notice that for NX = 30 and 300, the lowest test-data
perplexity for the PCFG models was less than that of the
Markov models.
For error rates, we again find that overall the PCFG
models worked at most equivalently and often worse than

Figure 8: Test-data perplexities and error rates for PCFG
models trained with the GS algorithm with random and
HMM initialisations (Billboard data, NΩ = 20 + 1). The
results for HMMs trained with the GS algorithm are shown
as references. See the caption to Fig. 2 for an explanation
of symbols.
the HMMs with NΓ = N∆ . Even for the case of NX =
300 and 1531, where the overfitting effect was not eminent
for the test-data perplexity, we see considerable variations
among error rates for different N∆ s. This is partly due
to the limited number of random number seeds and also
reflects weak correlation between test-data perplexities and
error rates. The lowest error rate for the PCFG models was
less than that of the Markov models for NX = 30 and they
were even for NX = 300.
Let us look more closely at the influence of initialisation
and compare the EM and GS algorithms (Fig. 9). For random initialisations, the overall tendency is similar as the
case of HMMs: variations in test-data perplexities tend to
decrease for larger NX and, for NX = 300 and 1531, the
optimal model according to training-data perplexity also
has low test-data perplexity (except for the cases of overfitting). For NX = 300 and 1531 and N∆ ≥ 5, there is a
clear tendency that the GS algorithm often found better
parameters than the EM algorithm, even though the number of samples were limited and there were exceptions. On
the other hand, for HMM initialisations, differences between the results for the GS and EM algorithms were often negligible. These results again confirm that there are
15

16

EM
GS
Lowest train-data perp.
HMM init. (EM)
HMM init. (GS)

18
16

Test-data perplexity P

12

12

Perplexity

14

NX =30

10
8

10
8
6

NX =1531

6
4

Pφ + Pπ
Pφ
Model with the lowest
test-data perplexity

14

4

1

10

2

PCFG model size N∆

1

(a) NX = 30 and 1531.
3.5

EM
GS
Lowest train-data perp.
HMM init. (EM)
HMM init. (GS)

Test-data perplexity P

3

14

Perplexity

12
10

NX =300

8
6
4

2.5
2

Pπ
Model with the lowest
test-data perplexity

1.5

1

PCFG model size N∆

10

1

1

100

(b) Pπ .

Figure 9: Test-data perplexities for PCFG models with
different initialisations (Billboard data, NΩ = 20 + 1). See
the caption to Fig. 2 for an explanation of symbols.

100

Perplexity/variety

many local optima for the PCFG-model parameters and
it is difficult to find good optima by random initialisation
even with the GS algorithm.

5.4.1

10

HMM size NΓ

(b) NX = 300.

5.4

100

(a) Pφ and Pφ + Pπ .

18
16

10

HMM size NΓ

Structure of Self-Emergent Grammar

PΓ
V
Model with the lowest
test-data perplexity
10

Information-Theoretical Analysis

Let us now examine the learned HMMs using the techniques developed in Sec. 4.3. The quantities Pφ (average
perplexity of output probabilities), Pπ (average perplexity of transition probabilities), Pφ + Pπ , PΓ (perplexity
of the stationary distribution), and V (average state association variety) calculated for HMMs learned with the
J-pop data with NΩ = 20 + 1 and NX = 3000 are shown
in Fig. 10. First, the perplexity of the stationary distribution PΓ increases roughly proportionally to the model
size. This means that all latent states are exploited almost
uniformly to capture the sequential dependence.
Second, the average perplexity of output probabilities
Pφ tends to decrease with increasing model sizes, which
confirms our expectation that the larger-size models have
sparser output probabilities to increase the predictive
power. At the same time, Pφ remains larger than 2 even for
NΓ = 100. As the HMMs could have Pφ = 1 for NΓ ≥ 21,

1

1

10

HMM size NΓ

100

(c) PΓ and V.

Figure 10: Information-theoretical analyses of HMMs
trained with the J-pop data (NΩ =20+1, GS, NX =3000).
See Sec. 4.3 for the definition of the information measures
Pφ , Pπ , PΓ , and V.
this is an evidence that the HMMs automatically captured
syntactic similarities among different chord symbols, which
was preferred in order to reduce the perplexity. On the
other hand, the increase of the average state assign variety
V indicates that the same chord symbols are represented
by multiple latent states depending on their context, which
has a similar effect as symbol refinement for natural lan16

State ID

(Interpretation)

0.38
Prob

1

Tonic

Am Em Am7 Em7

F

Dm7 Dm FM7 D7 Gsus4 D Bm7−5 Fm

G

E7

A7

1

C7 CM7 Other B

C

Tonic

Prob

2

Subdominant

A

Dominant

G7

E7

C

B

Other

Fm Em Em7 Other E

B

Others

G

A

A

F

D

E

Cm D

State ID

4

Others 1

5

Others 2

F

Dm Other Dm7 FM7 F/C Am

G

Other

Em CM7 C/E

A

A7

Prob

B

D

D7 Fadd9 G/B

Prob

F7

G7 Gsus4 Dm7

add9 G7sus4 NC Fm E7
G7sus4

G/C FM9

Prob

3

Dominant

G

Am7 Em7 C7

Am

NC Fadd9 FM7 E7 Cadd9 C5

C

CM7

A

F

Fm

Dm Cm7 C/G

Prob

2

Subdominant

Prob

(Interpretation)

C

Other

Cm

B

A

E

Cm7 C5

D

A M7 C1

(b) PCFG model (NΩ = 50 + 1, N∆ = 5, GS, NX = 1531).

Figure 11: Output probabilities of (a) an HMM learned
with the J-pop data and (b) a PCFG model learned with
the Billboard data. States are manually sorted for clearer
illustration and for each state 12 chord symbols with highest output probabilities are shown.
guage models [37, 38].
Lastly, the perplexity of transition probabilities Pπ only
slowly increases for larger model sizes compared to the increase of PΓ . This is another evidence that the larger-size
models tend to capture longer-range sequential dependence
by using different states to describe different contexts. We
also see that the sum Pφ + Pπ , which serves as an approximation of the training-data perplexity, decreases for
larger model sizes. Among models with different initialisations, there is a tendency that the one that minimises
the test-data perplexity has low Pφ + Pπ , which indirectly
implies that the optimal model is determined by balancing
low perplexities of output and transition probabilities.
5.4.2

3

Dominant

4

Others

ities do not strictly categorise or discriminate chord symbols, we can identify each latent state or nonterminal as a
‘category’ that associates with some chord symbols more
strongly than others. Fig. 11 shows examples of output
probabilities for an HMM trained with the J-pop data and
a PCFG model trained with the Billboard data. The first
three states/nonterminals in these examples can be identified as the tonic, subdominant, and dominant functions.
We can see the general tendency that extended chords belong to the same category as a triad with the same root
note (e.g. G and G7, and Dm and Dm7), and that chords
in the parallel relation (e.g. C and Am, and F and Dm)
belong to the same category. An interesting observation
is that Em is more strongly associated to the ‘tonic’ category than the ‘dominant’ category in both cases. Note
that these categories are obtained without referring to the
constituent pitches of chord symbols. In both examples,
there are one or two categories that contain the symbols
Other and NC as well as relatively rare chords in C major
key (e.g. B[, A[, Cm).
Fig. 12 shows the transition probabilities of the 4-state
HMM trained with the J-pop data (same as in Fig. 11(a)),
where only those probabilities larger than 0.05 are shown
for clarity. This Markov model can be interpreted as a
probabilistic representation of the familiar cadence structure, with strong dominant motion (‘Dominant’ to ‘Tonic’)
and rare motion from ‘Dominant’ to ‘Subdominant’.
Depending on initial parameters, for both HMMs and
PCFG models, the learned categories were not always related to traditional harmonic functions as clearly as the
above examples. This can be understood by the existence
of multiple local optima in the parameter space as discussed in Sec. 5.2.2. This may also explain the report [43]
that HMMs did not yield chord categories consonant with
traditional harmonic functions for some data, since the influence of initialisation was not studied there. It would be
interesting to study further how the traditional harmonic
functions can be characterised as a chord categorisation
shceme based on information measures and whether alternative forms of chord categorisation are possible.

(a) HMM (NΩ = 50 + 1, NΓ = 4, GS, NX = 3000).

1

0.14

0.96

Figure 12: Learned transition probabilities of the HMM in
Fig. 11(a). Only probabilities > 0.05 are indicated.

Gm7 Csus4 Gm

C7

0.33

2

Subdominant

0.54

0.82

Prob

4

G

0.12

0.45

B G7sus4

Prob

3

Tonic

0.16

Examples of Learned Chord Categories

An interesting feature of the unsupervised learning of
HMMs and PCFG models is that the self-emergent grammars often represent harmonic functions of chords in traditional harmony theories. Although the output probabil-

5.5

Discussion

As is confirmed for speech recognition, the perplexity of
the language model in a recognition system is directly re17

Chord
categories

zn−1

π

zn

zn+1

xn

xn+1

φ
Chord
symbols

xn−1

χ

xn

Chord
symbols

xn+1

xn−1

ω
Features

yn−1

ω
yn

yn+1

Features

(a) Chord Markov model + feature model

yn−1

yn

yn+1

(b) Chord HMM + feature model

Figure 13: Integration of chord models with feature models for music processing.
lated to the recognition accuracy [31]. In our case, the
proposed models are expected to improve accuracies for
chord recognition and melody harmonisation. In these applications, the chord model is usually combined with an
acoustic model or a melody model, which are described as
output probabilities ωxy = P (yn = y | xn = x) of some feature y (e.g. acoustic features or melodic notes) conditioned
on a chord symbol. If one uses Markov models as the chord
model, the total model becomes an HMM. If one uses our
HMMs as the chord model, the total model is still an HMM
but with an additional layer of latent variables (Fig. 13).
To obtain the optimal chord sequence given the observed
features y = y1:N , the following probability is used:
P
P (z, x, y)
P (x, y)
= z
.
(41)
P (x|y) =
P (y)
P (y)

improves the predictive power compared to the case without introducing categories for some data. Thus our result
implies the need for introducing chord categories in harmony models from the informatics viewpoint.
6

Conclusion

We have studied self-emergent HMMs and PCFG models
as generative statistical models for chord symbols whose latent grammatical structure can be learned unsupervisedly
from data. We have investigated the influence of trainingdata size and the model size and found that these models
outperform Markov models in terms of predictive power
particularly for small training-data sizes and small numbers of possible chord symbols, which are typical situations
in applications for music processing. By analysing selfemergent grammars of the models, we have confirmed that
they tend to capture syntactic similarities of chords and
the learned categories often corresponded to traditional
harmonic functions when the model size is small.
These models can be useful as a prior language model of
chord sequences for music processing such as chord recognition and melody harmonisation, as they are expected
to improve the accuracy compared to conventionally used
Markov models. In addition, as the notion of harmonic
functions have been used by musicians to learn harmonic
grammar and to create/analyse music, these models can
also be useful for automatic composition/arrangement and
computational music analysis. For example, the learned
model could indicate users for possibilities for chord substitutions in the process of creating chord progressions.
It is certainly interesting to apply the models for other
data in different styles and to compare their syntactic
structures. Because it is based on unsupervised learning,
the proposed modelling framework can be even more useful for uncommon music styles for which data exists but no
expert knowledge is developed. It would lead to the development of style-specific harmony theories and quantitative
and objective comparisons among them. This would open
a new way of data-driven computational musicology based
on the information-scientific aspect of data.

This probability can be computed by GS or other sampling methods. The situation is similar when one uses our
PCFG models as the chord model. Application for melody
harmonisation using the PCFG models extended with a
rhythmic model has been studied in Ref. [49].
With regard to computational cost, the PCFG models
require time complexity of square order with respect to
the sequence length and cubic order with respect to the
number of nonterminals, which is typically much larger
than that of HMMs. Given the result that PCFG models
did not much improve the predictive power of chord sequences compared to HMMs, there would be little interest
in using PCFG models for simple recognition tasks. Nevertheless, they may have importance in music analysis and
composition/arrangement tasks [19, 23, 27, 49], where the
tree structured grammar is considered useful.
Finally, we discuss implications of our result compared
to those in previous studies [42, 43] that also investigated
unsupervised learning of chord categories. These studies
showed that chord categories similar to traditional harmonic functions can be learned from data and our result is
a confirmation of this result. A difference is that whereas
in previous studies the use of chord categories was assumed
and the model learning was based on comparison between
different possibilities of categorisation, in our study we
showed that the introduction of chord categories in fact
18

S

Acknowledgment

z2

E. Nakamura is supported by the JSPS Postdoctoral Research Fellowship. This work is in part supported by JSPS
KAKENHI Nos. 24220006, 26280089, 26700020, 15K16054,
16H01744, 16H02917, 16K00501, and 16J05486 and JST
ACCEL No. JPMJAC1602.
A

Embedding of HMMs in PCFG Models

0
πzw
+ πz0 end

w∈Γ

(∀ z ∈ Γ).

1=

X

PHMM (x1:N end, z1:N ).

(43)

B

0
θz→ũw = δzu πzw
,

ψz̃→x = φzx ,

ψz→x =

X

πzini πz0 end φzx ,

z∈Γ
πz0 end φzx ,

z̃3

···

z̃N −1

zN

x1

x2

x3

···

xN −1

xN

(48)

Normalised Evidence Probability for PCFG
Models

As mentioned in Sec. 4.2, we need to appropriately normalise the evidence probability P (X) for PCFG models
to define the test-data perplexity. Here, we consider the
case of a single sequence x, which is to be applied for each
sequence in the test data. The evidence probability for
PCFG models is normalised as
X
P (y) = 1,
(49)

(44)

This can be easily derived from Eqs. (9), (10), and (43).
Now let us construct a PCFG model with 2NΓ nonterminals (plus the start symbol S), where NΓ = #Γ is the
state-space size of the HMM. For each z ∈ Γ, we label half
of the nonterminals with z and the other half with z̃. We
set the probabilities for this PCFG model as follows (here,
z, w, u ∈ Γ and x ∈ Ω):
ψS→x =

z̃2

which is what we wanted to obtain.
The resemblance between Eqs.(23)–(26) and Eqs.(45)
and (46) is manifest, where the probability πz0 end is simplified as κ and κ0 in the former case. Notice that we
needed twice as many nonterminals to constrain the generated derivation trees to those equivalent to the linear-chain
grammar. This significantly increases the computational
cost for learning and inference as we discussed in Sec. 5.5.
On the other hand, as we argued in Sec. 4.1, when we initialise PCFG model parameters by HMM parameters, it
is necessary to relax strict constraints to allow model parameters to explore a wider space for learning. This is the
practical reason why we used the approximate relation for
initialising PCFG models.

N =1 x1 ,...,xN ∈Ω
z1 ,...,zN ∈Γ

0
θS→z̃w = πzini πzw
,

z̃1

PHMM (x1:N end) = PPCFG (x1:N ),

The difference to the HMMs used in the main part is the
existence of πz0 end , which controls the expected length of
generated sequences, and the following normalisation condition holds for this new HMM:
∞
X

zN −1

probability PPCFG (x) equal to that of the HMM. That is,
for any positive integer N and any sequence x1:N , we have

That is, we formally add the symbol ‘end’ at the end of each
sequence. As a model that generates this type of data, the
state space of an HMM is now extended to Γ0 = Γ ∪ {end},
where the added state ‘end’ only outputs ‘end’. The initial
and output probabilities are defined without any changes,
but the transition probabilities are extended as a function
π 0 : Γ × Γ0 → R≥0 satisfying
X

···

···

Figure 14: Derivation tree of a PCFG model that mimics
an HMM.

In this appendix, we explain a strict embedding of HMMs
in PCFG models, which is mentioned in Sec. 4.1. The
result here validates the approximate embedding as in
Eqs. (23)–(26) and explains why we use such approximations in our study.
To show the embedding relation, we reformulate HMMs.
Instead of Eq. (1), we represent a data sequence in the
following form:
x1 · · · xN end .
(42)

1=

z3

y∈Ω∗

where the sum is taken over all sequences of any length
(Ω∗ = {y1:N |N ≥ 1}). What we want is an evidence probability P̄ (x) that is normalised in the set of sequences of a
fixed length, as in the case for Markov models and HMMs.
Let N (x) denote the length of a sequence x. We want

(45)
(46)

X

(47)

P̄ (y) = 1,

(50)

y∈ΩN (x)

and other components that are not displayed are all zero. It
can be easily checked that this PCFG model only generates
derivation trees of the form in Fig. 14 and has the evidence

where ΩN denotes the set of all sequences of length N .
Now, let P (N ) be the probability that any sequence of
19

length N would be generated, which is given as
X
P (y).
P (N ) =

For HMMs, let us notate the forward and backward variables as
(51)
αn (zn ) = P (zn , x1:n ),

y∈ΩN

P (zn−1 , zn , x) = P (zn−1 , zn , x¬n )P (xn |zn ),

P (zn−1 , zn , x¬n )

= P (zn−1 , x1:n−1 )P (zn |zn−1 )P (xn+1:N |zn ),

x∈Ω

X

m−1
X

(56)

which are obtained by the forward and backward algorithms [30]. With the Markovian assumptions, we have

Then P̄ (x) is given as P (x)/P (N (x)).
It remains to calculate P (N ) for a given N . We define
modified inside variables B̄nm (z), which are inductively
defined through the following modified inside algorithm:
X
B̄nn (z) =
ψz→x ,
(52)
B̄nm (z) =

βn (zn ) = P (xn+1:N |zn ),

(57)
(58)

and hence the following formula:
θz→zL zR B̄nk (zL )B̄(k+1)m (zR ).

zL ,zR ∈∆ k=1

P (xn = y |x¬n ) ∝

(53)

=

We then have P (N ) = B̄1N (S).

X

P (zn−1 , zn , x)
xn =y

zn−1 ,zn ∈∆

X

πxn−1 y φzn y αn−1 (zn−1 )βn (zn )

zn−1 ,zn ∈∆

C

Formulas for Symbol-Wise Prediction
Probabilities

for n > 1. The case for n = 1 is similar.
For PCFG models, the inside and outside variables are
defined as

Here, we derive formulas for the symbol-wise prediction
probability P (xn = y |x¬n ) in Eq. (30) for Markov models, HMMs, and PCFG models. Since we consider this
probability as a function of y, we have
P (xn = y |x¬n ) ∝ P (xn = y, x¬n ),

Bnm (z) = P (z ⇒ xn:m ),

Anm (z) = P (S ⇒ x1:n−1 z xm+1:N ),

P (xn |x¬n ) ∝

χxn0 −k:n0 −1 xn0

Ann (z) = P (S ⇒ x1:n−1 z xn+1:N )

(55)

and hence the following formula:
X
P (xn = y |x¬n ) ∝
ψz→y Ann (z).

n0 =n

for k + 1 ≤ n ≤ N − k, which will be a general formula.
For samples near the beginning or the end of the sequence,
we need a modification to this formula. For example, for
a first-order Markov model,
P (xn = y |x¬n ) ∝ χxn−1 y χyxn+1

P (x1 = y |x¬1 ) ∝

χini
y χyx2 ;

(61)

(62)

z∈∆

References

(2 ≤ n ≤ N − 1);

[1] J.-P. Rameau, Treatise on Harmony, Dover Publications, 1971. (The original French version was first
published in 1722).
[2] H. Riemann, Harmony Simplified: Or the Theory of
the Tonal Functions of Chords (2nd ed.), Augener,
1896.
[3] A. Schönberg, Structural Functions of Harmony (revised ed.), W. W. Norton & Company, 1969.
[4] W. Maler, Beitrag zur Durmolltonalen Harmonielehre
I (7th ed.), F. E. C. Leuckart, 2007.
[5] S. Kostka, D. Payne, and B. Almén, Tonal Harmony
with an Introduction to Twentieth-Century Music (7th
ed.), McGraw-Hill, 2013.
[6] J. Mulholland and T. Hojnacki, The Berklee Book of
Jazz Harmony, Berklee Press, 2013.

P (xN = y |x¬N ) ∝ χxN −1 y ,
and for a second-order Markov model,
P (xn = y |x¬n ) ∝ χxn−2 xn−1 y χxn−1 yxn+1 χyxn+1 xn+2
P (x1 = y |x¬1 ) ∝

(60)

where the symbol ‘⇒’ means that the right symbol/sequence is derived from the left symbol/sequence.
These variables can be computed by the inside and outside algorithms [29]. Especially, we have

(54)

which will be used to derive the following formulas.
For a k th-order Markov model, we have
n+k
Y

(59)

(3 ≤ n ≤ N − 2);

ini(2)
χini
y χyx2 χyx2 x3 ;

P (x2 = y |x¬2 ) ∝ χini(2)
x1 y χx1 yx3 χyx3 x4 ;

P (xN −1 = y |x¬(N −1) ) ∝ χxN −3 xN −2 y χxN −2 yxN ;

P (xN = y |x¬N ) ∝ χxN −2 xN −1 y .

The formulas for higher-order models are similar.
20

[7] H. Papadopoulos and G. Peeters, “Large-Scale Study
of Chord Estimation Algorithms Based on Chroma
Representation and HMM,” in Proc. CBMI, pp. 53–
60, 2007.

[23] M. Rohrmeier, “Towards a Generative Syntax of
Tonal Harmony,” J. Math. Music, 5(1), pp. 35–53,
2011.
[24] W. Granroth and M. Steedman, “Statistical Parsing
for Harmonic Analysis of Jazz Chord Sequences,” in
Proc. ICMC, pp. 478–485, 2012.
[25] S. Fukayama, K. Yoshii, and M. Goto, “ChordSequence-Factory: A Chord Arrangement System
Modifying Factorized Chord Sequence Probabilities,”
in Proc. ISMIR, pp. 457–462, 2013.
[26] D. Quick, “Composing with Kulitta,” in Proc. ICMC,
pp. 306–309, 2015.
[27] D. Quick, “Learning Production Probabilities for Musical Grammars,” J. New Music Res., 45(4), pp. 295–
313, 2016.
[28] D. Conklin, “Chord Sequence Generation with Semiotic Patterns,” J. Math. Music, 10(2), pp. 92–106,
2016.
[29] C. Manning and H. Schütze, Foundations of Statistical
Natural Language Processing, MIT Press, Cambridge,
1999.
[30] L. Rabiner, “A Tutorial on Hidden Markov Models
and Selected Applications in Speech Recognition,”
Proc. IEEE, 77(2), pp. 257–286, 1989.
[31] S. F. Chen and J. Goodman, “An Empirical Study of
Smoothing Techniques for Language Modeling,” Computer Speech and Language, 13, pp. 359–394, 1999.
[32] C. M. Bishop, Pattern Recognition and Machine
Learning, Springer, 2006.
[33] C. Harte, M. B. Sandler, S. A. Abdallah, and
E. Gómez, “Symbolic Representation of Musical
Chords: A Proposed Syntax for Text Annotations,”
in Proc. ISMIR, pp. 66–71, 2005.
[34] K. Yoshii and M. Goto, “A Vocabulary-Free InfinityGram Model for Nonparametric Bayesian Chord Progression Analysis,” in Proc. ISMIR, pp. 645–650,
2011.
[35] M. Johnson, “Why Doesn’t EM Find Good HMM
POS-Taggers?,” in Proc. EMNLP-CoNLL, pp. 296–
305, 2007.
[36] S. Goldwater and T. Griffiths, “A Fully Bayesian Approach to Unsupervised Part-of-Speech Tagging,” in
Proc. ACL, pp. 744–751, 2007.
[37] M. Johnson, “PCFG Models of Linguistic Tree Representations,” Computational Linguistics, 24, pp. 613–
632, 1998.
[38] T. Matsuzaki, Y. Miyao, and J. Tsujii, “Probabilistic
CFG with Latent Annotations,” in Proc. ACL, pp.
75–82, 2005.
[39] M. Johnson, T. Griffiths, and S. Goldwater, “Bayesian
Inference for PCFGs via Markov Chain Monte Carlo,”
in Proc. HLT-NAACL, pp. 139–146, 2007.

[8] M. Mauch and S. Dixon, “A Discrete Mixture Model
for Chord Labelling,” in Proc. ISMIR, pp. 45–50,
2008.
[9] R. Scholz, E. Vincent, and F. Bimbot, “Robust Modeling of Musical Chord Sequences Using Probabilistic
N-Grams,” in Proc. ICASSP, pp. 53–56, 2009.
[10] N. Boulanger-Lewandowski, Y. Bengio, and P. Vincent, “Audio Chord Recognition with Recurrent Neural Networks,” in Proc. ISMIR, pp. 335–340, 2013.
[11] M. McVicar, R. Santos-Rodrı́guez, Y. Ni, and T. De
Bie, “Automatic Chord Estimation from Audio: A
Review of the State of the Art,” IEEE/ACM TASLP,
22(2), pp. 556–575, 2014.
[12] J. F. Paiement, D. Eck, and S. Bengio, “Probabilistic Melodic Harmonization,” in Proc. Conf. Canadian
Soc. for Comp. Studies of Intelligence, pp. 218–229,
2006.
[13] D. Morris, I. Simon, and S. Basu, “MySong:
Automatic Accompaniment Generation for Vocal
Melodies,” in Proc. CHI, pp. 725–734, 2008.
[14] S. A. Raczynski, S. Fukayama, and E. Vincent,
“Melody Harmonization with Interpolated Probabilistic Models,” J. New Music Res., 42(3), pp. 223–235,
2013.
[15] R. Groves, “Automatic Harmonization Using a Hidden Semi-Markov Model,” in Proc. AIIDE, pp. 48–54,
2013.
[16] M. Allan and C. Williams, “Harmonising Chorales by
Probabilistic Inference,” in Proc. NIPS, pp. 25–32,
2005.
[17] R. P. Whorley and D. Conklin, “Music Generation
from Statistical Models of Harmony,” J. New Music
Res., 45(2), pp. 160–183, 2016.
[18] T. Winograd, “Linguistics and the Computer Analysis
of Tonal Harmony,” J. Music Theory, 12(1), pp. 2–
49, 1968.
[19] M. J. Steedman, “A Generative Grammar for Jazz
Chord Sequences,” Music Perception, 2(1), pp. 52–
77, 1984.
[20] C. Raphael and J. Stoddard, “Functional Harmonic
Analysis Using Probabilistic Models,” Comp. Music
J., 28(3), pp. 45–52, 2004.
[21] D. Temperley, “A Unified Probabilistic Model for
Polyphonic Music Analysis,” J. New Music Res.,
38(1), pp. 3–18, 2009.
[22] H. Kaneko, D. Kawakami, and S. Sagayama, “Functional Harmony Annotation Database for Statistical
Music Analysis,” in Proc. ISMIR, 2010.
21

[40] F. Lerdahl and R. Jackendoff, A Generative Theory
of Tonal Music, MIT Press, Cambridge, 1983.
[41] N. Chomsky, Syntactic Structures, Mouton & Co.,
1957.
[42] M. Rohrmeier and I. Cross, “Statistical Properties of
Tonal Harmony in Bach’s Chorales,” in Proc. ICMPC,
pp. 619–627, 2008.
[43] N. Jacoby, N. Tishby, and D. Tymoczko, “An Information Theoretic Approach to Chord Categorization
and Functional Harmony,” J. New Music Res., 44(3),
pp. 219–244, 2015.
[44] P. Mavromatis, “Minimum Description Length Modelling of Musical Structure,” J. Math. Music, 3(3),
pp. 117–136, 2009.
[45] P. Mavromatis, “Exploring the Rhythm of the Palestrina Style: A Case Study in Probabilistic Grammar
Induction,” J. Music Theory, 56(2), pp. 169–223,
2012.
[46] S. A. Abdallah and N. E. Gold, “Comparing Models
of Symbolic Music Using Probabilistic Grammars And
Probabilistic Programming,” in Proc. ICMC/SMC,
pp. 1524–1531, 2014.
[47] K. Lari and S. J. Young, “The Estimation of Stochastic Context-Free Grammars Using the Inside-Outside
Algorithm,” Computer speech & language, 4(1),
pp. 35–56, 1990.
[48] J. A. Burgoyne, J. Wild, and I. Fujinaga, “An Expert
Ground Truth Set for Audio Chord Recognition and
Music Analysis,” in Proc. ISMIR, pp. 633–638, 2011.
[49] H. Tsushima, E. Nakamura, K. Itoyama, and
K. Yoshii, “Function- and Rhythm-Aware Melody
Harmonization Based on Tree-Structured Parsing and
Split-Merge Sampling of Chord Sequences,” in Proc.
ISMIR, pp. 502–508, 2017.

22

