Recursive Bayesian Networks: Generalising and Unifying Probabilistic Context-Free Grammars and Dynamic Bayesian Networks

Probabilistic context-free grammars (PCFGs) and dynamic Bayesian networks (DBNs) are widely used sequence models with complementary strengths and limitations. While PCFGs allow for nested hierarchical dependencies (tree structures), their latent variables (non-terminal symbols) have to be discrete. In contrast, DBNs allow for continuous latent variables, but the dependencies are strictly sequential (chain structure). Therefore, neither can be applied if the latent variables are assumed to be continuous and also to have a nested hierarchical dependency structure. In this paper, we present Recursive Bayesian Networks (RBNs), which generalise and unify PCFGs and DBNs, combining their strengths and containing both as special cases. RBNs define a joint distribution over tree-structured Bayesian networks with discrete or continuous latent variables. The main challenge lies in performing joint inference over the exponential number of possible structures and the continuous variables. We provide two solutions: 1) For arbitrary RBNs, we generalise inside and outside probabilities from PCFGs to the mixed discrete-continuous case, which allows for maximum posterior estimates of the continuous latent variables via gradient descent, while marginalising over network structures. 2) For Gaussian RBNs, we additionally derive an analytic approximation, allowing for robust parameter optimisation and Bayesian inference. The capacity and diverse applications of RBNs are illustrated on two examples: In a quantitative evaluation on synthetic data, we demonstrate and discuss the advantage of RBNs for segmentation and tree induction from noisy sequences, compared to change point detection and hierarchical clustering. In an application to musical data, we approach the unsolved problem of hierarchical music analysis from the raw note level and compare our results to expert annotations.


Introduction
Long-term dependencies with a nested hierarchical structure are one of the major challenges in modelling sequential data. This type of dependencies is common in many domains, such as natural In contrast, DBNs are sequential models with a fixed set of random variables that reoccur at each time step [15,16]. The variables at each time step may be discrete or continuous, latent or observed, and may have an arbitrary non-cyclic dependency structure among each other, with additional links from the previous and to the next time slice. They comprise important model classes as special cases, such as hidden Markov models (HMMs) if there is only a single discrete latent variable or linear dynamical systems if all dependencies are linear Gaussians [17,18,16]. However, DBNs only allow for a fixed chain of Markov dependencies between time slices and cannot represent nested hierarchical structures.
In this paper, we present Recursive Bayesian Networks (RBNs), a novel class of probabilistic models that combines the strengths of PCFGs and DBNs by allowing for nested hierarchical dependencies in combination with arbitrary discrete or continuous random variables ( Figure 1). Our main contributions are as follows: 1. With RBNs, we provide a unified theoretical framework for a large class of important sequence models, including PCFGs and DBNs. 2. We generalise inside and outside probabilities from PCFGs to continuous latent variables, allowing for maximum posterior (MAP) inference in arbitrary RBNs. 3. For Gaussian RBNs, we derive an analytic approximation for the marginal likelihood and marginal posterior distribution, allowing for robust parameter optimisation and Bayesian inference. 4. We provide a quantitative evaluation on synthetic data and an application to the challenging task of hierarchical music analysis.

Related Work
PCFGs have a long tradition for modelling nested hierarchical dependencies in symbolic data with a variety of parsing algorithms for inferring the structure and variables' values [13,6,14]. Beyond their application to sequential data, PCFGs have been generalised to graph structures [19,20,21,22], which readily transfers to applications of RBNs. Latent vector grammars (LVeGs) [23] are an extension of latent variable grammars (LVGs) [24,25] with continuous latent states. As for RBNs, approximate parsing is possible in the Gaussian case. However, both LVGs and LVeGs are special cases of RBNs and do not draw the connection to graphical models. More recently, the availability of automatic differentiation libraries, such as PyTorch [26], has lead to a number of applications where gradients are propagated through the entire parsing process [27,28,29,30].
x z x x y z=N z=T Figure 2: RBN in Chomsky normal form. The non-terminal transition p N (x , x | x) and terminal transition p T (y | x) are grouped into an RBN cell with a structural distribution p S (z | x). We use gates [58] to describe structural distributions and extended factor graph notation [black squares; 59] for conditional joint distributions. Considering all possible ways how n observations can be generated by recursively applying the RBN cell produces an RBN chart, as shown in Figure 3.
Orange and red boxes, respectively, indicate the subsets X 0:3 and Y 0:3 of latent non-terminal and observed terminal variables generated from x 0:3 .
The process of parsing a PCFG or RBN can be formally rewritten as a sum-product network [SPN; 31,32,33]. Factor graph grammars [FGGs; 34] generalise PCFGs, case-factor diagrams [35] and SPNs by using a hyperedge replacement graph grammar [19] to describe a distribution over graph structures that is more general than that of RBNs (not only trees). However, none of the approaches addresses the problem of inference with continuous variables that we are facing in RBNs (exponentially many terms with exponentially many nested integrals).
A wide range of probabilistic and neural models operate with a fixed graphical structure and are loosely related to RBNs. Hidden tree Markov models [36,37] generalise HMMs from chain to fixed tree structures. They model data at each node as observations of a latent Markov process on the underlying tree, which is part of the input data. Additionally estimating the underlying tree structure has been addressed in [38,39,40]. Recursive neural tensor networks [41] use a PCFG for parsing a given sequence of symbols to obtain a tree structure, which is then fixed and used as the backbone for a neural network. More generally, there is a number of methods for inferring a fixed structure for graphical models [42,16,43,44], SPNs [45,46,47,48,32,49], or graph neural networks [50,51,52]. All these methods have in common that a fixed structure is either given or estimated but not treated in a probabilistic Bayesian manner.
Some approaches attempt a Bayesian treatment of the unknown structure of a graphical model or SPN via dynamic programming [53,54] or Markov chain Monte-Carlo sampling [55,56,57]. However, the structure is assumed to be independent of the latent variables (they only become dependent conditional on the data) and the latent variables cannot be used to control the structure, as it is the case in RBNs. The challenge of continuous variables also remains unsolved.

