Bayesian inference of phylogeny has great impact on evolutionary biology. It is believed that all the species are related through a history of a common descent BIBREF0, that is to say, the reason we have various wildlife, including human beings, is because of evolution. We can show the process of evolution and solve the problems, like what is the phylogeny of life, by showing a phylogenetic tree (see Figure FIGREF1).As a matter of fact, any research of DNA sequences or protein pattern taken from the species or an individual wildlife can start with a phylogenetic analysis. Language change, which is regarded as just like life evolution, can similarly borrow the phylogenetic analysis to discover the process in which how languages change. It is accepted that all human languages start from the proto-language which is the ancestor of modern languages. For example, Germanic languages (like English, German, Dutch), Romance languages (like Spanish, French and Portuguese), Slavic languages (like Russian, Bosnian and Polish) are from Indo-European languages, which was developed from proto-Indo-European languages. Therefore, we can borrow the computational method from biology and evolution to apply to the language topology.Not only Natural Language Processing (NLP), more and more linguistic branch disciplines begin to have digital datasets and use computational methods to solve some problems. Historical linguistics recently have dramatically increasing digital datasets. The availability of new data and researches of different language from different language families start to challenge the traditional way to carry on the study of historical linguistics. The comparative method has been the core method for linguistic reconstruction for the past 200 years, and is based on manually identifying systematic phonetic correspondences between many words in pair of languages BIBREF2, because different culture has different writing system and different way to spell their words. That is why International Phonetic Alphabet was created to help linguists analyze different languages in parallel. However, there are too few scholars, i.e., historical linguists, to analyze the world's over 7500 type of languages BIBREF2, BIBREF3, including thousands of languages that have not been studied and are facing extinction. Thereafter, computational methods can help people do researches on unstudied languages faster and give a possible solution. Phylogenetic inference of human languages task is composed with two parts: cognate set detection and phylogenetic tree construction. Cognate set detection automatically assists people put language words with similar or possible evolutionary patterns to one cluster. The phylogenetic tree construction task build trees given the information from the clusters. In the following, I will divided the whole paper into two main steps: the way I implement cognate detection would be discussed in section 2. After that, I will use the cluster data to carry on the phylogenetic inference program and build the relationship trees, which I will describe the way that how I finished this part in section 3. I will show the results and evaluate them in section 4, and make a conclusion in the last section 5.A great number of algorithms and mechanisms to antomatic cognate detection which could be applied in historical linguistics have been used and tested if they are working by many linguists and computer scientists BIBREF2, BIBREF4, BIBREF5, BIBREF6, BIBREF7. In detail, many of these works are very similar to each other, which consist of two main stages. For the first stage, they first extract the words with same meaning from the wordlists of different languages, either same or different language families, and compare them and use the distance calculation matrix to compute how similar they are. Regarding the second stage, a flat cluster algorithm or a network partitioning algorithm is used to partition all words into cognate sets, and also take the information in the matrix of word pairs as basis BIBREF4, BIBREF5. However, the methods in which those researchers use to compare the word pairs are totally different in that people could use different methods to pre-process their language datasets, or even use different algorithms to finish the comparison and clustering task. For example, intuitively, people will start the automated word comparison by computing the distance between the words, such as word embedding in NLP, GloVe BIBREF8, which computes the semantic similarities between the two words. In computational historical linguistics, phonetic segments are used to calculate how close the two words are instead, because the semantics of a single word is not changing as easy as phonetics is. The problem is since the program involves the data pre-processing, then the whole dataset would be traverse twice and the computation would be a big problem when the dataset is about to be over 100 languages. Consequently, people began to think about a faster method to help.The first linear time method was proposed by BIBREF9, and later modified by BIBREF6. The algorithm compares the word pairs by their consonant class. A consonant class is hereby understood as a rough partitioning of speech sounds into groups that are conveniently used by historical linguists when comparing languages BIBREF4. Table TABREF3 shows some of the international phonetic alphabet (IPA) (for more details about IPA, please refer to BIBREF10). After getting the IPA of the word pairs from the wordlist, the algorithm is to determine if they are cognate by judging their first two consonants class match each other or not. However, since the algorithm only compares the first two consonant classes, the accuracy of this algorithm is relatively low. I have two reasons for this: (a) In linguistics, the number of possible sounds in human languages in the world, excluding artificial languages, amounts to the thousands BIBREF5. It is unrealistic to enroll all the sounds into the system. If we enroll all the sounds in the algorithm to simulate the language change process, we need to get a new matrix in which the probabilities of one sound switching to the other sound will be calculated, which is very time-consuming. (b) comparing the first two consonant classes are not sufficient to determine if the two words in pair are derived from the same cognate word.The Edit Distance approach is to take the normalized Levenshtein distance BIBREF11, which is a concept in information theory. It aims to measure the difference of two sequences by calculating the minimum number of character or string edits, such as insertion, deletion, which are coincidentally two basic language phonetic change. The distance could be used as a probability to estimate how possible one word changes to the other one.This algorithm is for pairwise and multiple alignment analysis BIBREF2. It not only takes an expanded sound class into account, but it considers the prosodic aspect of each word. As a result, it is able to align within morpheme boundaries instead of the sound segments, suppose the morpheme information is the prior knowledge and we already have it.The previous three methods use the same strategy to put the words from different languages into clusters, i.e., UPGMA clustering algorithm, while LexStat uses language-specific scoring schemes which are derived from a Monte-Carlo permutation of the data BIBREF5. The word pairs from different languages are first aligned and scored, and the MC permutation shuffled the word pairs according to their scores. The scores could be calculated by the frequencies of daily use by native speakers. Thus, a distribution of sound-correspondence frequencies is generated. Then, the distribution is used to compare with an attested distribution and then converted into a language-specific scoring scheme for all word pairs.Following the algorithms above, with the consideration of both the advantages and disadvantages of them, in this project, I am going to use a modified method: sound-class based skip-grams with bipartite networks (BipSkip). The whole procedure is quite straightforward and could be divided into three steps. First step: the word pair and their skip-grams are two sets of the bipartite networks. The second step is optional, which is to refine the bipartite network. Before I run the program, I will be asked to input a threshold, which determines if the program should delete the skip-gram nodes linked to fewer word nodes than the threshold itself. According to the experiment, even though I did not input any threshold as one of the parameters, the algorithm could still give the same answer but with more executing time. In the last step, the final generated bipartite graph would be connected to a monopartite graph and partitioned into cognate sets with the help of graph partitioning algorithms. Here I would use Informap algorithm BIBREF12. To make a comparison to this method, I am using CCM and SCA for distance measurement in this experiment, too. UPGMA algorithm would be used accordingly in these two cases.Methods for Bayesian phylogenetic inference in evolutionary biology and historical linguistics are all based on the following Bayes rule BIBREF4, BIBREF13:orWhere $f$ means the probability density function, $\Lambda $ consists of the tree topology $\tau $, the branch length vector of the tree $T$ and the substitution model parameter $\theta $; $X$ is a data matrix with dimension $N*K$, within which there are $N$ rows indicating $N$ different words from different kinds of languages and they can be put into $K$ cluster in a language family. Figure FIGREF8 shows an example of the matrix. As we can see, the data matrix is a binary matrix with every element $i_{ij}$. 1 means language $i$ could be classified into cognate $j$, while 0 means language $i$ does not belong to cluster $j$. Based on the shape of the data matrix, to get the parameter ($\Lambda = \tau , X, \theta $) for tree generation, we need to sum over the whole matrix and will have to compute all possible topologies of $\frac{(2N-3)!}{2^{N-2}(N-2)!}$.This method is practical only for a small amount of dataset and the calculation is not easy once it is over 5 languages BIBREF13. In addition, based on the Bayesian formula, if we use prior probability $Pr(tree)$ and likelihood $Pr(Data|Tree)$ to produce the posterior probability, it seems that the posterior probability is very easy and convenient to formulate based on the Bayesian formula, but to compute this probability and get the greatest estimate of the final tree, the machine has to compute and compare all the possible trees and in each tree, the machine will compute all the combination of branches with different length.Metropolis-Hastings (MH) algorithm BIBREF14 , one of the Markov Chain Monte Carlo techniques, would be a good tool to overcome this computation difficulty by avoiding summing over all of the topologies by evaluating the posterior probability $f(\Lambda |X)$. The likelihood in this case from one data to the next parameter is calculated by the prunning algorithm. Prunning algorithm, or we can call it K81 model BIBREF15, is a Markov model of DNA evolution, and this model is a description of the DNA in evolution as a string of four discrete state, i.e., G, C, A, T. Fortunately, language model is very similar to DNA model in that both of them are discrete models , in the language model, we only apply this model in the binary dataset.MCMC method is helpful here since it generates a probability distribution $\pi =\lbrace \pi _{i}\rbrace , i=0,1,...,$. We construct a Markov chain with the stationary distribution being posterior probability distribution of the parameters. A new tree will be proposed by stochastically perturbing the current tree and the new tree would be accepted or rejected. If the new tree is accepted, then repeated the algorithm, which is subjected to the future perturbation. If this Markov chain is constructed and configured properly, the proportion of the time that any of the trees is visited is a valid approximation of the posterior probability of the tree BIBREF16. Besides, $\pi _{i}$ is sometimes still not easy to calculate, but we can get the function $\pi = \frac{\pi _{j}}{\pi _{i}}$ directly instead. MH algorithm is the method that could help us solve this function. This algorithm tells us that given sample states $i,j\in \Omega $, and every state $i\in \Omega $ has a distribution $\pi =\lbrace \pi _{i}\rbrace $. There is a potential transition from state $i$ to $j$ and the transition probability is $g_{ij}$, while the acceptance transition probability, meaning a move from state $i$ to $j$, is $\alpha _{ij}$. We can easily get the properties as shown below:Therefore, the chain $\textbf {P}=\lbrace p_{ij}\rbrace $ has $\pi $ as its stationary distribution and it satisfiesA sufficient solution for this chain to generate $\pi $ as its stationary distribution is that $\lbrace g_{ij}\rbrace $ be irreducible and aperiodic BIBREF13, BIBREF17. MH algorithm aims to sample parameters from the posterior distribution, which could help us generate the posterior distribution of phylogenetic trees and each state of this Markov Chain is probably labelled as cognate set. We can get the simple relationship as $\pi _{i}=f(\tau =i|X)$. To calculate $\pi = \frac{\pi _{j}}{\pi _{i}}$, we have ratio as shown below:We can have $\alpha _{ij}$We put $\alpha _{ij}$ the method and have algorithm 1. Under this method, the Markov chain can finally work and the time is shortened greatly, but due to the large scale of the dataset, the new problem is the chain has a large probability that it can get stuck in a local maxima if the posterior probability density function has multiple peaks.To solve this problem, BIBREF18 proposed Metropolis-coupled Markov Chain Monte-Carlo methods (MC3) aiming to solve the problem when Markov chain has to move across a distribution with multiple peaks, which requires $n$ chains run simultaneously. In the iterations, we can get $n$ different stationary distributions $\pi _{n}$, but only the first one $\pi _{1}$ would be the target density and the rest of them would be treated for improving mixing. To 'heat' some of these chains, the posterior probability of the corresponding chains are raised to the power of $\beta $. For instance, if the probability density function is $f(\Lambda |X)$, then the heated distribution is $f(\Lambda |X)^{\beta }$, where $\beta \in (0,1)$. Once we get a heated version of Markov chain, this technique actually make new states easily get accepted. In comparison, we can have the heated version of the chain:In the ratio, if $f(\Lambda |X) > f(\Lambda ^{\prime }|x)$, this will trigger the increase of each state to the power of $\beta $. As a result, the heated chain is more likely to accept new states than other cold chains. The MC3 algorithm needs very expensive computation, since it has to work with multiple CPU cores to run at the same time.Based on the information above, to avoid both shortcomings with cheaper computation, I am going to use MH algorithm with simulated annealing BIBREF19. This method is one of the modified version of MH algorithm with change of the acceptance ratio, similar to MC3 but only with single CPU. We are going to use the equation, where $T_{i}$ represents temperature.In this method, I will have an initial temperature $T_{0}$, which was set very high at the beginning, indicating heated Markov chain and could be decreased in the following steps until it is set to be $T_{i}\rightarrow 0$. Regarding the change of $T_{i}$ value, i.e., to cool down the Markov chain, a cool-down algorithm is inserted in the method. According to BIBREF19, a linear cool-down algorithm in which the $T_{0}$ initial value was extremely high at the very beginning and the algorithm decreases it in every iteration in $T_{i} = \lambda T_{i-1}$ with $\lambda $ being set between $[0.85,0.96]$ BIBREF20. In the experiment, I reduce the temperature to cool down the chain until $T_{i} = 10^{-5}$. The final state of this whole Markov chain is regarded to be the maximuma-posterior (MAP) estimate of the inference. The Bayesian analysis in this project uses binary datasets with all states 0 and 1. The Generalized Time Reversible Model (GTR model) BIBREF21, BIBREF22 is employed to compute the transition probabilities between states. When I build the phylogenetic trees, the tree branch length and also the shape of the tree are continuous variables, a move satisfying exponential distribution would be used as the potential transition function. When the program is launched, a uniform move will randomly choose two states in $\Omega $ and propose a new frequency with the summation of the frequencies of the states that are not changed. I use the Subprunning and Regrafting move and Nearest Neighbor Interchange BIBREF23 to operate on the nodes and the leaves to build the new trees.It is not easy to figure out which kind of skip-gram and sound-class system would generate the satisfying results, so I design my experiment as follows. (1) During the training process, I would use datasets proposed by BIBREF7, I will give the training dataset in table TABREF14, and compute skip-gram by length 4. Although there are languages repeated in both training language dataset and testing language dataset, such as Chinese in Sino-Tibet and Austronesian, I manually tick out the repeated data from both dataset and also the concept from both datasets are not overlapping each other. (2) Regarding sound-class, I would refer to SCA sound class BIBREF7 (3) Set the threshold in step 2 as 0.2 (20%). (3) The evaluation matrix is B-Cubes BIBREF24. The F-score based on the configuration above is, when I use connected components partitioning, 85.4%, and 85.2% when I use Infomap partitioning algorithm.To evaluate the cognate detection, I choose to use B-Cubed algorithm proposed by BIBREF25. The algorithm is completely based on clusters in way of precision, recall and F1-score. The precision and recall score for each cluster are calculated separately and the final recall and precision score for the entire output are from the combination of each cluster/entity in the output. This evalution matrix has been implemented by LingPy.To test if the phylogenetic tree generation module is working good or not, I utilize Generalized Quartet Distance (GQD) matrix, which compares the tree from historical linguistic experts and the one generated by machine. A quartet distance is used to measure how similar the two phylogenetic trees are, which is defined to count the number of quartet that differ between two trees. A quartet is defined as a set of four leaves selected from a set of leaves without replacement BIBREF26, i.e., suppose I have $n$ leaves in the tree, then I need to compare $\binom{n}{4}$ pairs of quartets in both trees. The only 4 topologies are shown in figure FIGREF16.Given a tree $\tau $ with $n$ leaves, I are able to divide $\tau $ into sets of stars $S(\tau )$ and sets of butterflies $B(\tau )$. Now I can define the similarities between the two trees GQD as follows, where tree $\tau _{g}$ is the ground truth tree from historical linguists theoretically.Practically, I used golden standard tree from Glottolog BIBREF28. Glottolog provides a detailed list about understudied languages, including their language families. This dataset collects 7,592 understudied spoken L1 language. Linguists began to use this dataset to build phylogenetic trees with higher accuracy, since the cognate sets are annotated manually. For example, BIBREF4 uses this as a ground truth to compare with the phylogenetic tree constructed with the cognate sets generated automatically. Their result shows their tree has higher quality if I just use annotated cognate sets, which support the motivation of the automated cognate detection of the unstudied language, helping linguists to study the less studied language more efficiently.Cognate Detection Table TABREF19 shows the statistics in test data developed by BIBREF2. The result of the BipSkip approach for cognate detection is shown in the table TABREF20 as well as the result of using SCA and CCM. As shown in the tables, we can see that the BipSkip approach is not the quickest method, although it is more user-friendly. CCM is surprsingly fastest method with slightly higher precision than BipSkip approach, especially in Austronesian and Sino-Tibetan. Indo-European languages stays almost the same, my guess is because the languages in Indo-European language family are more phonetically similar to each other than other languages in their corresponding language families, as a result of which the three methods would perform almost the same. SCA algorithm is not recommended here in that it costs the longest time and the largest spece since I need to prepare the expanded sound class and also some morphological features.Phylogenetic Inference The result of phylogenetic inference in modified MH algorithm is shown in table TABREF22. I designed a branch of experiments, changing the settings and get some good results. I set the initial temperature as $T_{0} = \lbrace 10, 20, 30, 40, 50, 60, 70, 80, 90, 100\rbrace $. During the project, I will take down the number of iteration, the time it takes to finish running the program. Table TABREF21 is the ground truth results, testing on golden standard cognates from Glottolog. It is hard to determine which one is outperforming the other two. Overall, Indo-European costs the shortest time and fewer iterations, since in the three methods this language family always has the highest the accuracy. In addition, the cognate sets from the automatic approaches is easier to build phylogenetic trees than the manually annotated standard cognate sets, from the result, the automatic methods obviously shorten the time of automatic phylogenetic inference.Obviously, the result is not surprisingly good. We can observe from the data, the accuracy for each language, some of them are only slightly over 50%. Among the five language families in the testing data, Indo-european has more accuracy than the other four language families, due to the similar phonetic features. Also, some places where native people used to live and was conquered by immigrates, for example languages in the islands in south pacific, Hawaii, Indonesia, etc, their accuracy is obviously high and easy to cluster by cognate and find their relationship. Some native languages in Australia, Pama-Nyungan language family whose languages are mainly on Autralian continent is surprisingly lower than any other southern pacific islands languages.From this exmperiment, we can obviously use this method to help historcial linguists to make an analysis of language development and change, but the result is not very accurate basically. How do we solve the problem? The current datasets only includes the main class of phonetic alphabet, I think it is necessary to enroll some language phonetic change background knowledge to let the machine recognize the probability of change from phoneme $A$ to $B$, such as Great Vowel Shift, etc.In this appendix, I am going to show the part of resulting tree generated from Indo-european language testing dataset. The whole tree structure is very big which involves 34 trees totally. I am only showning the first four trees. The number on the edge is the probability that they are related. The labels at the end of the tree are language type: such as poli represents polish, norw represents Norwegian; lati represents Latin. The number behind the language type is the index of the word in the wordlist.

Q: Is the proposed method compared to previous methods?

Q: What metrics are used to evaluate results?