Attribution for the method described on this page
The methodology described here is Takuya Takahashi's, the first author of the paper. The design of the mutation model representing accentual class merger, the formulation that embeds it in a Bayesian phylogenetic framework with spatial diffusion, and the C++ implementation described below are all his work. Ayaka Onohara took part in the collaboration as second author.
This page was written by Onohara to connect that methodology to Perspective 3 (Bayesian) of the seven perspectives used in her teaching. She is not the originator or implementer of the method itself. For citation, please refer to the original paper (Takahashi, Onohara & Ihara 2024) and the Zenodo code repository listed at the end.
1What the method is for
Borrowing is the hardest problem in linguistic phylogenetics. Accentual class merger is one of the few phenomena that sidesteps it.
When phylogenies are inferred from vocabulary, there is no way to tell whether two similar words descend from a common ancestor or were borrowed through contact. That ambiguity pulls the inferred tree away from genealogy and towards a map of contact.
Japanese pitch accent has a property that avoids this. It is accentual class merger. Several accentual classes that were distinguished at an earlier stage collapse into one in a given dialect and lose their distinction. Once lost, the distinction is not restored by borrowing words from elsewhere. Which classes have merged with which is therefore a phylogenetic signal that contact does not easily scramble.
Bimoraic nouns, for instance, have five accentual classes. In the Tokyo dialect, classes 2 and 3 have merged, and so have classes 4 and 5; this is written 1/23/45. Kyoto is 1/23/4/5, and Ibukijima is 1/2/3/4/5 — no classes merged at all. Traditional dialectology has long argued about phylogeny from the distribution of such merger states, but it had rarely been written down as a statistical model. This method formulates merger as a stochastic process and puts it inside a Bayesian phylogenetic framework.
2The data: what counts as an observation
All that is observed is the surface accent pattern and the merger state readable from it.
The scope is mainland dialects. Ryukyuan and Hokkaido are excluded. Fifteen taxa — modern dialects plus accent systems recorded in historical documents — are placed at twelve mainland locations and form the leaves of the tree. Geographical space is divided into 691 cells of roughly 20 km (0.25° longitude by 0.167° latitude).
Accent descriptions are split into six word categories \(C_1,\dots,C_6\), which differ in number of morae, part of speech and conjugation and are therefore not pooled. Writing \(l(i)\) for the number of accentual classes in \(C_i\): \(l(1)=3,\ l(2)=5,\ l(3)=2,\ l(4)=3,\ l(5)=2,\ l(6)=2\).
Two things are observed for each taxon. One is the surface accent pattern \(D_i\), the pattern actually heard, written as a sequence over H/L/M/F/R. The other is the surface merger state \(M_i\), written with slashes as in 1/23/45, saying which classes have collapsed together. Their combination \(Y_i = D_i \cap M_i\) is the data given to the model.
3Model I: why latent variables are introduced
Making the observation itself a Markov process would blow up the state space. The design that keeps it small is the heart of the model.
Each cell is assigned two latent variables per word category.
\(N_i\) is the surface merger state 1/234/5 plus the information of which class overrode the others. It is written '1','234','5', where the first character of each string is the leading class, the one that overrode the rest. '1','234','5' and '1','324','5' are the same 1/234/5 on the surface but are distinct extended merger states. Since which class swallowed which cannot be observed, this is a latent variable.
\(b_i\) is the “underlying” pattern assigned to each accentual class, taking one of the patterns empirically documented for that word category.
The surface pattern follows from these two. Within a merged block, the latent pattern of the leading class is what is observed as the surface pattern of every class in that block. Take bimoraic nouns. If the extended merger state is \(N_2 = \) '1','23','54' and the latent patterns are \(b_2 = (\mathrm{HHM},\ \mathrm{HLL},\ \mathrm{HHL},\ \mathrm{LHL},\ \mathrm{LLH})\), the surface is
In block '23' class 2 leads, so class 3 is also HLL; in block '54' class 5 leads, so class 4 is also LLH.
Two properties this design captures
First, merged classes change together. A shift of the surface from (LHH, HLL, HLL) to (LHH, LHL, LHL) is explained by a single change in the leading class's latent pattern, not by two mutations. Second, merged classes always show the same pattern. Both are definitional of merger, and neither can be expressed by the standard phylogenetic models that line up binary characters.
State space: why 196 is enough
Bimoraic nouns have five accentual classes. There are 52 ways to partition a set of five (the Bell number), but an extended merger state additionally distinguishes the leading class of each block, so each partition contributes the product of its block sizes. Counting by partition type: 5 gives \(1\times5=5\); 4+1 gives \(5\times4=20\); 3+2 gives \(10\times6=60\); 3+1+1 gives \(10\times3=30\); 2+2+1 gives \(15\times4=60\); 2+1+1+1 gives \(10\times2=20\); and 1+1+1+1+1 gives \(1\times1=1\). The total is 196. These are all the values \(N_i\) can take.
Had the surface pattern \(D_i\) been made a Markov process directly, the state space would run to thousands or tens of thousands and the matrix powers inside the pruning algorithm would not finish in practical time. The latent variables exist to solve that computational problem.
4Model II: building the tree prior out of space
The prior over phylogenies is constructed from the probability of inheriting from each cell, weighted by population size and geographical distance.
In this discrete-time model, at the start of each time step every cell inherits its entire pitch-accent system from exactly one cell. Taking some variables from one cell and the rest from another is not allowed. This guarantees that the history is represented as a phylogenetic tree rather than a network.
\(a_{ij}\) is the probability that cell \(P_i\) inherits its accent system from cell \(P_j\). \(\pi_j\) is the population size of \(P_j\), \(d_{ij}\) the great-circle distance between the two cells, and \(\sigma\) the standard deviation of the Gaussian kernel. \(K_i\) is no more than the constant that makes \(\sum_j a_{ij}=1\).
The equation says two things. Cells learn more readily from nearby cells (the Gaussian kernel), and populous cells exert more influence (the factor \(\pi_j\)). The consequence is that trees in which geographically close taxa group together score a higher prior probability. The population data are built from provincial estimates of the Heian period, distributed in space using a 1 km mesh population grid from 1995.
Storing the source cell for the past \(\tau\) steps in a \(\tau \times n\) matrix \(G\) lets the tree \(T\) over the taxa be cut out of \(G\). The prior probability of \(G\) is a product of terms of the form (2.1), so the pattern of spatial interaction among cells becomes the tree prior directly.
\(\sigma\) is not estimated
Leaving \(\sigma\) free prevented the MCMC from converging, so this analysis fixes \(\sigma = 70\,\mathrm{km}\). That is an analytical setting, not a result. Runs at other values of \(\sigma\) are reported as sensitivity analyses in the electronic supplementary material, where the tree topology is said to be broadly unchanged.
5Model III: splitting mutation into two kinds
Merger (\(q\)) and pattern substitution (\(\alpha\)) are written as two independent processes.
After inheriting, each cell may mutate. Mutations are assumed independent across word categories: the evolution of \(N_i,b_i\) is independent of \(N_j,b_j\) for \(i \neq j\).
\(u_i\) is the number of blocks in \(N_i\). Every pair of blocks is assumed to merge with the same probability \(q\), so the number of available pairs \(\binom{u_i}{2}\) multiplies it. When two blocks merge, the strings representing them are concatenated in random order, so the leading class of one overrides the classes of the other. For example \(N_i=\) '1','23','45' may mutate into '123','45', '231','45', '145','23', '451','23', '1','2345' or '1','4523', each with probability \(q/2\).
Latent accent patterns change element by element, moving with equal probability to any other pattern empirically documented for that word category. This is a discrete-time version of the Mk model, in which a trait moves among \(k\) possible states at a common rate. Monomoraic nouns have seven documented patterns — HL, LH, HM, HH, LL, LM, MH — so the probability that a given class's latent pattern changes is \(6\alpha\).
Merger is irreversible; substitution is reversible
This asymmetry is the crux of the model. The \(q\) process merges blocks in one direction only, and blocks once joined never come apart. The \(\alpha\) process moves back and forth among patterns. The phylogenetic signal is carried mainly by the former; the latter accounts for differences among dialects that share the same merger state.
6The posterior, and how the likelihood is computed
The joint posterior is written out, and the tree likelihood is computed with Felsenstein's pruning algorithm.
\(B^A=(b_1^A,\dots,b_6^A)\) holds the latent accent patterns at the root of the tree. No classes are merged at the root, so its latent patterns coincide with its surface patterns. \(P(q)\) and \(P(\alpha)\) are both uniform. \(P(G)\) is the product of the \(a_{ij}\) of equation (2.1). \(P(b_i^{A})\) is uniform: for monomoraic nouns, three classes drawn without repetition from seven patterns give \(7\times6\times5=210\) possibilities, each with probability \(1/210\).
The term \(P(M_i \mid N_i)\) acts as a filter. It is 1 when the extended merger state \(N_i\) is consistent with the observed surface merger state \(M_i\), and 0 otherwise. The MCMC therefore only has to explore the space of \(N_i\) that do not contradict the observed \(M_i\).
Pruning: treating non-leading classes as missing data
The last term \(P(D_i \mid T,\alpha,b_i^{A},N_i)\) is computed by running Felsenstein's pruning algorithm in an irregular way. The point is that a class which is not leading under \(N_i\) does not appear on the surface — that is, the observation is missing at that leaf.
Putting 1 at a leaf where the class is not leading is how missingness is handled: it means “whatever that leaf is, nothing is contradicted”, and the usual recursion then goes through unchanged.
Implementation: the two-state transition probability is written in closed form
In Takahashi's C++ implementation, the two-state case of pruning is computed not by matrix powers but by a closed-form expression. For a discrete-time two-state Markov chain with rates \(q_{12}, q_{21}\), setting \(\lambda = q_{12}/(q_{12}+q_{21})\) gives
prob_two in pruning.cpp)It follows immediately from the eigenvalues of the transition matrix being \(1\) and \(1-q_{12}-q_{21}\). Avoiding numerical matrix powers makes the inner MCMC loop fast and robust to rounding. The implementation uses long double throughout, and explicitly allocates the dynamic-programming arrays over the \(2M-1\) nodes on the way from leaves to root.
7Time calibration and MCMC settings
Branch lengths come out in time steps and have to be mapped onto real time. One time step is taken to be 25 years, roughly a human generation. An MCMC assuming 10 years per step was also run, and the resulting phylogeny is reported to be almost unchanged.
The observation year of the taxa (the leaves) was set to AD 1950, reflecting the survey and publication dates of the descriptions of modern dialects. Accent systems recorded in historical documents are placed at their respective dates. The prior on the age of the root — the common ancestor — is uniform over AD 450 to 1100. The upper bound of AD 450 reflects the date of the split between Ryukyuan and mainland Japanese, which recent work places around the 3rd to 7th centuries; since Ryukyuan has accentual classes not found on the mainland, the common ancestor of mainland accent systems should be younger than that split.
| MCMC setting | Value |
|---|---|
| Independent chains | 3 |
| Iterations per chain | 4 × 106 |
| Discarded as burn-in | 1 × 106 |
| Sampling interval (thinning) | 103 |
| Total sample size | 9,000 |
| Gaussian kernel \(\sigma\) (fixed) | 70 km |
| Spatial cells | 691 |
Two further MCMC runs were carried out from different starting trees (different initial values of \(G\)) as a check. The MCC tree was produced with TreeAnnotator, distributed with BEAST, and visualised in FigTree.
8What could be concluded
Three clades, the homeland of the common ancestor, and the posteriors of the mutation rates.
The inferred MCC tree consists of three clades — Tokyo type, Keihan type and N-kei type (Kyushu) — which agrees with the traditional classification of pitch-accent systems. What is notable is that Tokyo-type accent, distributed in two separated regions in the west (Hiroshima, Oita) and the east (Tokyo, Hirosaki) of the archipelago, is suggested to form a monophyletic group. It grouped together even with spatial diffusion in the prior.
| Time to the most recent common ancestor | 95% credibility interval |
|---|---|
| Common ancestor of all modern dialects (root) | AD 450 – 825 |
| Tokyo type | AD 1000 – 1550 |
| Keihan type (including Ruiju Myōgishō) | AD 825 – 1100 |
| N-kei type | AD 1050 – 1700 |
The root falls between the middle Kofun and the early Heian period. As noted under “Limits” below, however, this posterior depends strongly on the prior.
The homeland of the common ancestor
Because the matrix \(G\) records which cell each lineage occupied at every past step, a posterior distribution over the geographical location of the MRCA is available. The common ancestor of all taxa most probably lies around present-day Kyoto, Osaka, Nara and Kobe, and the MRCA of the Keihan type falls in the same region. The MRCA of the N-kei type most probably lies in northern or central Kyushu. The MRCA of the Tokyo type is spread comparatively widely, with the highest peak around Kinki and a second peak around Tokyo; Tohoku and Kyushu are unlikely.
This inference carries a structural bias. Equation (2.1) rewards inheritance from populous cells, so the MRCA tends to be placed where the population was large. The cell containing Kyoto has the highest posterior probability both because it is the most populous cell in Kinki and because the accent patterns of Ruiju Myōgishō, the source closest to the root, were assigned to it.
Mutation rates
The posterior median of the merger rate \(q\) is \(3.37 \times 10^{-4}\) per annum (95% credibility interval \(2.25 \times 10^{-4}\) – \(4.86 \times 10^{-4}\)), which means that a given pair of unmerged blocks merges about once every 2,967 years (\(1/3.37\times10^{-4}\)). The posterior median of the substitution rate \(\alpha\) is \(9.89 \times 10^{-5}\) per annum (95% credibility interval \(6.97 \times 10^{-5}\) – \(1.39 \times 10^{-4}\)); for monomoraic nouns the rate of leaving a pattern is \(6\alpha\), since there are six alternatives among the seven documented, so a latent pattern changes about once every 1,685 years.
Where this conflicts with established hypotheses
The result does not sit comfortably with several hypotheses from dialectology. First, the gairin subtype (Oita, Hirosaki) was assumed by Kindaichi to be phylogenetically distant from the chūrin subtype (Tokyo, Hiroshima), but the inferred tree suggests that Tokyo-type accent is monophyletic, with gairin and chūrin appearing as paraphyletic or polyphyletic groupings that share no immediate common ancestor. Second, the prominent hypothesis deriving the chūrin type from the system recorded in the Bumōki was not supported. Third, the modern Kyoto system is suggested to descend not from the Bumōki system but from an undocumented lineage surviving from the Muromachi period with merger state 1/2/3/4/5. On this reading, the lineages leading to modern Kyoto and to the Bumōki each underwent the merger of classes 2 and 3 independently.
The result also disagrees with the east–west split of mainland dialects found in lexically based work such as the Bayesian phylogenetic analysis of Lee & Hasegawa. That suggests the accent systems and the lexicon of mainland dialects may have different histories; it does not by itself show that either is wrong.
9Limits
The dates depend strongly on the prior
As the paper itself repeatedly cautions, the posterior on the age of the common ancestor depends heavily on the prior. The interval “AD 450–825” given above holds under a uniform prior on the root over AD 450–1100, and rests on the assumption that the MRCA of mainland accent does not predate AD 450. In contrast to the tree topology, which was stable across the sensitivity analyses, the dates are a weak conclusion.
Trimoraic nouns are excluded
The number of possible extended merger states grows steeply with the number of accentual classes, so trimoraic nouns, with six or seven classes, were excluded for computational reasons. Several earlier studies reconstruct phylogenies from trimoraic nouns, and this analysis may miss the signal they carry. How much is lost is not obvious, however, since in many dialects the accent patterns of trimoraic nouns are not independent of those of bimoraic nouns.
The data are assembled from several sources
No published database records the accent patterns of Japanese dialects exhaustively under a consistent standard, so the patterns were collected from a number of different studies. Inconsistency in descriptive criteria across the dialects included is unavoidable, and the results may be biased by the choice of sources.
Independence across word categories does not hold in reality
The model assumes mutations occur independently in each word category, yet in the data the accent patterns of class 1 of bimoraic nouns and of verbs are the same in most modern dialects. Correlated change of this kind requires a double mutation under the model and could bias the tree. A sensitivity analysis using only the noun subset (\(C_1, C_2\)) was run, and the topology is reported to be broadly unchanged.
The latent variables are a heuristic
Assuming that latent accent patterns follow a Markov process is a device that makes computation possible, not a claim about the substance of language change. The paper itself describes the setting as heuristic and somewhat difficult to interpret. What this method infers is a posterior under this model, not language change itself.
10Sources and code
The methodology and the implementation are Takuya Takahashi's
The design of the mutation model, its integration into a Bayesian phylogenetic framework, and the whole of the C++ implementation described on this page are the work of Takuya Takahashi, the first author. Ayaka Onohara took part in the collaboration as second author. Please cite the original paper and the code repository below directly.
- Takahashi, T., Onohara, A., & Ihara, Y. (2024) Bayesian phylogenetic analysis of pitch-accent systems based on accentual class merger: a new method applied to Japanese dialects. Journal of Language Evolution, 8(2), 169–191. doi:10.1093/jole/lzae004 (published 22 July 2024)
- Takahashi, T., Onohara, A., & Ihara, Y. Data and code for “Bayesian phylogenetic analysis of pitch-accent systems based on accentual class merger”. Zenodo. doi:10.5281/zenodo.7983118 (licence: CC BY 4.0)
Structure of the implementation
The MCMC is written in C++, with each part of the model in its own file. The division of responsibilities maps directly onto the structure of the model.
main.cpp — the MCMC itself (the sampling loop)parameter_set.cpp — proposal and acceptance of parameter setsgenealogy.cpp — management of the matrix \(G\): initialisation by simulation, cutting out the tree (cut_tree_exact), proposing changes to \(G\) (mutate_diagram), and the prior of equation (2.1) (log_prior_from_matrix)group_integration.cpp — representation of the extended merger state \(N_i\) and merger mutationpruning.cpp — Felsenstein pruning; the two-state transition probability in closed form, with dynamic programming in long doubletree.cpp — tree structure (children and branch lengths)transmission_rate.cpp — computation of the inheritance probabilities \(a_{ij}\)lexical_class.cpp, location.cpp, id_converter.cpp, csv_input.cpp — word categories, locations, ID conversion, data input
References for the method
- Felsenstein, J. (1981) Evolutionary trees from DNA sequences: a maximum likelihood approach. Journal of Molecular Evolution, 17, 368–376. (the pruning algorithm)
- Lewis, P. O. (2001) A likelihood approach to estimating phylogeny from discrete morphological character data. Systematic Biology, 50(6), 913–925. (the Mk model)
- Bouckaert, R. et al. (2012) Mapping the origins and expansion of the Indo-European language family. Science, 337, 957–960. (spatial models in linguistic phylogenetics)
- Lee, S. & Hasegawa, T. (2011) Bayesian phylogenetic analysis supports an agricultural origin of Japonic languages. Proceedings of the Royal Society B, 278, 3662–3669.
Which perspective this method belongs to
The method treats the tree \(T\), the mutation rates \(q, \alpha\) and the root state \(B^A\) all as random variables, and builds a posterior from a prior and a likelihood. The distinguishing feature is that the tree itself carries a posterior, so competing phylogenetic hypotheses can be compared as probabilities. This is Perspective 3 (Bayesian) of Seven Perspectives in Equations. It can also be read against ICAR in Perspective 3, which likewise puts spatial structure into a prior, and against Neighbor-Net in Perspective 1, which takes the opposite route of drawing a structure from distances.
This page was written by Onohara, but the method it describes was devised and implemented by Takuya Takahashi. Every figure is quoted from the published paper; no estimated or approximate values are given.
The whole picture of where each method sits is in Seven Perspectives in Equations, and which teaching session falls under which perspective is in Methods & Mathematics.