Recursive Bayesian Networks
RBNs are template-based graphical models that define a joint distribution over network structures and variables' values. The number of template variables is fixed, but the number of instantiated variables, their connectivity and values are governed by the joint distribution. As a rough analogy, RBNs can be thought of as DBNs that can not only be connected linearly to form a chain but also hierarchically to form a tree structure. Alternatively, they can be thought of as a PCFG in which each symbol is a (possibly continuous) random variable.

Definition
RBNs have three types of template variables: 1) latent non-terminal variables (discrete or continuous), 2) observed terminal variables (discrete or continuous), and 3) latent structural variables (always discrete). In the simplest case, illustrated in Figure 2, an RBN has one template variable of each type. Formally, an RBN is defined as follows: Definition 1 (Recursive Bayesian Network). An RBN is a tuple (X , Y, Z, T , S, p P ) with X : a set of latent non-terminal template variables (1) Y : a set of observed terminal template variables (2) Z : a set of latent structural template variables, paired up with the non-terminal variables (3) T : a set of transition distributions p(v 1 , . . . , v η | x) from a single non-terminal variable x ∈ X to a set of non-terminal and/or terminal variables v 1 , . . . , v η ∈ X ∪ Y (4) S : a set of structural distributions p(z | x), one for each non-terminal/structural pair (5) p P : a prior/start distribution for exactly one non-terminal variable.
The cardinality of a structural variable z ∈ Z corresponds to the number of possible transitions from the associated non-terminal variable x ∈ X ; η in (4) is called the arity of the transition.
Generating with an RBN is straightforward. We start by sampling the value of the first non-terminal variable x from the prior distribution p P (x) and then repeat the following steps until no unprocessed non-terminal variables are left: 1. sample the value of the associated structural variable from p(z | x) 2. choose a transition distribution p(v 1 , . . . , v η | x) based on the structural variable's value 3. sample the variables v 1 , . . . , v η from the transition distribution 4. for all newly generated non-terminal variables, go to step 1.
The major challenge and focus of this paper is to perform joint inference over the latent structure and non-terminal variables' values conditional on a given set of observations.
Chomsky Normal Form: In the simplest non-trivial case, an RBN has one latent non-terminal, one observed terminal, and one latent structural template variable, with one non-terminal transition of arity η = 2 and one terminal transition of arity η = 1, as illustrated in Figures 2 and 3. It is defined by four distributions In analogy to PCFGs, we call this the Chomsky normal form (CNF). Any RBN may be rewritten in CNF (see Appendix A.1 for details).
RBN Chart: During inference, we will make use of an RBN chart, similar to the parse chart for PCFGs [6]. Each non-terminal variable is associated to a layer in the chart. For discrete variables, they store the actual distributions, while for continuous variables they either hold the point estimate (for MAP inference) or the parameters of the approximate distributions (for inference in Gaussian RBN). Different instances of the same template variable are identified by a subscript indicating the span of data generated from them, which also corresponds to their position in the chart (see Figure 3). Sets of variables that are generated from a specific latent non-terminal variable x i:k are denoted by a bold capital letter with a corresponding subscript (X i:k , Y i:k , Z i:k ); omitting the subscript refers to all variables (X, Y, Z); for X and Z this also includes the root variables x 0:n and z 0:n , respectively. The subscripts are to be interpreted as time intervals, that is, Y i:i is empty, Y 0:1 = y 1 is the first observation, Y n−2:n = (y n−1 , y n ) are the last two observations etc.
Comparison to PCFGs: Any PCFG can be rewritten as an RBN in two different ways, which we call abstraction and expansion (see Appendix A.2 for details). Abstraction of a PCFG produces a discrete RBN with one latent non-terminal and one observed terminal variable. The resulting RBN is exactly equivalent to the original PCFG but describes the same relations in a more abstract and compact way. In contrast, expansion of a PCFG considers the symbols of the grammar as random variables in their own right, thereby endowing them with additional (possibly continuous) degrees of freedom. The resulting RBN is therefore more powerful than the original PCFG. A PCFG is abstracted to a discrete RBN by defining the start/prior, transition, and structural distributions (7-10) as where S is the grammar's start symbol, A, B, C are non-terminal symbols, b is a terminal symbol, X is any right-hand side of a rule, z=N and z=T indicate a non-terminal and terminal transition, respectively, W · → · is the weight of the corresponding rule, and rules that do not exist in the original PCFG are taken to have zero weight. In expansion, the PCFG is only used to define a "skeleton" for the RBN, while the specific random variables and the concrete transition distributions need to be additionally specified. This means that the resulting RBN model is more powerful than the original PCFG, as the symbols may, for instance, be expanded to continuous random variables.

