Modelling birdsong transmission with methods from molecular sequence analysis

TL;DR

This paper introduces a novel Markov model inspired by molecular sequence analysis to study birdsong transmission, fitting it with advanced algorithms and analyzing Java sparrow songs, revealing limitations in predictive power due to biological factors.

Contribution

The paper develops a new Markov model for birdsong transmission based on molecular sequence analysis methods, applying a novel fitting algorithm to real bird song data.

Findings

Limited predictive power of the model for Java sparrow songs.

Possible reasons include high song learning fidelity and short song length.

Model provides insights into song transmission dynamics.

Abstract

In many species of songbirds, juvenile males learn their songs from adult male tutors. In this paper we formulate a novel Markov model for birdsong transmission developed by analogy with models used in biological sequence analysis. We fit the model using the recently developed Interacting Particle Langevin Algorithm (IPLA) of Akyildiz et al. (arXiv:2303.13429) and analyse a collection of songs from Java sparrows (Lonchura oryzivora) originally recorded and studied by Masayo Soma and her collaborators. The model proves to have limited predictive power for a number of natural problems associated with song transmission in Java sparrows and we propose reasons for this, including the well-established faithfulness of song-learning and the comparative brevity of Java sparrow songs.

Tables

Table 1. Table 1: Formulae for stick-breaking coordinates ϕ {\boldsymbol{\phi}} and their logit-transformed partners 𝒖 {\boldsymbol{u}} in terms of the simplex coordinates 𝒑 {\boldsymbol{p}} . These changes of coordinate allow us to set up the SDEs for the IPLA using unconstrained variables.

Coordinate		Domain	Definition
Simplex	$𝒑=(p_1,\dots ,p_d)$	$𝒑\in {\Delta }_d$
Stick-breaking	$\varphi =({\varphi }_1,\dots ,{\varphi }_{d-1})$	$0\le {\varphi }_j\le 1$	${\varphi }_j=p_j/\left({\sum }_{i=j}^d{p}_i\right)$
Logit-transformed	$𝒖=(u_1,\dots ,u_{d-1})$	$u_j\in ℝ$	$u_j=\log \left({\varphi }_j/\left(1-{\varphi }_j\right)\right)$

Table 2. Table S1: A table summarising the frequency with which the various notes defined in [ 17 ] are sung. The rows of the table correspond to notes and are ordered according to decreasing frequency of occurrence in the data set. The first column gives the note’s name as specified by Lewis et al., the second gives the number of occurrences of the note across all songs used in our study, the third gives the character used to represent the note in alignments such as that illustrated in Figure 4 of the paper and the fourth column gives the note’s number as used in the transmission matrix T T : note that the 8 least commonly-sung notes, which account for around 2.25% of all notes recorded, are amalgamated into a single class. The rightmost two columns give, respectively, the proportion of all notes that are of the given type and the cumulative proportion for the given note and those sung more commonly.

	Number of	Note	Note	Proportion	Cumulative
Name	Occurrences	Char.	Class	of all notes	Proportion
Curve	8470	A	1	0.3687	0.3687
Slope	5855	B	2	0.2549	0.6236
Chip	2946	C	3	0.1282	0.7518
Chirp	1954	D	4	0.0851	0.8369
Stack	1741	E	5	0.0758	0.9127
Curve/Stack	991	F	6	0.0431	0.9558
Whine	310	G	7	0.0135	0.9693
Curve/Slope	188	H	8	0.0082	0.9775
Stack/Curve	143	I	9	0.0062	0.9837
U_Stack	132	J	9	0.0057	0.9895
NLP_Slope	101	K	9	0.0044	0.9939
Noisy_Chirp	97	L	9	0.0042	0.9981
Wavy	23	M	9	0.0010	0.9991
Wavy_Short	10	N	9	0.0004	0.9995
Stack/Slope	9	O	9	0.0004	0.9999
NLP_Other	2	P	9	0.0001	1.0000

Equations

$$\Delta_{d}=\left\{(p_{1},\dots,p_{d})\in\mathbb{R}^{d}\;\biggl{|}\;p_{i}\geq 0\;\;\mbox{and}\;\;\sum_{i=1}^{d}p_{i}=1\right\}.$$

$$P(\mathbf{y}|\mathbf{p},N)\equiv \mathrm{Mult}(\mathbf{y}|\mathbf{p})=\frac{N!}{{\prod }_{i=1}^d{y}_i!}\left(\prod _{i=1}^d{p}_i^{y_i}\right).$$

$$P(\mathbf{p}|\mathbf{\alpha })\equiv \mathrm{Dir}(\mathbf{p}|\mathbf{\alpha })=\frac{\Gamma (\alpha )}{{\prod }_{i=1}^d\Gamma ({\alpha }_i)}\left(\prod _{i=1}^d{p}_i^{{\alpha }_i-1}\right),$$

$$L=\prod _j\mathrm{Mult}({\mathbf{x}}_j|{\mathbf{p}}_j)\times \mathrm{Mult}({\mathbf{y}}_j|{\mathbf{q}}_j)$$

$${\mathbf{q}}_j=T{\mathbf{p}}_j.$$

$$T_{rs}=P(\mbox{pupil sings note $r$}\,|\,\mbox{tutor sings note $s$}).$$

$$\sum_{r=1}^{d}T_{rs}\;=\;\sum_{r=1}^{d}P(\mbox{pupil sings $r$}\;|\;\mbox{tutor sings $s$})\;=\;1.$$

$$\mathcal{L}=\log (L)=\sum _j\log \left(\mathrm{Mult}({\mathbf{x}}_j|{\mathbf{p}}_j)\right)+\log \left(\mathrm{Mult}({\mathbf{y}}_j|T{\mathbf{p}}_j)\right),$$

$$P(D)=\int P(D|\mathbf{\theta })P(\mathbf{\theta })d\mathbf{\theta }.$$

$$P(T|D)=\int P(T,\mathbf{p}|D)d\mathbf{p}=\frac{1}{P(D)}\int P(D|T,\mathbf{p})P(T,\mathbf{p})d\mathbf{p}.$$

$$\int P(D|T,\mathbf{p})P(T,\mathbf{p})d\mathbf{p}.$$

$$U(\theta ,x)\equiv -\log (p_{\theta }(x,y)),$$

$$k(\theta )\equiv p_{\theta }(y)=\int p_{\theta }(x,y)dx=\int e^{-U(\theta ,x)}dx$$

$$d{\theta }_t^N$$

$$d{X}_t^{i,N}$$

$${\theta }_{n+1}$$

$$X_{n+1}^{i,N}$$

$$\mathbb{E}{\left[\Vert T_n-T^{\star }{\Vert }^2\right]}^{1/2}$$

$${\nabla }_{T}k(T)={\nabla }_T\int e^{-U(T,\mathbf{p})}d\mathbf{p}=-\int {\nabla }_{T}U(T,\mathbf{p})e^{-U(T,\mathbf{p})}d\mathbf{p}.$$

$$U(T,\mathbf{p})$$

$$=-\sum _j\log \left(\mathrm{Mult}({\mathbf{x}}_j|{\mathbf{p}}_j)\right)+\log \left(\mathrm{Mult}({\mathbf{y}}_j|T{\mathbf{p}}_j)\right)$$

$$-\sum _{i=1}^d\log \left(\mathrm{Dir}({\mathbf{T}}_{\star i}|{\mathbf{\alpha }}_T)\right)-\sum _j\log \left(\mathrm{Dir}({\mathbf{p}}_j|{\mathbf{\alpha }}_p)\right).$$

$${\boldsymbol{\alpha}}_{p}=\underbrace{(0.5,\dots,0.5)}_{\mbox{$d$ times}}\qquad\mbox{and}\qquad{\boldsymbol{\alpha}}_{T}=\underbrace{(1.1,\dots,1.1)}_{\mbox{$d$ times}}.$$

$${\varphi }_1=\frac{p_1}{p_1+\cdots +p_d},{\varphi }_2=\frac{p_2}{p_2+\cdots +p_d},\dots ,{\varphi }_j=\frac{p_j}{{\sum }_{i=j}^d{p}_i},$$

$${\varphi }_1=\frac{p_1}{1},{\varphi }_2=\frac{p_2}{1-p_1},\dots ,{\varphi }_j=\frac{p_j}{1-{\sum }_{i=1}^{j-1}{p}_i}.$$

$$p_1$$

$$p_d$$

$$u_j=\mathrm{logit}({\varphi }_j)=\log \left(\frac{{\varphi }_j}{1-{\varphi }_j}\right),$$

$${\varphi }_j=\frac{e^{u_j}}{1+e^{u_j}}.$$

$$U(T,\mathbf{p})$$