Inference
The two main goals of inference in RBNs are to 1) train model parameters by maximising the marginal data likelihood and to 2) compute posterior distributions or maximum posterior (MAP) estimates of the network structure and non-terminal variables. In PCFGs, both is achieved by computing inside and outside probabilities [14], which will be the starting point for our generalisation to continuous variables.
Inside and Outside Probabilities: We define inside and outside probabilities, β and α, for RBNs in analogy to how they are defined for PCFGs, the only difference being that the variables may be continuous. We thus have where n is the length of the sequence and Y is fixed (and therefore omitted as argument on the left-hand side). That is, β(x i:k ) is the marginal likelihood of generating the sub-sequence Y i:k conditional on the respective non-terminal variable x i:k , while α(x i:k ) is the marginal likelihood of generating the two sub-sequences Y 0:i and Y k:n as well as the non-terminal variable x i:k . In both cases, β and α are functions of the corresponding non-terminal variable with the structure and the remaining variables being marginalised out. Based on the inside and outside probabilities, the marginal data likelihood and the marginal posterior distributions over non-terminal variables are p(Y) = β(x 0:n ) p P (x 0:n ) dx 0:n (17) and respectively. p(x i:k | Y) is an unnormalised probability distribution that specifies the probability of x i:k to exist via the normalisation constant p(x i:k | Y) dx i:k , while the normalised version corresponds to the marginal posterior distribution of x i:k for the case that it does exist.
Inside probabilities are recursively computed bottom-up. For an RBN in CNF we start with the base case (19) for single observations and then iterate (20) to the top of the RBN chart Outside probabilities are recursively computed top-down, while making use of the inside probabilities α(x 0:n ) = p P (x 0:n ) α(x j:k ) = p S (z j:l =N | x j:l ) p N (x j:k , x k:l | x j:l ) α(x j:l ) β(x k:l ) dx j:l dx k:l .
As for PCFGs, the two terms in (22) correspond to the possibility of x j:k being generated as the right or the left child, respectively. The main conceptual difference to PCFGs is that we treat the discrete structural part (marginalised out by the sums) separately from the potentially continuous variables (marginalised out by the integrals). For RBNs that are not in CNF, the equations have to be adapted accordingly (see Appendix A.3 for the general case).
Marginalisation: Computing the marginal data likelihood (17) and the marginal posterior distributions over non-terminal variables (18) requires to solve an exponential (w.r.t. the length n of the sequence) number of nested integrals in (19)(20)(21)(22), which is generally intractable. However, for the special case of Gaussian RBNs, we provide an adaptive closed-form approximation in Section 2.3. Moreover, marginalising only over the network structure for a fixed assignment of the non-terminal variables X is straight forward and allows for maximum posterior (MAP) inference in general RBNs.
Maximum Posterior Inference: For a fixed assignment of all non-terminal variables X, we can compute the joint marginal likelihood p(X, Y) over observed terminal and latent non-terminal variables by only marginalising over the structure. This follows the same principle as above but uses the modified joint inside and outside probabilities β i:k := p(X i:k , Y i:k | x i:k ) (23) and α j:k := p(X 0:j , Y 0:j , x j:k , X k:n , Y k:n ) , (24) where all variables are fixed (and therefore omitted as arguments on the left-hand side). Analogously, the joint marginal likelihood and the marginal posterior probability of x i:k to exist then are p(X, Y) = β 0:n p P (x 0:n ) (25) and where p i:k is the probability of x i:k to exist for this specific assignment of X. The corresponding equations for the recursion differ from (19)(20)(21)(22) only in that they do not integrate out the latent nonterminal variables (see Appendix A.3.3). As before, all computations can be efficiently performed via dynamic programming. Gradients w.r.t. the variables and/or parameters are readily obtained from libraries such as PyTorch [26]. Optimising the values of the latent non-terminal variables X via gradient descent yields maximum posterior (MAP) estimates, while the structure is marginalised out. MAP estimates for the structure (i.e. the best tree) conditional on an assignment for X can be computed (as for PCFGs) by replacing summation with maximisation [13,60].
There are two caveats: First, due to marginalising over multiple (exponentially many) network structures, p(X, Y) may be highly non-convex and optimising X via gradient descent is not guaranteed to find the global optimum. This is even the case for purely Gaussian RBNs, for which p(X, Y) is a mixture of Gaussians (one for each structure). Second, we can optimise X while marginalising out the structure and we can optimise the structure for a fixed assignment of X. However, successively optimising X and the structure is not equivalent to jointly optimising both and the maximum of p(X, Y) may be unrelated to the maximum of the best structure (also see Figure 4). This means that generally, exact joint MAP inference over the latent variables and the structure is hard. For Gaussian RBNs, we provide an approximate solution below.

Gaussian RBNs
In a Gaussian RBN (GRBN), the prior, non-terminal, and terminal distributions are linear Gaussians and the termination probability (structural distribution) is constant [termination/structural] For clarity, we will show all derivations for GRBNs in this basic form. For our evaluations and the application to music, we use a slightly extended version that includes linear transformations, mixtures of Gaussians, and multi-terminal transitions (Section 2.3.1). The derivations do not fundamentally change for the extended case (see Appendix A.4). In Appendix B, we show all calculations on a simple example.
Adaptive Approximation: If the structure of a GRBN was fixed, all variables would be jointly Gaussian distributed as in a conventional Gaussian Bayesian network [16]. However, due to the mixture max approximation components Figure 4: Three Gaussians components, the resulting mixture (green), maximum (orange), and moment-matching single Gaussian approximation (blue). Note that the maximum of the mixture (), the best component (), and the approximation () may be unrelated.
x λ W x x τ y n n z=N z=T z Figure 5: Graphical model of the Gaussian RBN for modelling music. The additional transposition variable τ is marginalised out during inference; the number of jointly generated observations n is uniquely determined by the location in the parse chart.
unknown structure, we effectively have a mixture of exponentially many Gaussians, one for each possible structure. While in principle all integrals can be solved analytically, the exponential growth makes exact inference intractable. Therefore, our goal is to derive a parsing strategy that retains tractability by adaptively applying local approximations to the Gaussian mixtures occurring in each recursion step. We will here focus on the simplest case of approximating the mixtures with a single Gaussian (illustrated in Figure 4, details in Appendix A.4.2), which can be efficiently computed in closed form [18,61]. The inside and outside probabilities are thus represented by a simple Gaussian and this form is reestablished in each iteration by approximating the occurring mixtures. Consequently, the marginal posterior distributions over latent variables (18) are also simple Gaussians and the marginal data likelihood (17) can be computed in closed form. This approximation scheme can be extended and refined by using existing methods for approximating each Gaussian mixture by one with fewer components [62,63].
Marginalisation: In (20) and (22), we have to integrate over products of Gaussian distributions to marginalise out the latent variables. To solve these integrals, we make use of the fact that the product of two Gaussians over a variable x can be rewritten as [see e.g. 64] withΣ Hence, when integrating over x, only the first term on the rhs. of (33) remains. A detailed step-by-step derivation of all results can be found in Appendix A.4.1. With the latent variables being marginalised out, (20) and (22) become simple mixtures of Gaussians that can be easily approximated to retain the simple analytic form of the inside and outside probabilities.
Tree Induction: As described above, exact joint MAP inference over the continuous latent variables and the structure is generally intractable. Moreover, the maximum of the approximate posterior does not necessarily coincide with the maximum of the exact posterior or that of a particular structure (see Figure 4). Thus, first optimising X (based on the approximation) and then estimating the structure (conditional on the picked value of X) may lead to arbitrarily bad results for tree induction. Therefore, we leverage the adaptive character of our approximation scheme to compute local structure estimates in each step, before loosing relevant information due to further approximations. Specifically, during the bottom-up pass for computing inside probabilities, all structures are scored by the maximum of their marginal likelihood, based on its current approximation (31). The best overall structure is then selected (as usual) in a top-down pass (see Appendix A.4.3 and our example in Appendix B).