$$-\sum _{r=1}^d\sum _{i=1}^d({\alpha }_T-1)\log \left(T_{ri}\right)$$

$$=-\sum _j\sum _{i=1}^d({\alpha }_p-1+x_{ij})\log (p_{ij})+y_{ij}\log \left(\sum _{r=1}^d{T}_{ir}{p}_{rj}\right)$$

$$-\sum _{r=1}^d\sum _{i=1}^d({\alpha }_T-1)\log \left(T_{ri}\right)$$

$$\left(\prod _j\prod _{i=1}^d{p}_{ij}\right)\times \left(\prod _{r=1}^d\prod _{i=1}^d{T}_{ri}\right)$$

Full text

Keywords birdsong  $\cdot$ social learning  $\cdot$ maximum marginal likelihood estimation  $\cdot$ bridge sampling  $\cdot$ stochastic gradient ascent

1 Introduction

In many species of songbirds, juvenile males learn their songs from adult male tutors. In this paper we ask whether tutor-to-pupil song transmission can be studied by analogy with molecular sequence evolution. Where songs consist of a sequence of notes, genomes are sequences of nucleotides. Like nucleotides, notes can change from one type to another due to learning errors (i.e. point mutations). Individual notes or even song segments can be inserted, deleted, or duplicated, just like nucleotides and motifs in genomes. Hence, tools from molecular phylogenetics may provide a foundation for modelling birdsong transmission. However, there are important differences between genomes and birdsong. For example, individuals have a single genome, but a bird may sing many different songs. Although there are only four nucleotide bases, birds can sing many different note types. Moreover, molecular phylogenetics is an advanced field where estimates of mutation rates and substitution matrices are available, while birdsong transmission is an emerging subject where such tools are not yet available. These differences necessitate the development of novel statistical approaches.

1.1 Java sparrow songs

Java sparrows (Lonchura oryzivora) are a highly vocal estrildid finch native to the island of Java, Indonesia ([13]). Male birds typically sing between 2 and 8 note types and their songs often include repetitive syllables, trills, and non-trill notes. Previous studies have established that male juveniles (pupils) learn their songs from their so-called social fathers, whom we will call tutors in this paper. ([26]). Learning fidelity is high: pupils learn their tutors’ note repertoires, note sequences (song structure) and temporal features ([17, 16]). The strong effect of social learning provides the rationale for creating a data-driven, model for the transmission of note sequences.

1.2 Overview of the data

We study a dataset published in ([17]) consisting of whole song recordings from a captive population of Java sparrows, reared by the Soma lab at the University of Hokkaido. The birds have been selectively tutored, so that almost all pupil-tutor relationships are known. Recordings were made in a soundproof chamber with each bird singing in isolation and all recordings for a given bird were made within a single week (time point). For each bird, when more than 10 full song recordings were available from a single time point, 10 were randomly selected for inclusion. If fewer than 10 recordings were available at all time points, all songs from one time point were used, prioritising recordings made when the bird was approximately 2–5 years old. Full details of the data collection procedures are provided in Lewis, Soma et al. (2021). In total, there are 676 songs recorded from 73 individuals (mean 9.3 songs per individual, range 3 to 10), including 58 pupil-tutor pairs. The birds came from 10 different social lineages, named with arbitrary colours by the Soma lab. A social lineage is a chain of birds all of whose members can trace their songs back to the same founder: see Figure 1. Lineages varied from 1 to 4 generations in length.

Figure˜2 is a visual representation of one song in the Java sparrow dataset. The coloured waves are the individual sound elements (i.e. notes) which are separated by intervals of near-silence. [17] inspected such spectrograms and excised notes using the software tool Koe ([8]). They then manually classified the notes based on their shape and spectral features: a full discussion of the note classification procedure is provided in [17]. In total, Lewis et al. classified 22,972 notes into 16 note classes: see Figure˜3 for typical examples of the note types. There is now software to computationally excise and classify notes with high reliability using tools such as TweetyNet ([4]). By excising and classifying the notes, Lewis et al. converted each song into a note sequence. We generated multiple sequence alignments of note sequences from tutors and pupils using profile hidden Markov models ([15]: in preparation).

1.3 Reduction of songs to count data

Figure˜4 illustrates the kind of data that we analyze. It consists of a set of aligned note sequences from a tutor and a pupil. We restricted attention to those columns of the alignment where both birds have notes in at least half of their songs (as opposed to gap characters). Next, we reduced the alignment to a pair of matrices of counts:

[TABLE]

For example, we write ${\boldsymbol{x}}_{j}=(x_{1j},\dots,x_{dj})$ to indicate the vector of the tutor’s counts at the $j$-th aligned position and ${\boldsymbol{y}}_{j}$ for the pupil’s counts. We will do the same thing for all pupil-tutor pairs and so the index $j$ will run across all aligned positions, both within the alignment for a single known pupil/tutor pair and across all such pairs. In Section 2 below, we develop a probabilistic model for the generation of such counts.

1.4 Outline of the paper

The remainder of the paper is organised as follows. Section 2 lays out the mathematical foundations of our model and the associated inference problems. In particular, in Section 2.1 we develop a probabilistic model for birdsong transmission using an analogy to models used in molecular phylogenetics. A key ingredient is the tranmission matrix $T$, whose element $T_{rs}$ is the probability that a pupil will sing note $r$ given that his tutor sang note $s$. In Section 2.2 we formulate a Bayesian inference problem for the model’s parameters and in Sections 2.3—2.4 we explain how to use the recently developed Interacting Particle Langevin Algorithm (IPLA) of [1]—which is essentially an optimisation algorithm—to estimate $T$.

One of the attractions of the IPLA is that, subject to certain restrictions on the optimisation target, it comes with strong convergence guarantees. In Section 3 we ask whether our optimisation target, a certain marginal likelihood, satisfies these conditions and conclude that it does not: full proofs are supplied in an Appendix. To address this problem, we did numerical experiments in which we applied the IPLA to synthetic data generated from our model. We report on these in Section 4.1, where we find that despite the lack of theoretical guarantees, the inference problem shows no evidence of multimodality and we are able to recover the transmission matrix. Bolstered by these observations, we fit a transmission matrix for the birdsong data: these results are reported in Sections 4.2 and 4.3. All these calculations rely on codes that we developed ourselves and ran in R [20], accelerated with the Rcpp package [7, 6].

Armed with the fitted transmission matrix, we then used the Hamiltonian Monte Carlo package Stan [3] and the bridgesampling R package [11] to estimate the posterior probabilities that each possible pair of birds had a pupil-tutor relationship: these results are summarised in Section 4.4. The paper concludes with a discussion in Section 5 and an Appendix providing detailed proofs.

2 Notions and notation

Here, we define terms and introduce tools that we will use to specify our probabilistic model.

The $d$-dimensional simplex $\Delta_{d}$ is the set

[TABLE]

One can think of points in the interior of $\Delta_{d}$ as discrete probability distributions over the numbers $1,\dots,d$.

The multinomial distribution is a probability distribution over vectors of counts ${\boldsymbol{y}}\in\mathbb{N}^{d}$ that is parameterised by the total number of counts $N=\sum_{i=1}^{d}y_{i}$ and a vector of probabilities ${\boldsymbol{p}}\in\Delta_{d}$. Its probability mass function is

[TABLE]

The Dirichlet distribution is a probability distribution over $\Delta_{d}$ and is parameterised by a vector of positive shape parameters ${\boldsymbol{\alpha}}\in\mathbb{R}_{+}^{d}$. Its probability density function is

[TABLE]

where $\alpha=\sum_{i=1}^{d}\alpha_{i}$ is the sum of the shape parameters.

2.1 A probabilistic generative model

If we then introduce the following unobserved (and, indeed, unobservable) quantities

[TABLE]

then the likelihood for the count data is:

[TABLE]

where the product ranges over all aligned positions for all pupil-tutor pairs. Here ${\boldsymbol{p}}_{j}=(p_{1j},\dots,p_{dj})$ is the tutor’s vector of note-usage probabilities at the $j$-th aligned position and ${\boldsymbol{q}}_{j}$ is the pupil’s.

In the spirit of substitution models used in phylogenetics, we assume that the pupil’s probabilities ${\boldsymbol{q}}_{j}$ are related to the tutor’s ${\boldsymbol{p}}_{j}$ by a single, fixed transmission matrix $T$ that is shared across all aligned positions, so that

[TABLE]

Here $T$ is a $d\times d$ matrix and

[TABLE]

As the pupil has to sing something, the columns of $T$ must satisfy:

[TABLE]

Hence, we obtain a new expression for the log-likelihood:

[TABLE]

where, as above, the sum ranges over all aligned positions for all pupil-tutor pairs.

2.2 Bayesian inference

Given we have the aligned note sequences and a transmission model, we need only specify a set of priors to set up a Bayesian inference problem. Let $D$ denote the data (the counts ${\boldsymbol{x}}_{j}$ and ${\boldsymbol{y}}_{j}$ from the alignment) and ${\boldsymbol{\theta}}$ be the parameters ($T$ and the latent probabilities ${\boldsymbol{p}}_{j}$). Then we write

$P(D\,|\,{\boldsymbol{\theta}})$

the likelihood. $P(D\,|\,{\boldsymbol{\theta}})=\prod_{j}\mathop{\textrm{Mult}}({\boldsymbol{x}}_{j}\,|\,{\boldsymbol{p}}_{j})\times\mathop{\textrm{Mult}}({\boldsymbol{y}}_{j}\,|\,T{\boldsymbol{p}}_{j})$.

$P({\boldsymbol{\theta}})$

the prior on the parameters. We use $d$-dimensional Dirichlet distributions here: one for every column of $T$ and one for each ${\boldsymbol{p}}_{j}$.

$P({\boldsymbol{\theta}}\,|\,D)$

the posterior distribution over ${\boldsymbol{\theta}}$, $P({\boldsymbol{\theta}}\,|\,D)\propto P(D\,|\,{\boldsymbol{\theta}})P({\boldsymbol{\theta}}).$

$P(D)$

the marginal likelihood. This is a value from the marginal distribution obtained by integrating ${\boldsymbol{\theta}}$ out of $P(D,{\boldsymbol{\theta}})$

[TABLE]

2.3 Maximal marginal likelihood estimation of $T$

It would be straightforward to use Hamiltonian Monte Carlo (HMC) to draw samples from the posterior over $T$ and the ${\boldsymbol{p}}_{j}$. But as we are only interested in $T$, the ${\boldsymbol{p}}_{j}$ are nuisance parameters. Ideally, we would compute or sample from the marginal distribution obtained by integrating them out of the posterior:

[TABLE]

where ${\boldsymbol{p}}$ is a vector of all the latent probabilities of note-usage ${\boldsymbol{p}}_{j}$ and $T$ is the transmission matrix.

If we could compute $P(T\;|\;D)$, we could then maximise it with respect to the elements of $T$ to obtain a maximal marginal likelihood estimate (MMLE) for $T$. And since the marginal likelihood $P(D)$ does not depend on $T$, it would be sufficient to maximise

[TABLE]

Unfortunately the integral in (9) is intractable and so we are forced to resort to numerics, using a recently-developed algorithm described below.

2.3.1 MMLE with stochastic differential equations

Recent work of [1], inspired by the methods of [14], suggests a diffusion-based approach to maximal marginal likelihood estimation, the interacting particle Langevin algorithm (IPLA), that is a member of a family of algorithms that generalise the EM-algorithm using ideas that originated in [18]. The setup includes observed data denoted by $y$ (for us this would be the count data), parameters that one wishes to estimate denoted by $\theta$ and latent variables, denoted by $x$, that one wishes to marginalise-out. They begin by defining the negative log-likelihood,

[TABLE]

where $p_{\boldsymbol{\theta}}(x,.)$ is the joint probability distribution of the latent variables $x$ and observed data $y$ given the parameters $\theta$. The goal of IPLA is to maximise the marginal likelihood:

[TABLE]

[1] propose to maximise $k(\theta)$ via gradient ascent, approximating $\nabla k(\theta)$ by averaging over $N$ of independent realisations of a stochastic differential equation (SDE) on the full parameter space:

[TABLE]

for $i=1,\dots,N.$ Here $X^{i,N}_{t}$ are the $N$ particles used for estimating the gradient and $\{(B_{t}^{i,N})_{t\geq 0}\}^{N}_{i=0}$ is a family independent Brownian motions. In practice, Akyildiz et al. use an Euler-Maruyama discretisation to generate approximate realisations. Given $(\theta_{0},X_{0}^{1:N,N})\in\mathbb{R}^{d_{\theta}}\times(\mathbb{R}^{d_{x}})^{\otimes N}$ and for any $n\in\mathbb{N}$

[TABLE]

where $\gamma>0$ is a step size and $\xi_{n}^{i,N}=B_{n\gamma}^{i,N}-B_{(n-1)\gamma}^{i,N}$ for $0\leq i\leq N$ are i.i.d. samples from the standard normal.

An advantage of the IPLA is that, subject to certain smoothness and convexity assumptions on $U$, one can obtain provable bounds on:

• •

the range of step sizes $\gamma$ that will yield convergence to $T^{\star}$, the MMLE for $T$;

• •

the expected deviation of $T_{n}$ from the MMLE, $T^{\star}$:

[TABLE]

The second set of bounds depends on the step-number $n$ and so also provides a bound on the rate of convergence of the algorithm. These smoothness and convexity assumptions, which are discussed—and, where possible, verified—for our application in Section 3 below, also ensure that

[TABLE]

2.3.2 The IPLA for birdsong transmission

For our application the role of Akyldiz’s data $y$ is played by the count data, the latent variables $x$ are our ${\boldsymbol{p}}_{j}$, the collection of latent note probabilities, and $\theta$ is the transmission matrix $T$. The optimisation target, a negative log-posterior, is

[TABLE]

Here ${\boldsymbol{T}}_{\star i}$ is the $i$-th column of $T$ and the shape parameters for the priors are

[TABLE]

This choice of ${\boldsymbol{\alpha}}_{p}$ corresponds to a distribution that favours ${\boldsymbol{p}}_{j}$ in which most of the entries are small, while the Dirichlet distribution with ${\boldsymbol{\alpha}}_{T}$ favours weakly $T_{*j}$ with uniform probabilities.

Unfortunately one cannot use the SDEs in Eqn. (12) directly for either the columns ${\boldsymbol{T}}_{\star i}$ or the vectors of probabilities of note selection ${\boldsymbol{p}}_{j}$. The issue is that these quantities must lie in the $d$-simplex, but the SDEs cannot enforce such constraints. Thus to apply the IPLA we must transform the ${\boldsymbol{p}}_{j}$ and columns of $T$ onto unconstrained scales: we use logit-transformed, stick-breaking coordinates, which are explained below.

2.4 Unconstrained coordinates

The problem of needing parameters to lie in a simplex when using an algorithm that requires unconstrained variables is well known and we adopted the approach used by the developers of the HMC package Stan [3]. Alternative coordinate systems exist, (e.g. softmax coordinates or the system described in [2]), but as we use Stan in Section 4.4, we chose to minimise the number of coordinate systems we needed to introduce.

The first thing to note is that the $d$-simplex is a $(d-1)$-dimensional object as, given the first $(d-1)$ of the $p_{j}$, we can solve the constraint $\sum_{j=1}^{d}p_{j}=1$ for the remaining component $p_{d}$. A point ${\boldsymbol{p}}=(p_{1},\dots,p_{d})$ in the $d$-simplex $\Delta_{d}$ will thus be associated with $d-1$ unconstrained coordinates. We begin by constructing stick-breaking coordinates ${\boldsymbol{\phi}}=(\phi_{1},\dots,\phi_{d-1})\in[0,1]^{d-1}$. To define these coordinates—and to understand the origin of their name—imagine starting with a stick of unit length and snapping-off a piece of length $p_{1}$, then one of length $p_{2}\;\dots$ and so on. The stick-breaking coordinate $\phi_{j}$ is then defined as the fraction of the remaining stick that gets snapped-off at step $j$ in this process. Figure 5 illustrates the construction for an example with ${\boldsymbol{p}}\in\Delta_{3}$, but in general we have

[TABLE]

where the denominator of $\phi_{j}$ is $p_{j}+\dots+p_{d}$, the length of the stick that remains after the first $j-1$ pieces have been broken off. In what follows it will prove convenient to rewrite the denominators in terms of the pieces already removed, so that the formulae become:

[TABLE]

Next we invert the transformation to give

[TABLE]

Although the sum of the stick-breaking coordinates is no longer constrained, the $\phi_{j}$ are still restricted to satisfy $0\leq\phi_{j}\leq 1$ and so are unsuited to use in the IPLA. We thus apply the logit transformation to the components of ${\boldsymbol{\phi}}$ to obtain ${\boldsymbol{u}}=(u_{1},\dots,u_{d-1})\in\mathbb{R}^{d-1}$ defined by

[TABLE]

a transformation that’s also simple to invert:

[TABLE]

Table 1 summarises the coordinate changes required to convert between ${\boldsymbol{p}}\in\Delta_{d}$ and the unconstrained variables needed for the IPLA. In the remainder of the article, we will use the following notation:

• •

$\Phi:\Delta_{d}\rightarrow[0,1]^{d-1}$ will be the map that sends a point ${\boldsymbol{p}}\in\Delta_{d}$ to its associated stick-breaking coordinates, as given by Eqn. (17).

• •

$\Lambda:[0,1]^{d-1}\rightarrow\mathbb{R}^{d-1}$ will be the map that sends stick-breaking coordinates to their logit-transformed counterparts, as given by Eqn. (19).

• •

${\boldsymbol{u}}_{j}=\Lambda(\Phi({\boldsymbol{p}}_{j}))$ will be the $(d-1)$-dimensional vector of logit-transformed, stick-breaking coordinates associated with the tutor’s note-usage probabilities ${\boldsymbol{p}}_{j}$;

• •

${\boldsymbol{\tau}}_{k}=\Lambda(\Phi({\boldsymbol{T}}_{\star k}))$ will be the $(d-1)$-dimensional set of logit-transformed, stick-breaking coordinates associated with the $k$-th column of the transmission matrix.

3 Establishing convergence

The proofs of convergence in [1] depend on three assumptions, the first two of which involve the function $U({\boldsymbol{u}},{\boldsymbol{\tau}})$ defined by transforming the log density in Eqn. (14) to a density over the unconstrained variables. In this section we introduce these assumptions and check whether our $U({\boldsymbol{u}},{\boldsymbol{\tau}})$ satisfies them. In brief, it doesn’t. Although we can arrange our implementation of the IPLA so that two of the three conditions are satisfied, we can also prove that the remaining one is not. In the rest of this section we give precise statements of the assumptions and remark briefly on the extent to which our system satisfies them.

Nevertheless, numerical results presented in Section 4 suggest that the IPLA does converge to a sensible MMLE for $T$ when applied to synthetic data and produces plausible results when applied to Java sparrow song. This is in keeping with the discussion in [1], where the authors say

$\dots$ our results are obtained in the case where $U$ is gradient Lipschitz and obeys a strong convexity condition, but we believe that similar results can be obtained under much weaker (nonconvex) conditions, [27].

3.1 Smoothness and regularity assumptions

After suppressing terms that do not depend on the ${\boldsymbol{p}}_{j}$ or $T$ and combining terms in $\log(p_{ij})$, the negative log-posterior function is:

[TABLE]

where the ${\boldsymbol{p}}_{j}$ and the columns of $T$ are constrained to lie in the simplex $\Delta_{d}$.

To transform this density to one over the unconstrained coordinates ${\boldsymbol{u}}_{j}$ and ${\boldsymbol{\tau}}_{k}$, we need to include a factor arising from the Jacobian of the change of variables. Calculations detailed in Section A.1 of the Appendix establish that the Jacobian factor is just the product

[TABLE]

and so the optimisation target for the IPLA becomes

[TABLE]

where the $p_{ij}$ and the $T_{ri}$ are to be regarded as functions of the unconstrained coordinates.

To state the conditions, we need the following definitions:

Definition** (Lipschitz function).**

A function $f:X\rightarrow Y$ between two Banach spaces is said to be Lipschitz continuous if, for all ${\boldsymbol{x}}_{1},{\boldsymbol{x}}_{2}\in X$ we have, for some fixed constant $K$, that

[TABLE]

Any such $K$ is called a Lipschitz constant for $f$ and when such a constant exists, we say that “$f$ is Lipschitz".

Definition** (Strong convexity).**

A differentiable function $f:\mathbb{R}^{n}\rightarrow\mathbb{R}$ is said to be strongly convex if, for all ${\boldsymbol{v}},{\boldsymbol{v}}^{\prime}\in\mathbb{R}^{n}$ we have, for some fixed constant $\alpha$, that

[TABLE]

where the angle braces on the left hand side denote a Euclidean inner product. Strong convexity means, roughly, that the growth of $f$ is at least quadratic, but there are a host of equivalent conditions: see, for example, [28].

The proofs of convergence in [1] as applied to our problem then depend on the following assumptions:

A1:

If the count data are held fixed, then $\nabla_{{\boldsymbol{\theta}}}U$ is Lipschitz in ${\boldsymbol{\theta}}=({\boldsymbol{u}},{\boldsymbol{\tau}})$ where ${\boldsymbol{u}}$ encompasses of all the ${\boldsymbol{u}}_{j}$. We prove this in Appendix A.

A2:

The function $U({\boldsymbol{u}},{\boldsymbol{\tau}})$ is strongly convex. This assumption is violated for our problem and we provide a proof in Appendix A.

A3:

If we fix $N$, the number of particles (that is, random walkers) in Eqns. (12) and (13) of the IPLA, then the distribution from which we draw the initial condition should be such that

[TABLE]

has bounded second moment. We have arranged for this to be true, as is also shown in Appendix A.

4 Numerical results

Here we report briefly on two numerical experiments, one with synthetic data generated by the model in Section 2.1 and another with the birdsong data from [17].

4.1 Synthetic data

We generated synthetic count data as follows:

• •

Fixed a number of notes $d=5$, a number of aligned positions $n=500$ and a number of recordings per bird $m=8$. The values of $d$ and $m$ are similar to those in our real data, but we chose $n$ rather smaller, so the computations would finish more quickly.

• •

Generated a $d\times d$ transmission matrix $T$ whose columns were drawn from a Dirichlet distribution with shape parameters $\alpha_{i}=\alpha_{T}=1.0$ for $i\in\{1,\dots,5\}$. This corresponds to a uniform distribution over $\Delta_{5}$.

• •

Generated note-usage probabilities for the tutors ${\boldsymbol{p}}_{j}\in\Delta_{5}$ for $j\in\{1,\dots,500\}$ by drawing them from a Dirichlet distribution with shape parameters $\alpha_{i}=\alpha_{p}=0.5$. This favours ${\boldsymbol{p}}_{j}$ similar to observed data in that some of the entries are much, much smaller than others (see Section 4.2 for an overview of the real data).

• •

Generate note-usage counts ${\boldsymbol{x}}_{j}$ for the tutors by drawing $(x_{1j},\dots,x_{5j})$ from the multinomial distribution with $m=8$ total counts and probabilities given by ${\boldsymbol{p}}_{j}$.

• •

Use $T$ and the ${\boldsymbol{p}}_{j}$ to compute note-usage probabilities for the pupils using the model ${\boldsymbol{q}}_{j}=T{\boldsymbol{p}}_{j}$.

• •

Generate note-usage counts $y_{ij}$ for the pupils by drawing from the multinomial distribution as sketched above.

The virtue of this approach is that we know the true value of $T$. Figure 6 summarises a numerical experiment to estimate $T^{\star}$ using the IPLA applied to count data generated as above and suggests that the estimates, though sometimes biased, are reasonably accurate, despite the lack of strong convexity in the negative log posterior $U({\boldsymbol{u}},{\boldsymbol{\tau}})$.

We also investigated this example for the possibility of multimodality in the marginal likelihood by repeating the IPLA estimation of $T^{\star}$ with $N=512$ particles 200 times, initialising the algorithm differently each time. We thus obtained, for each of the $d\times d$ entries $T^{\star}_{ij}$ in the transmission matrix, a collection of 200 estimates and we examined the distribution of these estimates, both visually (histograms are available in Figure S2 of the Supplementary Material) and with the dip test of unimodality developed in [12]. This is a hypothesis test whose null hypothesis is that the data are drawn from a distribution that has a continuous, unimodal density. In all cases, the estimates for the $T^{\star}_{ij}$ appeared sharply concentrated near a single value and the dip test suggested little evidence of deviation from unimodality ($p>0.7$ in all cases).

4.2 Data from recorded birdsong

Building on work reported by [17], we aligned note sequences derived from separately-recorded songs for pupil-tutor pairs and focused attention on positions where fewer than half of each bird’s aligned songs contained gap characters. Although [17] defined 16 note types, many of them are sung very infrequently (see Figure˜7 and Table S1 in the Supplementary Material), thus we amalgamated the eight rarest types into a single “Other" category to produce a total of nine note types. Next we used the IPLA to estimate a note transmission matrix $T^{\star}$, producing results illustrated in Figure˜8 and Figure˜9. For the most commonly-sung notes, $T^{\star}$ is very close to the identity, reflecting the faithfulness of song transmission: pupils typically learn to reproduce their tutors’ notes accurately. For the three rarest note types, including the amalgameted note class 9, pupils appear sometimes to replace the tutor’s note with a more commonly sung one.

4.3 Transmission Matrix