Gaussian RBNs for Music
For the application to music, we slightly extend the basic GRBN discussed so far by introducing transpositions and multi-terminal transitions (changes in the equations highlighted in blue). The corresponding graphical model of the RBN cell is shown in Figure 5. Furthermore, we describe how GRBNs can be applied to categorical data.
Transpositions: A transposition rotates the dimensions of the latent variable by a number of steps τ before generating the child. This is achieved by multiplying with an orthonormal transposition matrix T τ that corresponds to the identity matrix with cyclicly rearranged columns. For the prior distribution, we assume a uniform weighting of all possible transpositions where D is the dimensionality of the data (D = 12 for music in 12-tone equal temperament). For the non-terminal transitions, the probability for a specific transposition is determined by the weight parameter W Note that transpositions are only applied to the left child, because Western classical music is thought to be fundamentally goal directed [2,8,65]. This means that the character of a section is largely determined by how it ends (the right child), which should also be reflected in the value of the parent node. In contrast, the role of the left child is to harmonically prepare the ending (or prepare a preparation to the ending etc). We therefore allow for arbitrary transpositions in the left child and we will see below that our model indeed captures the most important type of preparation in Western classical music: the cadential dominant-tonic progression.

Multi-Terminal Transitions:
A multi-terminal transition generates multiple observed variables from a single latent variable. The variables are generated i.i.d. and their number is governed by a Poisson distribution with rate parameter λ Multi-terminal transitions do not conform to the CNF assumed so far and we need to add the term [for GRBNs, see (30)] (39) to (20) in order to account for the possibility to terminate from a higher-level variable. For k = i + 1, this term becomes the base case (19) of an RBN in CNF.
Multi-terminal transitions account for the situation where changes in the hierarchical structure occur at a lower rate than the time series is sampled. In between the structural changes, the data is assumed to be generated from the same model, which could also be more elaborate than i.i.d. samples, as long as the relevant model parameters are captured by the RBN's latent variables.
Categorical Data: The observed variables of a GRBN are unconstrained real-valued, which poses a problem if the data are categorical. This situation is comparable to using Gaussian processes (GPs) [66] for classification and can be approached with similar methods. In our application to musical data, we observe one or more notes being played at any particular time and normalise these counts to obtain observations that correspond to the parameter of a categorical distribution. The natural likelihood function for this type of observations is a Dirichlet distribution. Therefore, we adapt the approach suggested in [67] for GPs, who assume a Dirichlet likelihood, which is then approximated by a Gaussian likelihood in log-space. Since an observation from a Dirichlet distribution corresponds to a normalised sample from independent Gamma distributions, each Gamma distribution can be separately approximated by a log-normal distribution, which results in a diagonal covariance matrix for the Gaussian likelihood in log-space. Matching the first and second moment yields [67] y (l) where 0 < y (l) j < 1 is the l th element (normalised count) of the j th observation, y j is the corresponding mean of the approximate Gaussian likelihood in log-space, and Σ (j) ll is the l th element on the diagonal of the covariance matrix for the j th observation. We thus have to replace y j and Σ (j) for y j and Σ T in (37).

Experiments
We performed a quantitative evaluation on synthetic data and applied our model to hierarchical music analysis of Bach preludes. We show that RBNs are superior to change point detection (CPD) and hierarchical clustering (HC) for tree induction and our method is able to infer fundamental harmonic principles of Western classical music. Experiments were run on a 3.6 GHz Quad-Core Intel Core i7 processor with 32GB RAM. The model parameters were trained via gradient descent on the (approximate) marginal neg-log likelihood.

Quantitative Evaluation on Tree Induction
We performed a quantitative evaluation on synthetic data for the task of segmenting a noisy time series and inferring the underlying tree. For comparison, we used the best-performing change point detection (CPD) method from the ruptures library [72] for segmenting the time series, combined with bottom-up hierarchical clustering (HC) for inferring the tree structure ("HC/CPD"). For details of the methodology, see Appendix C.1.
The evaluation results in Figure 6(a) show that the RBN tree estimates (Section 2.3) consistently outperform the one from HC/CPD, in terms of both precision and recall (and thus also in F1 measure). The marginal node probabilities show an interesting performance pattern. They excel in terms of precision, which means that a node with high marginal probability is very likely to actually exist in the tree (low false-positive rate). However, they severely underestimate the overall node probabilities, which leads to recall falling far below the baseline. This means that a node with low marginal probability may in fact occur in the tree (high false-negative rate).
We think that the poor recall measure of the marginal probabilities is primarily due to (and the downside of) a fully Bayesian treatment that quantifies uncertainty. Even if the marginal probabilities have a maximum at the correct node location, probability mass will still spread around it and be allocated to a number of less probable locations. While this is the desired behaviour of a Bayesian method, it inevitably results in a lower recall value. The high precision value confirms that uncertainty is adequately quantified and not underestimated. That being said, the marginal probabilities provide an exceptionally rich basis for qualitative analyses. For instance, all ground-truth nodes are located at local maxima of the marginal probabilities and we can read off a number of other potential node locations, which essentially trace out the grid defined by the piece-wise constant segments (see Figure 8 in Appendix C.1).