We fitted a transmission matrix for the eight most commonly used notes, along with a ninth, amalgamated note-type for the most rarely-sung notes, using IPLA (Figure˜9). $T^{\star}_{ij}$ is the estimated probability that a pupil sings note $i$, given the tutor sings note $j$ in a given position. The diagonal elements corresponding to the six most commonly-sung notes are very close to one, suggesting that the learning fidelity is high overall. The exceptions are notes 7, 8 and 9, which have probabilities of 0.68, 0.8 and 0.3, respectively, of changing to one of the most commonly-sung notes, 1 or 2. Thus, we expect only small changes in the note sequences per generation, unless the tutor often sings notes of type 7, 8 or one of the rare notes aggregated into class 9.

4.4 Evidence

We used bridge sampling ([9, 10]) as implemented in the bridgesampling R package of [11], in combination with the HMC sampler Stan ([3]) to estimate the probability that any one bird tutored another. That is, for each potential (pupil, tutor) pair $(a,b)$, we aligned the birds’ songs and then estimated

[TABLE]

where here the index $j$ ranges over all aligned positions for the putative pupil-tutor pair, ${\boldsymbol{x}}_{j}$ is the putative tutor’s note usage counts at position $j$, ${\boldsymbol{y}}_{j}$ are those of the putative pupil and ${\boldsymbol{\alpha}}_{p}$ is as in Eqn. (15). We will refer to this quantity as the evidence for the pair and, as it depends on the number $n_{ab}$ of aligned positions, will work with a related quantity

[TABLE]

the mean log evidence per aligned site, to permit comparisons between different putative pupil-tutor pairs.

In addition to bridge sampling, we also estimated the evidence with a harmonic mean estimator (see [22] and [10]). The results typically agreed well with those from bridge sampling (differences of less that 2% in the estimates of $e_{ab}$ for 94% of putative pairs), but were more computationally expensive to obtain and had a larger variance than those from bridge sampling, so we haven’t reported them.

Figure˜10 shows log evidence per aligned site as specified by Eqns. (26) and (27). The blocks of high log evidence near the diagonal align broadly with the social lineages, suggesting that our method can pick up some social relationships among birds. However, the considerable variation within these blocks suggest that this is not always the case. Furthermore, the evident symmetry of the matrix suggests that our method may struggle to identify the direction of transmission. We were unable to compute the evidence for some pairs, either because we were unable to align their songs or because there were insufficiently many aligned positions. This occurs when the song sequences are too dissimilar (e.g. birds sing completely different notes most of the time).

First, we considered whether our method can accurately discern the direction of song tutoring within the 57 pupil-tutor pairings provided in the original data. We compared the log evidence for the true pairings with those in which the tutor and pupil swapped roles and performed a paired Wilcoxon signed-ranked test (0.05 significance) to determine whether the true pairings had a significantly higher log evidence than the reverse. Our results were insignificant ($p$-value 0.3217, V = 581), thus we could not correctly determine the direction of song tutoring. This is consistent with our informal observation of the symmetry of Figure˜10.

Next, we considered whether we could use the log evidence results to deduce pupil-tutor pairings, regardless of the direction. For every pupil in the original dataset, we took our predicted tutor to be the bird with the highest log evidence. We compared this against 57 known pupil-tutor pairs and found the method to be correct only 2 times. This is also consistent with Figure˜10 where the columns, representing the pupils, do not typically have a single pixel showing substantially higher log evidence than the rest. Instead, there are several pixels showing high evidence, indicating several potential song tutors.

Moving beyond individual pairings, we investigated whether our method could distinguish between first and second generation pupil-tutor relationships. For all cases where the comparison was possible, we compared the log evidence for the pupils and their tutors with that for the pupils and their grand tutors (tutor of their tutor). As above, we compared the two sets of values using a paired Wilcoxon signed-ranked test (0.05 significance). The results were insignificant (V = 166, $p$-value = 0.2221).

Finally, we investigated whether social lineages could be detected by our method. First, we classified each putative pupil-tutor pair according to whether both birds were part of the same lineage or not. We then compared the two sets of log evidence values—those from within-lineage pairs and those from cross-lineage pairs—graphically: see Figure˜11. Although the log evidence tended to be higher for the within-lineage pairs, there was some overlap in the inter-quartile ranges and outliers in both sets. We also tried to use log evidence to predict a bird’s lineage by considering all putative (tutor, pupil) pairs in which the given bird played a role, associating the pair with the lineage of the given bird’s partner. We then assigned the bird to the lineage whose pairs had the highest mean log evidence.

We noticed that the two Dark Pink birds (JS299, JS155) had very high log evidence values with most other birds, excluding the cases where the evidence could not be estimated. These two birds to sing mostly A’s, interspersed with some B’s. Since A is the most commonly sung note, this could explain why the birds gave unusually high log evidence values. When we remove these two birds from the analysis, then out of the 70 remaining birds, the approach sketched above assigned 39 to the correct lineage. This corresponds to a predictive accuracy of $\approx 55\%$, whereas a naive classifier using the most common class would have $\approx 21\%.$ We also evaluated our predictions using the Adjusted Rand Index (ARI), which measures the similarity between two clusterings while accounting for chance agreement. Unlike classification accuracy, ARI captures the overall grouping structure rather than point-by-point label matches. Our predictions achieved an ARI of 0.2763, indicating performance better than random labelling (ARI = 0) but well short of a perfect match (ARI = 1), consistent with the observed classification accuracy of 55%. Both the ARI and percentage accuracy are consistent with Figure˜10, where the blocks along the diagonal, which align broadly with the lineages, are not regions of consistently high log evidence. Thus our model made modestly accurate predictions of social lineage. This is consistent with Figure˜10, where the blocks along the diagonal, which align broadly with the lineages, are not regions of consistently high log evidence.

5 Conclusion

Using Java sparrows as a model organism, we have presented the first sequence evolution model of birdsong transmission. By applying the IPLA of [1], we estimated a note transmission matrix for the entire population. We regard this as a first approximation: the variation in note-usage pattern evident in Figure˜7 might suggest the use of a Bayesian hierarchical model in which the prior for the note usage probabilities varies from lineage to lineage. As a further simplification, we lumped the eight rarest note-types into a single catch-all category, numbered 9 in our scheme.

Generally, our estimated transmission matrix is consistent with the well-established faithfulness of song transmission: most diagonal elements $T^{\star}_{jj}$ are very close to 1.0. The three exceptions include the two comparatively rarely-sung note classes seven and eight—Whine and Curve/Slope, respectively—which pupils sometimes convert to the most commonly-sung note type, Curve. The catch-all class of very rare notes, nine, shows similar conversions to the two most commonly-sung note classes. This may be an artifact of the alignment process, during which phrases containing rarely-sung notes may be aligned with similar phrases that lack these notes. Ideally, one would incorporate the probabilistic models underpinning alignment along with those representing transmission into one grand estimation problem, but this would be computationally intractable.

Using the estimated note transmission matrix, we computed bridge sampling estimates for the log evidence for all potential pupil-tutor relationships. Although our model was not able to detect pairwise relationships (e.g. correctly match pupils to tutors) or the direction of song tutoring, we did achieve modest results in assigning pupils to their correct lineages. The key issue appears to be the faithfulness of song transmission and the comparative sparsity of data, particularly the song lengths. In bioinformatic contexts, where mistakes during, for example, DNA replication are also very rare, researchers are typically working with kilobases or even megabases worth of data. For example, the median length of a human gene is 24 kilobases ([24]) and the human genome is $\approx$ 3 billion base pairs long. Even with a very small mutation rate of $0.5\times 10^{-9}$ ([23]), there are still a substantial number of variants to study. By contrast, our songs are all less than 90 notes long and their collections of aligned phrases are even shorter, making observed note changes far rarer.

Furthermore, our sequence evolution model lies downstream of the preprocessing of the song recordings conducted by [16, 17]. There are two main preprocessing steps; namely the segmentation of notes from the song spectrograms and the classification of notes into note classes. The recordings were taken in controlled conditions with minimal background noise, so it was relatively simple to manually identify the start and end points of notes. After segmentation, [17] manually classified the notes using a suite of acoustic characteristics. They repeated the classification with a second observer, achieving a consensus of $97.5\%.$ Although our model is dependent on the work of [17], we are confident that this was done as accurately as anything in the field. Modern tools such as TweetyNet ([4]), which uses a neural network to preprocess song recordings from a range of birds could increase the pool of available data substantially and expand the scope for sequence evolution models to other species of songbirds.

Additionally, our sequence evolution model relies on the multiple sequence alignments of the pupil and tutor songs. We produced these alignments using profile hidden Markov models, which treat transmission at all sites as independent of each other. Arguably this assumption doesn’t hold as birds’ songs are known to contain motifs. For example, if a bird sings “ACACAC”, then the probability of singing a C depends on whether the bird has sung an A immediately before. The same challenge is found in bioinformatics where genetic data may contain a variety of motifs and repeated elements. Despite this problem, researchers have still been able to derive useful biological insights from their alignments ([21, 25]). In our case, we inspected the alignments to check that they identified regions of consensus: see Figure˜4 for an example. Whilst relaxing the independent sites assumption may be more accurate, it would involve formulating a new, complex approach to multiple sequence alignment that would be beyond the scope of this paper.

To establish that our approach would perform better given enough aligned motifs, we could in principle repeat our analysis with synthetic data. Currently the median number of aligned positions per pair is 18, but, using the profile hidden Markov models of [15], we could in principle generate songs with arbitrarily long motifs, but the computations reported here already lie near the limit of our resources.

Acknowledgements

The authors acknowledge the assistance provided by the Research IT team and the use of the Computational Shared Facility at The University of Manchester. We also acknowledge Dr. R. Tucker Gilman who was Kwong’s co-supervisor. Finally, we thank Dr. Rebecca N. Lewis, Professor Masayo Soma (University of Hokkaido) and the members of her lab, for giving us access to their data and sharing their expertise.

Funding

Kwong was supported by EPSRC Doctoral Training Programme grant EP/W524347/1.

Supplementary Material

A schematic diagram illustrating the coordinate systems used in the paper; a set of histograms illustrating the lack of multimodality in the IPLA estimates of the entries $T^{*}_{ij}$ in the transition matrix from Section 4.1 and a table listing the note types appearing in our data along with the frequencies with which they are sung.

Appendix A Smoothness and regularity conditions: proofs

In this section we prove the following:

Proposition 1**.**

When one transforms the negative log posterior density whose ${\boldsymbol{p}}$ and $T$-dependent terms are listed in Eqn. (21) to the unconstrained coordinates, the terms that depend on ${\boldsymbol{u}}$ and ${\boldsymbol{\tau}}$ are

[TABLE]

where $\Lambda^{-1}$ and $\Phi^{-1}$ are the inverses of the functions $\Lambda$ and $\Phi$ defined in Table 1, $\Phi^{-1}_{i}({\boldsymbol{\phi}})$ is the $i$-th component of $\Phi^{-1}({\boldsymbol{\phi}})$, ${\boldsymbol{u}}_{j}\in\mathbb{R}^{d-1}$ is the vector of unconstrained coordinates associated with the tutor’s note-usage probabilities at alignment position $j$ and ${\boldsymbol{\tau}}_{k}$ is the vector of unconstrained coordinates associated with the $k$-th column of the transition matrix.

Theorem 1**.**

The gradient of $U({\boldsymbol{u}},{\boldsymbol{\tau}})$ is Lipschitz and so assumption A1 of Akylidiz et al. is satisfied.

Theorem 2**.**

$U({\boldsymbol{u}},{\boldsymbol{\tau}})$* is not strongly convex and so assumption A2 of Akylidiz et al. is violated.*

Theorem 3**.**

Let $\tau_{ij,0}$ denote the initial condition for the unconstrained coordinate $\tau_{ij}$ in the IPLA and let $u_{ij,0,k}$ indicate the initial condition for the unconstrained coordinate $u_{ij}$ for the $k$-th particle. The distributions from which these initial conditions are drawn are such that there exists a constant $H$ such that, for all $k$,

[TABLE]

and so assumption A3 of Akylidiz et al. is satisfied.

A.1 Preliminaries: changes of coordinate

We will need to compute derivatives with respect to the unconstrained coordinates, so here we compute the Jacobians of the coordinate changes from stick-breaking to simplex coordinates and from logit-transformed variables to those constrained to lie in the unit interval. The first Jacobian we will need is $J_{{\boldsymbol{\phi}}}$, a $d\times(d-1)$ matrix that relates stick breaking coordinates ${\boldsymbol{\phi}}=(\phi_{1},\dots,\phi_{d-1})$ to points ${\boldsymbol{p}}=(p_{1},\dots,p_{d})$ in the $d$-simplex. It is given by

[TABLE]

where

[TABLE]

which follows from the expressions in Eqn. (18). In addition to $J_{\boldsymbol{\phi}}$, we will need $J_{{\boldsymbol{u}}}$, a $(d-1)\times(d-1)$ matrix that relates the stick-breaking coordinates ${\boldsymbol{\phi}}=(\phi_{1},\dots,\phi_{d-1})$ to their logit-transformed partners ${\boldsymbol{u}}=(u_{1},\dots,u_{d-1})$. As $\phi_{j}$ depends on $u_{j}$ alone, $J_{{\boldsymbol{u}}}$ is diagonal with entries given by

[TABLE]

where, to obtain the final line, we used

[TABLE]

In the arguments that follow, it will prove useful to have the product

[TABLE]

The $(d-1)\times(d-1)$ block consisting of the top $d-1$ rows of the product $J_{\boldsymbol{\phi}}J_{\boldsymbol{u}}$ is the Jacobian

[TABLE]

of the transformation from $u\in\mathbb{R}^{d-1}$ to the tuple $(p_{1},\dots,p_{d-1})$ obtained by projecting-out the first $d-1$ components of a point $p\in\Delta^{d}$ in the $d$-simplex. This Jacobian will play a role in the next section, about transforming densities.

A.2 Transforming the posterior density

Proof of Proposition 1.

We wish to transform the posterior likelihood

[TABLE]

which, up to normalisation, is the posterior density over $({\boldsymbol{p}},T)$, to the corresponding unnormalised posterior density over the unconstrained variables ${\boldsymbol{u}}$ and ${\boldsymbol{\tau}}$. The usual change of variables formula for densities says that

[TABLE]

where the final factor in the product above is the determinant of a block-diagonal Jacobian whose blocks are of two kinds

[TABLE]

Both are the Jacobians of transformations from $\mathbb{R}^{d-1}$ to tuples of coordinates such as $(p_{1},\dots,p_{d-1})$ with $p_{i}\geq 0\;\forall i$ and $\sum_{i=1}^{d-1}\leq 1$ and both have a form similar to that of $J_{{\boldsymbol{p}}{\boldsymbol{u}}}$ in Eqn. (56). In particular, both sorts of blocks are lower-triangular, and so their determinants given by the products of their diagonal elements:

[TABLE]

where $\phi_{ij}=\Lambda^{-1}(u_{ij})$ and the equality in the second line follows from Eqn. 18. Similar calculations lead to the conclusion that

[TABLE]

Combining Eqns. 57–69 leads to the conclusion that the transformed posterior likelihood is given by

[TABLE]

Taking logs and dropping terms that don’t depend on the ${\boldsymbol{u}}_{j}$ or ${\boldsymbol{\tau}}_{j}$ leads to

[TABLE]

where ${\boldsymbol{q}}_{j}({\boldsymbol{u}}_{j},{\boldsymbol{\tau}})=T({\boldsymbol{\tau}}){\boldsymbol{p}}_{j}$. Making the substitutions

[TABLE]

yields the formula in Proposition 1 and concludes the proof. ∎

A.3 Lipschitz gradients

Before turning to the log-posterior used in this paper, it’s helpful to assemble a few general facts about Lipschitz functions.

Lemma 1** (Linear combinations of functions with Lipschitz gradient).**

If all members of set of functions $f_{j}:\mathbb{R}^{n}\rightarrow\mathbb{R}$ with $j\in\{1,\dots,N\}$ have Lipschitz gradients, then so does any linear combination of the form

[TABLE]

where the $\beta_{j}\in\mathbb{R}$ are arbitrary real coefficients.

Lemma 2** (Bounded Hessians and Lipschitz gradients).**

If a continuous function $f:\mathbb{R}^{n}\rightarrow\mathbb{R}$ is twice continuously differentiable and its Hessian has bounded norm, then $f$ has a Lipschitz gradient.

Both these results follow readily from the definition of a Lipschitz function and from standard results found in, for example, Section 1.2.2 of [19].

A.3.1 Proof of Theorem 1 via two Lemmas

Proof of Theorem 1.

The negative log-posterior $U({\boldsymbol{u}},{\boldsymbol{\tau}})$ in Eqn. (28) has three kinds of terms:

those proportional to $\log\left(\Phi^{-1}_{i}(\Lambda^{-1}({\boldsymbol{u}}_{j}))\right)$, 2. 2.

those proportional to $\log\left(\Phi^{-1}_{i}(\Lambda^{-1}({\boldsymbol{\tau}}_{r}))\right)$ and 3. 3.