Hierarchical Music Analysis
Harmonies in Western classical music exhibit a nested hierarchical structure that can be modeled by PCFGs operating on abstract chord symbols [8,73,10,3]. While these grammars can be applied to expert annotations of a musical score, hierarchical music analysis from the raw note level is an unsolved problem. We trained a GRBN (Section 2.3.1) on the 24 major preludes of Johann Sebastian Bach's "Wohltemperiertes Klavier I & II" (see Appendix C.2 for technical details and complete results).
Our first major finding is that the prior mean, shown in Figure 6(b), corresponds to a major pitch profile (as could be expected from the training data) and is in excellent agreement with recent Bayesian estimates from the literature [68]. The fact that the major profile appears in the prior (i.e. as the continuous equivalent of a grammar's start symbol) shows that our model picks up fundamentally important structures from the musical data. Our second finding is that only two transpositions have non-zero weights: the identity with a weight of 78% and the fifth scale degree (7 semitones) with a weight of 22%. This corresponds to the left child being generated as the dominant of the parent and realises the most important harmonic preparation in Western classical music: the cadential dominant-tonic relation. A closer inspection of the expert analysis ( Figure 10 in Appendix C.2) reveals that when considering the possible surface patterns (raw notes) of the labeled chords, most non-identity transitions can indeed be explained as (noisy) fifth transpositions. The strong weight of fifth transpositions in our model is a highly non-trivial empirical confirmation of the established music theoretical insight that Baroque music is fundamentally driven by dominant-tonic relations. While the estimated tree in Figure 6(c) fails to reproduce the large-scale structure of the expert analysis (e.g. the separation into two main parts), it accurately captures the measure-wise harmonic changes on the bottom level.
On the one hand, we see considerable room for improvement by integrating more advanced concepts, such as different modes (major/minor), diatonic in addition to chromatic transposition, or balancing of trees. On the other hand, our model was able to capture fundamental properties of Western classical music based on only 24 pieces. We therefore think that Gaussian RBNs are a highly promising approach for hierarchical music analysis from the raw note level, which should be further investigated.

Conclusion
We introduced Recursive Bayesian Networks (RBNs), a novel class of probabilistic models that unifies the strengths of probabilistic context-free grammars (PCFGs) and dynamic Bayesian networks (DBNs), generalising both model classes. We defined RBNs as a joint distribution over tree-structured Bayesian networks and their (discrete or continuous) variables and described how to perform inference over both the model structure and the variables by leveraging parsing methods for PCFGs. The provided formalisation connects with the methods for formal grammar as well as with the versatile notation for graphical models. On two data sets, we demonstrated the potential of RBNs for modelling nested hierarchical dependencies in real-valued time series and musical data. The class of RBNs represents a substantial contribution to the machine learning toolkit by unifying two of the most important approaches for modelling sequential data and bears a large potential for further development and applications.

Appendix
Recursive Bayesian Networks: Generalising and Unifying Probabilistic Context-Free Grammars and Dynamic Bayesian Networks Transforming an RBN into CNF is done analogously to the procedure for PCFGs. We assume the original RBN does not contain any epsilon productions, that is, a non-terminal variable always produces one or more other non-terminal and/or terminal variables.

1) Eliminate terminal variables from mixed transitions: This is done by introducing intermediate non-terminal variables. For each transition
with non-terminal variables x (1) , x (2) , . . . ∈ X and terminal variables y (1) , y (2) , . . . ∈ Y, we introduce new non-terminal variables x y (1) , x y (2) , . . . and replace the transition by where the new non-terminals x y (i) replace the original terminals y (i) . We then add new deterministic transitions . . . that convert each non-terminal to its equivalent terminal variable. For the newly added non-terminals, there is only a single transition and hence a degenerate structural variable that can only take a single value.
2) Eliminate more than two latent non-terminal variables: This is done by introducing new nonterminals that capture combinations of multiple old non-terminals. Below, we show how the number of non-terminals can be reduced by one. Applying this procedure repeatedly allows for reducing the number of non-terminals from an arbitrary number down to two, as required for CNF. A transition that generates n non-terminals x (1) , . . . , x (n) is rewritten as where we introduced the new non-terminal variable x = (x (1) , x (2) , . . . , x (n−1) ) that stores all the information from the first n − 1 original non-terminals. The actual "work" is done by p(x , x (n) | x), which is the equivalent of the original n-fold transition. p(x (1) , x (2) , . . . , x (n−1) | x ) is a deterministic transition that just "unpacks" the information stored in x . Repeating this procedure to come to only pairwise transitions corresponds to a chain of these deterministic "unpacking" operations. As above, the newly added non-terminals have only a single possible transition.
3) Eliminate unary cycles: Unary cycles p cycle (x | x), where x and x are the same non-terminal template variable (x ≡ x ), are first transformed into unary transitions to a new non-terminal variable and then eliminated as described below. We define a new non-terminal variablex = (x , n), where n > 0 represents the number of steps taken in the cycle before exiting it and x is the value at the moment of exiting it. The transition distribution tox is where in the recursive case, the variables of all intermediate steps are successively marginalised out. In practical applications, if p(z=cycle | x ) < 1, the probability of remaining in the cycle decays exponentially and the recursion can be truncated after a number of steps. If on the other hand p(z=cycle | x ) ≈ 1 so that truncating is not possible, one can work with the stationary distribution of the resulting Markov chain (i.e. the Markov chain with transition distribution p cycle (x | x)).
The structural probability to take a transition from x tox is p(z=cycle | x), i.e. the probability of entering the cycle in the first place. The RBN cell ofx is identical to that of x, except for the transition into the cycle, which is eliminated (the structural distribution thus has to be renormalised for the remaining transitions). The transitions use only the x -component ofx, ignoring the n-component. In this way, we have expresses the state after an arbitrary number of steps in the unary cycle as a distinct value of the new non-terminal variablex.
The intermediate variable x and its RBN cell is eliminated if it was only reachable via x. The new transitions p 1 , . . . , p n are merged into the cell of x, while the original transition p unary to x is removed. This requires redefining the structural distribution p(z | x) such that the probability mass p(z=unary | x) that was formerly assigned to p unary is now split among the new transitions p 1 to p n according to the structural distribution p(z | x ) of x . Specifically, for a new transition p i , we define A.2 Relation to PCFGs As described in Section 2.1, a PCFG can be rewritten as an RBN by abstraction or expansion, where abstraction produces an equivalent RBN that describes the same relations in a more abstract and compact way, while expansion produces a more general RBN using the original PCFG as a skeleton.
We describe the two procedures in detail below and use the following definition of a PCFG: where A, B, C ∈ N are non-terminal symbols of the PCFG, b ∈ T is a terminal symbol, X ∈ N 2 ∪T is any right-hand side of a rule, z=N and z=T indicate a non-terminal and terminal transition, respectively, W A→X is the weight of the corresponding PCFG rule, and rules that do not exist in the original PCFG are taken to have zero weight.
To show equivalence, we need to prove that the transition probabilities from a given non-terminal symbol A are the same in the original PCFG and the new RBN.
Proof. In the RBN, the probability for a non-terminal transition A → B C is and that for a terminal transition A → b is which matches the corresponding probabilities in the PCFG, gained by normalising the respective weights.
Conversely, any discrete RBN can be rewritten as a PCFG. Theorem 2. A discrete RBN with n latent non-terminal template variables x 1 , . . . , x n , m observed terminal template variables y 1 , . . . , y m , and a prior p P (x 1 ) over x 1 can be rewritten as a PCFG with where · ⊕ · concatenates the value ranges of the respective variables, X ∈ x i denote that the value X is in the value range of the RBN variable x i , and the second case in (60) requires there be a transition p i (x 1 , . . . , x k | x i ) such that A ∈ x i and X ∈ x 1 ⊕ · · · ⊕ x k .