those proportional to $\log\left(\sum_{r=1}^{d}\Phi^{-1}_{i}(\Lambda^{-1}({\boldsymbol{\tau}}_{r}))\Phi^{-1}_{r}(\Lambda^{-1}({\boldsymbol{u}}_{j}))\right)$.

In light of Lemma 1, all we need prove is that terms of these three types have Lipschitz gradients. Those that involve either ${\boldsymbol{u}}_{j}$ or ${\boldsymbol{\tau}}_{r}$ alone turn out to be straightforward as they are covered by Lemma 3 below. The terms involving both ${\boldsymbol{u}}_{j}$ and ${\boldsymbol{\tau}}_{r}$ together require somewhat more involved calculations, but can be dispatched as is shown in the proof Lemma 4 below. Thus we are able to establish the first of the three properties required by Akyildiz et al.: our negative log-posterior has a Lipschitz gradient with respect to the unconstrained coordinates. ∎

Lemma 3** (Lipschitz gradient for sums of log-probabilities).**

If a function $f:\Delta_{d}\rightarrow\mathbb{R}$ has the form

[TABLE]

then the map $\hat{f}:\mathbb{R}^{d-1}\rightarrow\mathbb{R}$ defined by $\hat{f}=f\circ\Phi^{-1}\circ\Lambda^{-1}$ has a Lipschitz gradient

[TABLE]

Lemma 4** (Lipschitz gradients for terms involving both ${\boldsymbol{\tau}}_{r}$ and ${\boldsymbol{u}}_{j}$).**

Consider the set of $d\times d$ matrices whose columns lie in $\Delta_{d}$ and regard its elements as points in $\Delta_{d}^{d}$, the Cartesian product of $d$ copies of the $d$-simplex. Define a function $\hat{T}:\mathbb{R}^{(d-1)\times d}\rightarrow\Delta_{d}^{d}$ that maps $d$-tuples of logit-transformed stick-breaking coordinates to $\Delta_{d}^{d}$. If ${\boldsymbol{\tau}}=({\boldsymbol{\tau}}_{1},\dots,{\boldsymbol{\tau}}_{d})\in\mathbb{R}^{(d-1)\times d}$ is $d$-tuple of unconstrained coordinates and $T=\hat{T}({\boldsymbol{\tau}})$, then ${\boldsymbol{T}}_{\star j}$, the $j$-th column of $T$, is given by

[TABLE]

Now suppose that $f:\Delta_{d}\rightarrow\mathbb{R}$ is a function of the form

[TABLE]

Then the function $\tilde{f}:\mathbb{R}^{d\times(d-1)}\times\mathbb{R}^{d-1}\rightarrow\mathbb{R}$ defined by

[TABLE]

where $T=\hat{T}({\boldsymbol{\tau}})$ and ${\boldsymbol{p}}=\Phi^{-1}\left(\Lambda^{-1}({\boldsymbol{u}})\right)$, has Lipschitz gradients with respect to ${\boldsymbol{u}}$ and the ${\boldsymbol{\tau}}_{j}$.

A.3.2 A useful observation about $T_{ij}p_{j}/q_{i}$

As an ingredient in the proof of Lemma 4 we will need the following simple result:

Lemma 5** (Boundedness of ratios of the form $T_{ij}p_{j}/q_{i}$).**

If ${\boldsymbol{p}}$, $T$ and ${\boldsymbol{q}}$ are as in Lemma 4, then

[TABLE]

whenever $q_{i}>0$.

Proof of Lemma 5.

As ${\boldsymbol{T}}_{\star j}\in\Delta_{d}$, we know that $T_{ij}\geq 0$. Similarly, $p_{j}\geq 0$ and so

[TABLE]

If $q_{i}>0$, we can divide these inequalities by it to obtain

[TABLE]

which is the result we sought. ∎

A.3.3 Proofs of Lemmas 3 and

4

Proof of Lemma 3.

We can compute the gradient $\nabla\hat{f}$ using the chain rule and the Jacobians from Eqns. (29) and (30).

[TABLE]

Eqn. (46) gives the product $J_{\boldsymbol{\phi}}J_{\boldsymbol{u}}$ , which we can then act with from the right on $\nabla_{\boldsymbol{p}}f$ to obtain

[TABLE]

In general, the $j$-th component of $\nabla\hat{f}$ is given by

[TABLE]

We will now establish that the Hessian of $\hat{f}$ with respect to the $u_{j}$ has bounded norm. Begin by noting that $\partial\hat{f}/\partial u_{j}$ depends only on $u_{j}$, so the Hessian of $\hat{f}$, whose entries are $H_{jk}=\frac{\partial^{2}f}{\partial u_{j}\partial u_{k}}$, is diagonal. Calculations similar to those leading up to Eqn. (30) then establish that

[TABLE]

Now, $\phi_{j}(1-\phi_{j})$ is a quadratic with a maximum at $\phi_{j}=1/2$ and so the diagonal entries of the Hessian are bounded by

[TABLE]

which establishes that, for an arbitrary vector ${\boldsymbol{x}}\in\mathbb{R}^{d-1}$,

[TABLE]

As $\hat{f}$ is continuous and twice continuously differentiable, Eqn. (84), along with Lemma 2, establishes that $\nabla\hat{f}$ is Lipschitz with constant

[TABLE]

∎

Proof of Lemma 4.

We will need to refer to variously transformed versions of the unconstrained coordinates ${\boldsymbol{u}}\in\mathbb{R}^{d-1}$ and ${\boldsymbol{\tau}}_{k}\in\mathbb{R}^{d-1}$, so we introduce the following notation:

[TABLE]

where ${\boldsymbol{p}}\in\Delta_{d}$ is a point in the $d$-simplex and ${\boldsymbol{\phi}}\in[0,1]^{d-1}$ is the associated set of stick-breaking coordinates, while

[TABLE]

where ${\boldsymbol{T}}_{\star k}\in\Delta_{d}$ is both the $k$-th column of $T$ and a point in the $d$-simplex while ${\boldsymbol{\psi}}_{k}\in[0,1]^{d-1}$ is the associated a set of stick-breaking coordinates. As in the note transmission model in Eqn. (4), we will write

[TABLE]

and note that as we are interested here in values of ${\boldsymbol{p}}$, $T$ and ${\boldsymbol{q}}$ associated with finite unconstrained coordinates, we can be sure that, for all $j,k\in\{1,\dots,d\}$

[TABLE]

With these definitions, calculations similar to those in Eqns. (46)–(83) lead to

[TABLE]

and

[TABLE]

The components of the gradient of $\tilde{f}$ with respect to the unconstrained variables are thus

[TABLE]

and

[TABLE]

Our goal now is to show that this gradient is Lipschitz using Lemma 2. We’ll do this by computing the second partial derivatives of $\tilde{f}$ and demonstrating that they, and thus all entries in $\tilde{f}$’s Hessian, are bounded. This, in turn, will imply that the Hessian has bounded norm, completing the proof.

Consider first the terms $\phi_{j}T_{ir}p_{r}/q_{i}$ that appear in Eqn. (91). We will prove that they have bounded partial derivatives with respect to the unconstrained coordinates.

[TABLE]

The first of the terms at right in Eqn. (93) is

[TABLE]

We argued in the proof of Lemma 3 that $0\leq\phi_{j}(1-\phi_{j})\leq 1/4$ and we know from Eqn. (85) and Lemma 5 that $T_{irp_{r}}p_{r}/q_{i}\leq 1$ and so the whole term is bounded by 1.

The next term in Eqn. (93) is

[TABLE]

where we have used Eqn. (46). Here too arguments depending on the boundedness of $\phi_{k}$, and Lemma 5 establish that this term too is bounded by 1. Finally we turn to the third term in Eqn. (93)

[TABLE]

where I have used Eqn. (86). Arguments similar to those used for the previous two terms establish that this term is bounded by $(d+1)$, and thus so is the second partial derivative in Eqn. (93). Essentially the same reasoning also proves that

[TABLE]

and this observation, along with Eqn. (91), establishes that

[TABLE]

To complete our proof of the boundedness of the entries in $\tilde{f}$’s Hessain, we need to consider

[TABLE]

The analogues of Eqns. (94)–(96) are

[TABLE]

[TABLE]

and

[TABLE]

and, as above, consideration of the boundedness of stick-breaking coordinates, Eqn. (85) and Lemma 5 establish the boundedness of the derivative in Eqn. (98). Thus all entries in the Hessian of $\tilde{f}$ are bounded and so, using Lemma 2, we can conclude that the terms in $U({\boldsymbol{u}},{\boldsymbol{\tau}})$ that involve both ${\boldsymbol{\tau}}_{r}$ and ${\boldsymbol{u}}_{j}$ have Lipschitz gradients. ∎