A.2.2 Expansion of a PCFG
Expansion of a PCFG to an RBN uses the PCFG as a "skeleton" to define the number of template variables and the structural transitions. The domains and transitions for the variables need to be added, which results in an RBN that is more powerful than the original PCFG. Specifically, we have for the sets of latent non-terminal and observed terminal template variables and for the structural transitions. Additional, we have to define the domain for each of the non-terminal and terminal variables in X and Y, and for each rule A → X 1 X 2 . . . from the original PCFG, we have to define a concrete transition distribution . . ∈ X ∪ Y are non-terminal or terminal variables in the RBN, respectively, depending on whether X 1 , X 2 , . . . ∈ N ∪ T are non-terminal or terminal symbols in the PCFG).
Expansion of a PCFG into an RBN seems appealing if a simple PCFG can be used to describe the type of variables (as opposed to their values) as well as the structure of the generative process. The actual transitions on the variables' values may then take place on a sub-symbolic/continuous level, which cannot be described by a PCFG.

A.3.1 Inside Probabilities
The inside probability is the probability of generating the observed terminal variables Y i:k from the latent non-terminal variable x i:k . This means that we need to marginalise over all possible paths of generation. Transitions may directly generate observed variables, but they may also generate lower-level non-terminals, in which case we have to recurse using the respective inside probabilities from those variables.
Let T x ⊆ T be the set of possible transitions from the latent non-terminal template variable x ∈ X (of which x i:k is one specific instantiation), with p(z i:k =τ | x i:k ) being the probability for the transition τ ∈ T x to be selected. This constitutes the first sum in (64) below, which marginalises over the transitions.
The transition τ generates η new non-terminal and/or terminal variables, where η is the arity of τ . These may be located at different positions in the parse chart, depending on which part of the observed variables Y i:k is generated from them. That is, the variables' locations in the parse chart are not known during generation and are determined in hindsight once all observed variables are generated; thus, they are known for parsing. We denote the respective splitting points by j 1 , . . . , j η−1 (they have to fulfill the condition i < j 1 < . . . < j η−1 < k) and the respective variables by v i:j1 , . . . , v jη−1:k ∈ X ∪ Y. The second multi-sum in (64) aggregates the probabilities of the different splitting possibilities, that is, of all valid assignments of j 1 , . . . , j η−1 (η − 1 degrees of freedom). For instance, a transition of arity η = 2 has one free splitting point j 1 to sum over.
Some of the generated variables may be observed/terminal variables, for which nothing more needs to be done as they directly constitute the respective part of Y i:k . For the subset of non-terminal variables, which we denote by {v j:j ∈ X }, we need to insert their respective inside probabilities and marginalise them out. This constitutes the product and multi-integral in (64).
The general form of the inside probabilities then is The concrete RBNs considered in the paper have only two transitions, one non-terminal transition of arity two and one terminal transition of arity one (for CNF) or more (for the extended GRBNs used in the quantitative evaluation and for modelling music). For non-terminal transition of arity two, the multi-sum in (64) reduces to a single sum and the multi-integral to a double integral, which gives us (20). For the terminal transition, (64) simplifies to (19) or the extended version (38), respectively.

A.3.2 Outside Probabilities
The outside probability is the joint probability of generating the latent non-terminal variable x j:j as well as the prefix and suffix of observed terminal variables, Y 0:j and Y j :n , respectively. For this, we now have to consider all possible ways how x j:j as well as the prefix and suffix could have been generated from a parent non-terminalx (x andx may correspond to the same template variable or to two different ones).
Let T −1 x ⊆ T denote the set of transitions that include x as one of the generated variables. Importantly, if x appears multiple times in the generated variables of a particular transition, these different options of generating x are represented as multiple distinct entries in T −1 x , one for each occurrence. The first sum in (66) runs over these different possibilities of generating x.
For a transition τ ∈ T −1 x of arity η, let j 0 , . . . , j η be the splitting points, including the start and end point j 0 and j η of the parent variablex j0:jη , which have to fulfill the condition 0 ≤ j 0 < . . . < j η ≤ n (where n is the length of the sequence). One pair of adjacent splitting points (j m , j m+1 ) corresponds to the occurrence of x j:j , where m is the position (starting at zero) at which x appears in the generated variables of the particular transition τ . We therefore have the additional constraints j m = j and j m+1 = j , resulting in η − 1 remaining free indices to sum over (as for the inside probabilities above). This corresponds to the second multi-sum in (66).
The set of non-terminal variables generated from the parentx j0:jη , excluding x j:j , is denoted by {v l:l ∈ X } \ x j:j . Together with the directly generated terminal variables, these generate part of the prefix and suffix, Y j0:j and Y j :jη . The remaining prefix and suffix, Y 0:j0 and Y jη:n , are generated from the parent variablex j0:jη . For the parent, we recurse via its outside probability α(x j0:jη ), while for the newly generated non-terminal variables (except x j:j ), we have to use the respective inside probability β(v l:l ) in (66). Additionally, we have to marginalise out the parent (first integral) and the newly generated non-terminal variables (second multi-integral). For a non-terminal transition of arity two, as we have it in the paper, the multi-sum in (66) reduces to a single sum and {v l:l ∈ X } \ x j:j contains only a single non-terminal, the second child. Importantly, T −1 x has two elements, one for x j:j being generated as the right child and one for it being generated as the left child, which gives us (22).

A.3.3 Joint Inside and Outside Probabilities
The joint inside and outside probabilities (23) and (24) for an RBN in CNF are computed analogously to (19)(20)(21)(22) for the normal inside and outside probabilities, that is, This differs from (19)(20)(21)(22) only by dropping the integrals and dependencies on the non-terminal variables (as their values are now fixed). Joint inside and outside probabilities for the general case are obtained from (64) and (66) analogously, i.e. again by dropping the integrals and dependencies on the non-terminal variables.

A.4 Gaussian RBNs
In the following, we present derivations for the extended case of GRBNs, described in Section 2.3.1, which includes linear transformations T for the left child. For this, we will make use of the fact that a normal distribution over a transformed variable T x can be rewritten as where ||T || is the absolute value of the determinant of T and in (72) we made use of the fact that in our case, the transformation matrices are orthonormal, so that T −1 = T and ||T || = 1.
Note that for an implementation, some of the results should be rewritten in order to minimise the number of matrix inverses that need to be taken. In particular, the identity is useful for the implementation, but we omit it in our derivation for clarity.

A.4.1 Marginalisation
For the inside probability β(x i:k ), the integral in (20) is N (x i:j ; µ where in (74) we inserted (31) and (36); in (75) we used (33) twice to rewrite the pairwise products of Gaussians over x i:j and x j:k and marginalise them out; in (76) we used (72) to rewrite the transformation; and in (77) we used (33) a third time to rewrite the resulting product as a single Gaussian over x i:k .
For the outside probability α(x j:k ), the integrals in (22) for x j:k being generated as the right child are with and where in (81) we took the constant termination probability (30) out of the integral and dropped it; in (82) we inserted (31), (32) and (36); in (83) we applied (33) to marginalise out x i:j ; in (84) we used (72) to rewrite the transformation; in (85) and (86) we used (33) twice to marginalise out x i:k ; and in (87) we rewrote the final result using (91-93). Due to the asymmetric terms in the outside probabilities, the result is somewhat more complex than for the inside probabilities.
Analogously, the integrals in (22) for x j:k being generated as the left child are p N (x j:k , x k:l | x j:l ) α(x j:l ) β(x k:l ) dx j:l dx k:l (94)

A.4.2 Approximation
A Gaussian mixture distribution p(x) with normalised mixture weights c i , means µ i , and covariance matrices Σ i can be approximated with a single Gaussian as The approximation p(x) matches the first and second moments of p(x) and minimises the Kullback- [61,18]. This direction of the KLD is the one used e.g. in expectation propagation, not the one used in e.g. variational methods [18]. That means, p(x) will adequately represent the support and uncertainty of p(x) (e.g. it will be non-zero wherever p(x) is non-zero). On the other hand, a value of x may have a high probability in p(x) even though in p(x) it has not (also see Figure 4).

A.4.3 Tree Induction
Exact joint optimisation of the structure and the continuous latent variables is intractable. We therefore choose the best tree for a GRBN based the maximum of the (approximated) inside probability. Inserting (30) and (77) x 0:4 x 0:3 x 3:4 x 0:1 x 1:3 x 1:2 x 2:3 y 1 y 2 y 3 y 4 Figure 7: Parse chart for a sequence of length n = 4 with best tree estimate (see text for details).
which is maximised by taking the mode of the Gaussian and maximising over j and (if using If we have multi-terminal transitions (or more generally other possible transitions), we also have to maximise over the different possible transitions. For each non-terminal variable, we compute and store the best choice during bottom-up computations of the inside probabilities. Afterwards, we can construct the best tree by starting at the root node and recursively picking the best structure top-down.

B Example
In this section, we present the complete calculations for the inside probabilities, the tree estimate, and the marginal likelihood for a basic GRBN (no transpositions or multi-terminal transitions) on a simple example sequence of length n = 4 (also see Figure 7). We assume parameters in (27)(28)(29)(30) and a scalar sequence Y = (y 1 , y 2 , y 3 , y 4 ) = (0, 1, 2, 0) .

B.1 Preliminaries
Inside probabilities are approximated with a single Gaussian specified by c i:k , which are the relevant quantities to be computed. At the bottom level, we use (19) for the base case and insert (29) and (30) to obtain where we can directly read off c where we already inserted Σ NL = Σ NR = 1.
If the sum in (113) has only a single term, we immediately get If there is more than one term in the sum in (113), this means that there are multiple splitting options that are marginalised out and we therefore need to do two things.
First, we have to identify the best splitting option to be able to compute the tree estimate. This is done by using (109) and comparing the values of where Σ Second, we have to approximate the resulting mixture with a single Gaussian using (107), where the mixture weights have to be normalised. For the univariate/scalar case considered here, we then get Finally, the marginal likelihood (17) where we have used (33) to rewrite the product of Gaussians and inserted µ P = 0 and Σ P = 1.
The normal distribution is defined as where the second line is for the scalar case.

B.2 Calculations
We start with the inside probabilities at the bottom level for the latent variables x 0:1 , x 1:2 , x 2:3 , x 3:4 and from (112) we read off (without any approximations) Next, we compute the inside probabilities on the first level for the variables x 0:2 , x 1:3 , x 2:4 . The only possible splitting option is for j = i + 1 and from (117) we get (again without approximation) Turning to the values for x 0:3 and x 1:4 , we now have two terms in the sum in (113), which means that we need to evaluate the best split and approximate the mixture. The corresponding parameters of the mixture are given by (114-116) as that for x 0:3 both splits are equally well and for x 0:3 the split x 1:4 → (x 1:3 , x 3:4 ) at j = 3 is better. This is intuitively clear, since generating (y 2 , y 3 ) = (1, 2) from the same non-terminal variable x 1:3 = 1.5 is more likely than generating (y 3 , y 4 ) = (2, 0) from x 2:4 = 1, given that in both cases the values are generated from a Gaussian with variance 1.
We approximate the mixtures with a single Gaussian with parameters given by (119-121) as and we see that the split x 0:4 → (x 0:3 , x 3:4 ) for j = 3 is the best one. Intuitively, this makes sense because it splits between y 3 and y 4 , which is the biggest step. Not we have only minor differences between the split options, because for simplicity we have chosen our variance parameters with a value of 1, which is relatively large compared to the spread of the values. Choosing smaller variances would result in more prominent splitting preferences.
We can now construct the full tree by also picking the best split for x 0:3 , which is a tie between splitting at j = 1 and j = 2, so we can choose either one (in practice one might consider random tie breaking to avoid biases due to variable order). The resulting tree is shown in Figure 7.
The parameters for the inside probability of x 0:4 , given by approximating the three Gaussian mixture components, are

C.1 Details for Quantitative Evaluation
We performed a quantitative evaluation on synthetic data for the task of segmenting a noisy time series and inferring the underlying tree. We used the Gaussian RBN for music (Section 2.3.1) with  some simplifications: 1) data were continuous and not categorical, 2) they had only three dimensions instead of twelve, 3) the prior distribution did not have any transpositions, 4) the left child could only be transposed by zero or one step. The data were sampled from this model, using zero prior mean, Σ p = 1, Σ NL = Σ NR = 1 · 0.1 2 , λ = 5, equal weights W 0 = W 1 = 0.5 for transposition by zero or one, and Σ T = 1 · noise 2 with different noise levels {0.01, 0.05, 0.1, 0.15, 0.2, 0.25}. We used a terminal probability of p term = 0.6 for sampling and rejected any sequences with a length outside the range of 50-55. An example of the data is shown in Figure 8.
For comparison, we used the best-performing change point detection (CPD) method from the ruptures library [72] for segmenting the time series, combined with bottom-up hierarchical clustering (HC) for inferring the tree structure ("HC/CPD"). For each noise level, we selected the CPD method and parameters with best F1 score based on the ground-truth segments of 100 training sequences (see Figure 9). In HC, pairs of adjacent segments with the smallest Euclidean/L2 distance between their mean values were successively combined to construct the tree. For the RBN, all parameters were trained from scratch by minimising the marginal likelihood of the observations of only 10 training sequences (i.e. no ground-truth information and less training data than for the baseline was used), separately for each noise level.
The models were evaluated on 500 test sequences by computing their precision and recall w.r.t. the ground-truth trees. Each possible node was treated as a separate binary classification task and the results reflect the number of correctly predicted nodes. Note that due to the strongly unbalanced class distribution (many more possible node locations than actual nodes in the tree) precision and recall or the combined F1 score are the appropriate performance metrics (as opposed to e.g. accuracy). For the RBN they were computed in two different ways: 1) based on the best-tree estimate ("RBN max") and 2) based on the marginal node probabilities (18) For the single-tree estimates (baseline model and best-tree estimate from RBNs) we compared the ground-truth and estimated tree node-by-node to count correctly predicted nodes (TP), nodes that are in the prediction but not the ground-truth (FP), and nodes that are in the ground-truth but not the prediction (FN). For the marginal node probabilities, we computed the corresponding rates by counting all nodes in the ground-truth tree (TP+FN), summing the marginal probabilities over all