A.4 Strong convexity

Proof of Theorem 2.

Lemmas 6 and 7 below establish that none of the three kinds of terms (see the proof of Theorem 1) that contribute to $U({\boldsymbol{u}},{\boldsymbol{\tau}})$ are strongly convex and so neither is their sum. ∎

Lemma 6** (Linear combinations of log-probabilities are not strongly convex).**

If a function $f:\Delta_{d}\rightarrow\mathbb{R}$ has the form

[TABLE]

then the map $\hat{f}:\mathbb{R}^{d-1}\rightarrow\mathbb{R}$ defined by $\hat{f}=f\circ\Phi^{-1}\circ\Lambda^{-1}$ is such that for every real number $\alpha>0$ and displacement ${\boldsymbol{\delta}}\in\mathbb{R}^{d-1}$, there exists ${\boldsymbol{u}}\in\mathbb{R}^{d-1}$ such that

[TABLE]

Lemma 7** (Terms involving both ${\boldsymbol{\tau}}_{k}$ and ${\boldsymbol{u}}_{j}$ are not strongly convex).**

If ${\boldsymbol{p}}$, ${\boldsymbol{u}}$, ${\boldsymbol{\tau}}$, ${\boldsymbol{T}}$ and $\hat{T}$ are as in Lemma 4 and $f:\Delta_{d}\rightarrow\mathbb{R}$ is a function of the form

[TABLE]

then the function $\tilde{f}:\mathbb{R}^{(d-1)\times d}\times\mathbb{R}^{d-1}\rightarrow\mathbb{R}$ defined by

[TABLE]

where $T=\hat{T}({\boldsymbol{\tau}})$ and ${\boldsymbol{p}}=\Phi^{-1}\left(\Lambda^{-1}({\boldsymbol{u}})\right)$, is such that for every real number $\alpha>0$ and all displacements $\delta{\boldsymbol{\tau}}\in\mathbb{R}^{(d-1)\times d}$ and $\delta{\boldsymbol{u}}\in\mathbb{R}^{d-1}$, there exist ${\boldsymbol{\tau}}\in\mathbb{R}^{(d-1)\times d}$ and ${\boldsymbol{u}}\in\mathbb{R}^{d-1}$ such that

[TABLE]

Proof of Lemma 6.

Assume for contradiction that there is some $\alpha>0$ such that, for all ${\boldsymbol{u}},{\boldsymbol{\delta}}\in\mathbb{R}^{d-1}$,

[TABLE]

Eqn. (83) tells us that components of $\nabla\hat{f}$ are

[TABLE]

so that the components of $\nabla\hat{f}({\boldsymbol{u}})-\nabla\hat{f}({\boldsymbol{u}}-{\boldsymbol{\delta}})$ are

[TABLE]

Substituting this into the inner product in Eqn. (102) yields

[TABLE]

For fixed ${\boldsymbol{\delta}}$, we can make the terms in this sum arbitrarily small by sending $u_{j}$ to $-\infty$ and so the inner product itself can be made arbitrarily small, contradicting the inequality in Eqn. (102). ∎

Proof of Lemma 7.

Assume for contradiction that there is some $\alpha>0$ such that, for all ${\boldsymbol{\tau}},\delta{\boldsymbol{\tau}}\in\mathbb{R}^{(d-1)\times d}$ and ${\boldsymbol{u}},\delta{\boldsymbol{u}}\in\mathbb{R}^{d-1}$,

[TABLE]

Eqns. (91) and (92) tell us that the components of $\nabla\tilde{f}$ are

[TABLE]

and

[TABLE]

As in the proof of Lemma 6, we can make the contributions to the inner product from those terms proportional to $e^{u_{j}}/(1+e^{u_{j}})$ or $e^{\tau_{jk}}/(1+e^{\tau_{jk}})$ arbitrarily small by sending the $u_{j}$ or $\tau_{jk}$ toward $-\infty$, but the terms in the first sum, which are of the from $\beta_{i}T_{ij}p_{j}/q_{i}$, require more careful treatment.

First, note that Eqns. (18), (20) and (31) imply that

[TABLE]

and

[TABLE]

As we are interested in gradients with respect to the unconstrained coordinates, the terms we need to consider involve $T_{ij}p_{j}/q_{i}$ with $1\leq i,j<d$ and so we can write them as follows:

[TABLE]

The exponential factor $e^{\tau_{ij}+u_{j}}$ can be made arbitrarily small by sending an appropriate unconstrained coordinate—$\tau_{jk}$ or $u_{k}$—to $-\infty$. Meanwhile, the quantity in parentheses is bounded above by 1, so an argument similar to that in the proof of Lemma 6 establishes that the contributions to the inner product in Eqn. (104) from terms of the form $\beta_{i}T_{ij}p_{j}/q_{i}$ can be made arbitrarily small, no matter the values of $\delta{\boldsymbol{u}}$ and $\delta{\boldsymbol{\tau}}$. This provides the contradiction we sought—the inner product cannot be both arbitrarily small and bounded from below by a quadratic—and so completes the proof of Lemma 7. ∎

A.5 Bounded second moments

Proof of Theorem 3.

Previous work suggests that pupils learn their tutors’ songs faithfully, so we initialise the ${\boldsymbol{\tau}}_{j}=(\tau_{1j},\dots,\tau_{dj})$ with unconstrained coordinates that correspond to a constant matrix close to the the $d\times d$ identity,

[TABLE]

Once $\epsilon$ is fixed (we used $\epsilon=0.1$ in the computations that lead to Figure 9) this is a constant matrix and so the second moments are of the $\tau_{ij,0}$ are trivially bounded by the squares of their values.

We initialise the ${\boldsymbol{u}}_{j}$ by computing the unconstrained coordinates corresponding to a value ${\boldsymbol{p}}_{j}$ drawn from a Dirichlet distribution with shape parameters

[TABLE]

where the $x_{ij}$ give the tutor’s note-usage counts at aligned position $j$ and $\alpha_{p}=1$ is the parameter of the prior on ${\boldsymbol{p}}_{j}$.

As we draw the ${\boldsymbol{u}}_{j,0}=(u_{1j,0},\dots,u_{dj,0})$ independently, it is sufficient to establish that the unconstrained coordinates corresponding to a ${\boldsymbol{p}}_{j}$ drawn from a Dirichlet distribution have bounded second moments. Lemma 8 below does this, establishing that $\mathbb{E}[u_{ij,0}^{2}]$ is bounded by a quantity that depends only on $\alpha_{ij}$, and so completes the proof of Theorem 3. ∎

Lemma 8** (Second moments of unconstrained coordinates for Dirichlet ${\boldsymbol{p}}$).**

Let ${\boldsymbol{u}}=(u_{1},\dots,u_{d-1})\in\mathbb{R}^{d-1}$ be the unconstrained coordinates corresponding to ${\boldsymbol{p}}=(p_{1},\dots,p_{d})\in\Delta_{d}$. If ${\boldsymbol{p}}$ is drawn from a Dirichlet distribution with shape parameters ${\boldsymbol{\alpha}}=(\alpha_{1},\dots,\alpha_{d})$, then

[TABLE]

where

[TABLE]

Proof of Lemma 8.

If ${\boldsymbol{u}}$ and ${\boldsymbol{p}}$ are as in Lemma 8 then the distribution of ${\boldsymbol{u}}$ has density

[TABLE]

where

[TABLE]

and, using Eqn. (16),

[TABLE]

The leftmost factor in $f({\boldsymbol{u}})$ is the density of the Dirichlet distribution, while the remaining factor comes from the determinant of the Jacobians in Eqn. (46).

Rearranging the products in Eqn. (108) and using Eqn. (105) to write the density explicitly in terms of the $u_{j}$, we find

[TABLE]

Gathering factors that involve $e^{u_{j}}$ leads to

[TABLE]

where the $\gamma_{i}$ appearing in the exponents of the denominators are

[TABLE]

Each of the factors in the representation of $f({\boldsymbol{u}})$ in Eqn. (110) depends on only a single one of the $u_{j}$ and so if we define

[TABLE]

the fact that $f({\boldsymbol{u}})$ is properly normalised implies

[TABLE]

We can also calculate $I_{j}$ explicitly by introducing $\phi_{j}=e^{u_{j}}/(1+e^{u_{j}})$:

[TABLE]

where $B(p,q)$ is the Beta function.

Now express the second moment of $u_{j}$ in terms of the $I_{j}$:

[TABLE]

Using this, we can establish the our result by bounding the integral at right in the last line of Eqn. (112):

[TABLE]

Finally, by combining Eqns. (112) and (113) we obtain

[TABLE]

completing the proof of Lemma 8. ∎

Supplementary Material