C.2 Hierarchical Music Analysis
The scores were pre-processed by computing pitch-class distributions (PCDs), as used for the identification of musical keys [74,75,76], using the pitchscapes library [71]. We used a resolution of 70 equally spaced time slices per piece, resulting in sequences of 12-dimensional categorical distributions. The tree shown in Figure 6(c) for Johann Sebastian Bach's Prelude No. 1 in C major, BWV 846, corresponds to a harmonic expert analysis performed by the authors. Figure 10 shows the annotated tree with additional chord labels, which are provided in a simplified notation commonly used in Jazz lead sheets to be more accessible to a broad audience. In Table 1, we list the results for all 24 preludes. For a better interpretation of the model and the presented results, there are two relevant points to consider.

C.2.1 Chromatic versus Diatonic Transposition
It is interesting to look in more detail at what musical aspects the model can or cannot represent. In a nutshell, it can represent chromatic transposition but cannot represent diatonic transposition, which has a number of consequences, as described in the following.
The transpositions of the left child perform a cyclic rotation of the corresponding probabilities in the pitch-class distribution represented by the latent variable, which corresponds to a chromatic transposition. This determines not only which pitch classes have a significant probability to occur (the in-scale tones) but also the specific weights. For instance, the tonic and fifth scale degree typically have the highest weights. A transposition by 5 or 7 semitones from a current major key (say C major) thus corresponds to a modulation to the sub-dominant (F major) or dominant (G major) key, respectively. This includes adaptation of the fourth and seventh scale degree of the target key, respectively (B→B for F major; F→F for G major), as well as the correct assignment of strong weights to the tonic and fifth scale degree.
However, diatonic transposition cannot be represented in this way. For instance, to represent a modulation from C major to A minor, the model has two options that are both far from optimal. 1) It can choose not to apply a chromatic transposition, which ensures that all in-scale tones are correctly represented (i.e. they have significant weight). This, however, means that the relative weights are not appropriate for A minor. In particular, the strong weights on the tonic and fifth scale degree are not present and, instead, the third and seventh scale degree (C and G, the former tonic and fifth scale degree) have disproportionally strong weight. Correcting these weights has to occur through the Gaussian transitions, which can only be explained with a relatively high transition variance.
2) The second option would be to perform a chromatic transposition by 9 semitones, which ensures that the strongest weights remain on the tonic and fifth scale degree of the new key. However, three out-of-scale tones (C , F , G ) now have a high weight, while the respective in-scale tones do not. Again, this has to be corrected for by the Gaussian transition noise at a potentially even higher cost than in the first case. This is a highly plausible explanation for why we only see non-zero weights for the identity and (chromatic) transposition by a fifth in our experiments. Any diatonic modulations are best explained by reweighting using via Gaussian transition noise without a transposition, rather than by a chromatic transposition, which would require an even stronger reweighting (except for modulation to the sub-dominant and dominant key, which can be appropriately explained by a chromatic transposition).

C.2.2 Chord Labels
It is important to note that the chord labels in the expert annotation convey significantly more information than just what pitch classes can be expected to occur in the respective section. For example, the very same pitch-class distribution of G-C-E could amongst others be labeled as a C major chord in second inversion, a G major chord with 64-suspension (Gsus64), or an A minor seventh chord with omitted root, which might be easily confused by a musically untrained annotator. Which of these labels is correct depends in many cases on the context, such as how a chord resolves to the next one. While these differences are important from a musical perspective (they express a different experience of the same musical events), our model was trained to only predict pitch-class distributions. Therefore, in its current state, it cannot reproduce these distinctions, but we expect future versions to significantly improve in this respect. End of Table 1