text:
# Error in entropy properties in Mathematical Theory of Cryptography by Claude E. Shannon

I am reading this classic paper by Claude E. Shannon and I think there may be a couple of errors in his description of the properties of Entropy/Uncertainty. The screenshot shown at the bottom of this post is from page 9. Firstly, I believe there is a typo at the first red circle which should read $$10^{50}$$. Secondly, I think the equation shown for the value of the entropy $$H$$ of this problem is also incorrect. I believe the correct equation for $$H$$ is derived as follows:

events $$=\{p_1, p_2, ... , p_{10^{30}}, q_1, q_2, ..., q_{10^{50}}\}$$

where events $$p_i$$ each have a probability $$\frac{1}{2}(10^{-30})$$ of occurring, and events $$q_j$$ each have a probability $$\frac{1}{2}(10^{-50})$$ of occuring

\begin{align} H &= -\sum_{i=1}^{10^{30}} p_i\log{p_i} - \sum_{j=1}^{10^{50}} q_j\log{q_j} \\ &= -\sum_{i=1}^{10^{30}} \frac{1}{2}(10^{-30})\log{(\frac{1}{2}(10^{-30}))} -\sum_{j=1}^{10^{50}} \frac{1}{2}(10^{-50})\log{(\frac{1}{2}(10^{-50}))} \\ &= -\frac{1}{2}\log{(\frac{1}{2}(10^{-30}))} -\frac{1}{2}\log{(\frac{1}{2}(10^{-50}))} \\ &\approx 40 \end{align}

Hence, I think the equation for $$H$$ should read:

$$H = -\frac{1}{2}\log{\frac{1}{2}10^{-30}} -\frac{1}{2}\log{\frac{1}{2}10^{-50}}$$

I would really appreciate if someone could let me know if I'm correct or where I'm wrong. I am not doubting Claude, just wondering if the document itself has typos. Thanks!

• Yes, you are right, these are typos. Btw, you can simplify $-\frac12\log(\frac1210^{-30})-\frac12\log(\frac1210^{-50})$ to $\frac12\log(2\cdot10^{30})+\frac12\log(2\cdot10^{50})=40+\log2$ (I assume from the context these are base-10 logarithms). – Emil Jeřábek Mar 22 at 9:57
• Thanks, yes they're base-10. That's a cool simplification. – nuggimane Mar 22 at 11:01


text:
Change the chapter
Question
Hot air rises because it has expanded. It then displaces a greater volume of cold air, which increases the buoyant force on it.
1. Calculate the ratio of the buoyant force to the weight of $50.0^\circ\textrm{C}$ air surrounded by $20.0^\circ\textrm{C}$ air.
2. What energy is needed to cause $1.00 \textrm{ m}^3$ of air to go from $20.0^\circ\textrm{C}$ to $50.0^\circ\textrm{C}$?
3. What gravitational potential energy is gained by this volume of air if it rises 1.00 m? Will this cause a significant cooling of the air?
1. $1.10$
2. $2.79 \times 10^4 \textrm{ J}$
3. $12.7 \textrm{ J}$
Solution Video

# OpenStax College Physics for AP® Courses Solution, Chapter 14, Problem 75 (Problems & Exercises) (10:00)

#### Sign up to view this solution video!

Rating

No votes have been submitted yet.

## Calculator Screenshots

Video Transcript
This is College Physics Answers with Shaun Dychko. Some 50 degree air is surrounded by 20 degree air. And, our question is to find the ratio of the buoyant force on the hot air to its weight, the weight of the hot air. So, the buoyant force is equal to the weight of the fluid displaced by the thing that is submerged. So, you can imagine that this hot air is submerged so to speak in this cold air. And so, we need to find the weight of the cold air that was displaced by this hot air. So, that'll be the mass of the cold air displaced multiplied by G. Acceleration due to gravity or gravitational field strength. And, now we don't know what the mass is, but we can figure out the volume because we're given the volume of the 50 degree air, and that will be the same as the volume displaced. It'll displace an amount of 20 degree air equal in volume to the size of the thing, the thing being this hot air. So, the density of the 20 degree air is the mass of the 20 degree air divided by its volume. I did not put a subscript on the volume for the 20 degree air because the volume of the 20 degree air displaced is the same as the volume of the 50 degree air and so there's no need for a subscript. Now, we'll solve this for M 20, the mass of the 20 degree air by multiplying both sides by V, and that gives us this expression which we can substitute in for M 20. So, the weight of the fluid displaced by the hot air, and I mean fluid in a more general sense to include air in gases. So, the weight of the 20 degree air displaced is the density of the 20 degree air multiplied by its volume times G. So, this can all be written in place of W in this Archimedes' Principle formula to say that the buoyant force equals density of 20 degree air times its volume times G. Now, we want to take this buoyant force and divide it by the weight of the 50 degree air. So, let's find out the weight of the 50 degree air now. And, with the same sort of reasoning up here, it's going to be the density of 50 degree air times its volume times G. And, now we divide these two things. So, divide this by this and we get this expression and the volumes and G's cancel. And, we're left with the density of the 20 degree air divided by the density of the 50 degree air. And, now we don't know what these densities are but we do know temperatures. So, it turns out that this ratio is going to be the ratio of the temperatures. So, we can use this ideal gas law, which says pressure times volume equals the number of moles of the gas times universal gas constant R times its absolute temperature. And, we'll solve this for V and we'll substitute that into our density formula. Now this is density of 20 degree or 50 degree air. It is just density in general for a gas. Now, pressure we can assume is going to be constant in this question because it's all in the atmosphere here and if there is much pressure difference, then there will be wind blowing and high pressure would move towards low pressure and so for a situation like this to sustain for a little while, it's going to have to not have any pressure differences. So, assume pressure is constant. I mean, there are lots of assumptions in this question, really. Even assuming that G is constant is a bit of an assumption because Nah, never mind. That's a pretty good assumption because they're both at the same distance from the center of the earth. Never mind. Okay. So, let's not digress too much. Density is mass divided by volume. So, let's divide by this. Because we want to express things in terms of temperature since that's the information we know, we can't go look up the density of air in a data table because density of air, which is about 1.29 kilograms per cubic meter, this value is at standard temperature and pressure. And, these are not standard temperatures. This is 20 degrees Celsius and this is 50 degrees Celsius. Standard temperature is 0 degree Celsius. So, we can't resort to that. And so, we have to think of something different, which is to consider that density is mass divided by volume. And, we'll divide by this, which is the same as multiplying by its reciprocal. So, we have mass times pressure over number of moles times universal gas constant times temperature. And so, we can express density in this way in terms of T. Okay. So then, the mass, also, we want to get rid of because we don't know what that is either and it's going to be the number of moles of the gas times the molar mass of the gas. And, this molar mass is kind of a chemistry concept but think of it as the number of grams in a mole of this gas mixture, which will contain nitrogen, oxygen, little bit of CO2 and so on. Okay. So, we plug that in for M and we get this. And that's useful because now these number of moles cancels and we didn't know what that was anyway, so that's good that it cancels. And then we have all of this. Molar mass times pressure divided by the universal gas constant times temperature. So, now we can rewrite this buoyant force to weight expression as the ratio of densities, but instead of densities, we're going to substitute this. So, we have the density for the 20 degree air is molar mass times pressure divided by gas constant times the temperature of this one, the 20 degree air. And then, we'll divide that by the density of the 50 degree air but using this expression in place of density. Now, dividing by this fraction is the same as multiplying by its reciprocal. So, we're going to write R T50 over M P. And then a whole bunch of things cancel. And, we're left with T50 over T20. So, the ratio of the buoyant force to the weight, then, is going to be 50 degrees divided by 20 degrees but expressing both of those in absolute temperatures, by adding 273.15. And, we get 1.10 is the ratio of the buoyant force to the weight. Alright. Part B is asking what energy is needed to cause a cubic meter of this area to go from 20 to 50 degrees Celsius. So, the energy that will be absorbed is going to be the mass of the air times the specific heat of air times the change in temperature. And, we don't know what the mass is. We do know volume, and so we will express all of this in terms of density. So, density is mass divided by volume. And, we'll solve for M by multiplying both sides by V and then we will substitute this in place of M. So, the amount of energy that will be consumed to increase the temperature of the air will be its volume times the density times the specific heat times change in temperature. And, we are using this 1.29 kilograms per cubic meter because that's the best estimate we can get in this circumstance for the density. And, the actual density will be somewhere close to this, anyway. And so, we put in 1.29 there times a cubic meter times 721 Joules per kilogram per Celsius degree times the change in temperature, which is the difference between 50 and 20, and this gives 2.79 times 10 to the 4 Joules of thermal energy required to increase a cubic meter of air in temperature from 20 to 50. The amount of gravitational potential energy that a cubic meter of air would gain as it increases in height of one meter will be its mass times gravitational field strength times the change in height. And, mass is volume times density, as we said up here. And, we'll substitute that. And so, we have a cubic meter times 1.29 kilograms per cubic meter density times 9.81 newtons per kilogram gravitational field strength times one meter increase in height, which is 12.7 Joules. Now, this energy must come from somewhere because energy cannot be just created out of nothing. That's the conservation of energy law for which no exception has ever been noticed ever. So, this gravitational potential energy is going to come from thermal energy of the air because the air is going upwards because of temperature differences which are causing density differences and thereby causing their buoyant force. So, it turns out that it would not take very much thermal energy to impart this amount of gravitational potential energy because 12.7 Joules is much much less than 2.79 times 10 to the 4 Joules. And, there we go.

format: plain


text:
# Print R Matrices in LaTeX Math Environments

This is a quick and easy method to print R matrices in LaTeX math environments. It’s based on the R function found here. (The newline at the end of the bmatrix should be removed)

Use this function

bmatrix = function(x, digits=NULL, ...) {
library(xtable)
default_args = list(include.colnames=FALSE, only.contents=TRUE,
include.rownames=FALSE, hline.after=NULL, comment=FALSE,
print.results=FALSE)
passed_args = list(...)
calling_args = c(list(x=xtable(x, digits=digits)),
c(passed_args,
default_args[setdiff(names(default_args), names(passed_args))]))
return(cat("\\begin{bmatrix}\n",
do.call(print.xtable, calling_args),
"\\end{bmatrix}"))
}

And then insert the matrix as seen here

$$a = <<results = 'asis', echo=FALSE>>= bmatrix(a) @$$

Grab the .Rnw example file here. It should compile into something like this

Last updated 2014-09-09


text:
## Featured resource

### Numeracy: What Classroom Teachers Should Know

What is numeracy in an Australian context?

Members: $22.00 inc.GST Others:$ 27.50 inc.GST

Home > Make It Count 2 > Resources > Significant episodes

# Significant episodes

A Significant Episode is an event, or a small series of linked events, that has impacted on educators professionally and has helped develop or change the way they teach and think about teaching. These Significant Episodes provided rich data for the Make It Count Cluster Findings.

Below is a list of all of the Significant Episodes from the eight Clusters. Each story is also provided in PDF format for download and has focus questions to prompt discussions amongst educators.

## 8ways learning map

“In bestowing this name, it gave authenticity to our Make It Count project and enhanced community engagement.”

## Accelerating learning

"Each boy has had individual conversations with their teacher and both admitted that whilst they could have and should have worked harder first semester, they have both given a commitment to achieve a higher standard this semester."

## A change in how I taught

“It was scary to start doing that. It was completely different to any pedagogy or methodology I’d done before. On the other hand, I wanted to do something different because I found what I was doing with my 5-6-7 class last year wasn’t successful.”

## A shift from comfort

“What I found over time was that the development of this strong sense of number and an understanding of how to manipulate its parts, may have been the stumbling block for many of the middle primary students I used to teach."

## Big day out: Culture and maths come together I

"...think about making connections between maths and culture."

## Big day out: Culture and maths come together II

"I believe that with Indigenous ways learning, everything is tied together rather than the compartmentalised Western concept which puts maths here, and English over there..."

## Teacher-student connection

“It was not just the quality of the questions that impacted on learning but the spread of questions to children in the class and the types of questions being asked.”

## Chocolate cake

"The key outcome for me was the power of a meaningful and conceptually relevant context for investigating mathematical ideas."

## Collaborative planning

“Collaborative planning needs to be a focus for the future.”

## Commitment

"Students and parents gave positive feedback about the impact the project has had..."

## Education assistants appreciate maths training

“It was an aha moment when we realised how important it is that children can subitise and where partitioning belongs in the concept of place value.”

## Estimation Walks (a maths300 lesson)

“I never used to estimate. Now I estimate how big my brother’s room is, then I measure it.”

## Explicit learning intentions

“Making learning intentions and success criteria explicit helps Aboriginal students to better understand the expectations of the task and when they are successful with the task.”

## Aboriginal educators take First Steps

“One of the outcomes is the critical building of connections between the school and community...and the development of Aboriginal educators as leaders within the school and the community.”

“I let the class into the room and this boy who normally chooses to sit at the back with talkers/non-workers decided to sit front and centre, by himself.”

## Hands on activities I

“I believe that the hands-on approach has seen Indigenous students (and other students) engage with mathematics in a more positive way.”

## Hands on activities II

"We played for a bit longer, modelling the story then writing the equation and she got it."

## His confidence continues to grow

“Almost immediately Student A began to answer questions and/or give an explanation on how or why he had arrived at his particular solution.”

## How I look at maths

"Following the learning sequence has made a big difference on how I look at maths."

## How can it not be the best pedagogy?

"I have seen significant change in teaching."

## I barely had to say a word!

"I will never forget the light bulb I saw appear above her head! And I barely had to say a word."

“Devoting specific time in department meetings to focus on teaching and to share lesson ideas.”

## Leading from the middle I

“A simple yet powerful shift in the style of teaching now includes a questioning technique... as a way to properly and effectively engage students.”...

## Leading from the middle II

“During the implementation, staff followed an ‘I do’, ‘We do’, ‘You do’ philosophy to demonstrate the process in a variety of contexts to encourage students to at least attempt the questions.”


text:
In mathematics outside of school we use radians instead of degrees for convenience. I have found that many undergraduate students don’t understand what a radian is.

An angle of 1 radian is created when the length of the arc created by the angle is equal to the length of the radius. Thus, a radian is the angle that occurs when the arc and radius are equal in length. This can be seen in the image below, where we have an angle of 1 radian, and the arc and radius are of the same length.

Speaking more mathematically, if we have a radius length denoted $$r$$, then an angle of 1 radian occurs when that angle creates an arc length equal to $$r$$. Expanding upon this notion, we discover that an angle of 2 radians gives an arc length that is twice the size of the radius. By increasing the angle, we get a larger arc length and subsequently a larger radian value because the arc is getting larger in comparison to the radius which is always fixed. On the other hand, by decreasing the angle, we get a smaller arc length and subsequently a smaller radian value because the arc is getting smaller in comparison to the fixed radius. From this we can say that the arc length, $$\alpha$$, is equal to the angle in radians, $$\theta$$, multiplied by the radius, $$r$$.

But why is a 360 degree angle equivalent to $$2\pi$$ radians? Note that the circumference of a circle is equal to $$2\pi r$$. This means that the length of the arc of a full circle is $$2\pi r$$. That is “two multiplied by the value of pi multiplied by the value of the radius”. But since we know that the arc length is equal to the angle in radians multiplied by the radius, we can do some algebra to find what the angle will be when the arc length is equal to $$2\pi r$$.

As can be seen from above, we have proven that an angle of 360 degrees is equivalent to $$2\pi$$ radians and also gained some understanding and intuition about what radians really are.


text:
Learning Poisson Binomial Distributions

# Learning Poisson Binomial Distributions

MIT
costis@csail.mit.edu
Research supported by a Sloan Foundation Fellowship, a Microsoft Research Faculty Fellowship, and NSF Award CCF- 0953960 (CAREER) and CCF-1101491.
Ilias Diakonikolas
University of Edinburgh
ilias.d@ed.ac.uk
Research supported by a Simons Foundation Postdoctoral Fellowship. Some of this work was done while at Columbia University, supported by NSF grant CCF-0728736, and by an Alexander S. Onassis Foundation Fellowship.
Rocco A. Servedio
Columbia University
rocco@cs.columbia.edu
Supported by NSF grants CNS-0716245, CCF-1115703, and CCF-1319788 and by DARPA award HR0011-08-1-0069.
###### Abstract

We consider a basic problem in unsupervised learning: learning an unknown Poisson Binomial Distribution. A Poisson Binomial Distribution (PBD) over is the distribution of a sum of independent Bernoulli random variables which may have arbitrary, potentially non-equal, expectations. These distributions were first studied by S. Poisson in 1837 [Poi37] and are a natural -parameter generalization of the familiar Binomial Distribution. Surprisingly, prior to our work this basic learning problem was poorly understood, and known results for it were far from optimal.

We essentially settle the complexity of the learning problem for this basic class of distributions. As our first main result we give a highly efficient algorithm which learns to -accuracy (with respect to the total variation distance) using samples independent of . The running time of the algorithm is quasilinear in the size of its input data, i.e., bit-operations.111We write to hide factors which are polylogarithmic in the argument to ; thus, for example, denotes a quantity which is for some absolute constant . (Observe that each draw from the distribution is a -bit string.) Our second main result is a proper learning algorithm that learns to -accuracy using samples, and runs in time . This sample complexity is nearly optimal, since any algorithm for this problem must use samples. We also give positive and negative results for some extensions of this learning problem to weighted sums of independent Bernoulli random variables.

## 1 Introduction

We begin by considering a somewhat fanciful scenario: You are the manager of an independent weekly newspaper in a city of people. Each week the -th inhabitant of the city independently picks up a copy of your paper with probability . Of course you do not know the values ; each week you only see the total number of papers that have been picked up. For many reasons (advertising, production, revenue analysis, etc.) you would like to have a detailed “snapshot” of the probability distribution (pdf) describing how many readers you have each week. Is there an efficient algorithm to construct a high-accuracy approximation of the pdf from a number of observations that is independent of the population ? We show that the answer is “yes.”

A Poisson Binomial Distribution of order is the distribution of a sum

 X=n∑i=1Xi,

where are independent Bernoulli (0/1) random variables. The expectations need not all be the same, and thus these distributions generalize the Binomial distribution and, indeed, comprise a much richer class of distributions. (See Section 1.2 below.) It is believed that Poisson [Poi37] was the first to consider this extension of the Binomial distribution222We thank Yuval Peres and Sam Watson for this information [PW11]. and the distribution is sometimes referred to as “Poisson’s Binomial Distribution” in his honor; we shall simply call these distributions PBDs.

PBDs are one of the most basic classes of discrete distributions; indeed, they are arguably the simplest -parameter probability distribution that has some nontrivial structure. As such they have been intensely studied in probability and statistics (see Section 1.2) and arise in many settings; for example, we note here that tail bounds on PBDs form an important special case of Chernoff/Hoeffding bounds [Che52, Hoe63, DP09]. In application domains, PBDs have many uses in research areas such as survey sampling, case-control studies, and survival analysis, see e.g., [CL97] for a survey of the many uses of these distributions in applications. Given the simplicity and ubiquity of these distributions, it is quite surprising that the problem of density estimation for PBDs (i.e., learning an unknown PBD from independent samples) is not well understood in the statistics or learning theory literature. This is the problem we consider, and essentially settle, in this paper.

We work in a natural PAC-style model of learning an unknown discrete probability distribution which is essentially the model of [KMR94]. In this learning framework for our problem, the learner is provided with the value of and with independent samples drawn from an unknown PBD . Using these samples, the learner must with probability at least output a hypothesis distribution such that the total variation distance is at most , where are accuracy and confidence parameters that are provided to the learner.333[KMR94] used the Kullback-Leibler divergence as their distance measure but we find it more natural to use variation distance. A proper learning algorithm in this framework outputs a distribution that is itself a Poisson Binomial Distribution, i.e., a vector which describes the hypothesis PBD where .

### 1.1 Our results.

Our main result is an efficient algorithm for learning PBDs from many samples independent of Since PBDs are an -parameter family of distributions over the domain , we view such a tight bound as a surprising result. We prove:

###### Theorem 1 (Main Theorem).

Let be an unknown PBD.

1. [Learning PBDs from constantly many samples] There is an algorithm with the following properties: given and access to independent draws from , the algorithm uses

 ~O((1/ϵ3)⋅log(1/δ))

samples from , performs

 ~O((1/ϵ3)⋅logn⋅log21δ)

bit operations, and with probability at least outputs a (succinct description of a) distribution over which is such that

2. [Properly learning PBDs from constantly many samples] There is an algorithm with the following properties: given and access to independent draws from , the algorithm uses

 ~O(1/ϵ2)⋅log(1/δ)

samples from , performs

 (1/ϵ)O(log2(1/ϵ))⋅~O(logn⋅log1δ)

bit operations, and with probability at least outputs a (succinct description of a) vector defining a PBD such that

We note that, since every sample drawn from is a -bit string, for constant the number of bit-operations performed by our first algorithm is quasilinear in the length of its input. Moreover, the sample complexity of both algorithms is close to optimal, since samples are required even to distinguish the (simpler) Binomial distributions and , which have total variation distance Indeed, in view of this observation, our second algorithm is essentially sample-optimal.

Motivated by these strong learning results for PBDs, we also consider learning a more general class of distributions, namely distributions of the form which are weighted sums of independent Bernoulli random variables. We give an algorithm which uses samples and runs in time if there are only constantly many different weights in the sum:

###### Theorem 2 (Learning sums of weighted independent Bernoulli random variables).

Let be a weighted sum of unknown independent Bernoullis such that there are at most different values among Then there is an algorithm with the following properties: given and access to independent draws from , it uses

 ˜O(k/ϵ2)⋅log(n)⋅log(1/δ)

samples from , runs in time

 poly(nk⋅ϵ−klog2(1/ϵ))⋅log(1/δ),

and with probability at least outputs a hypothesis vector defining independent Bernoulli random variables with such that where .

To complement Theorem 2, we also show that if there are many distinct weights in the sum, then even for weights with a very simple structure any learning algorithm must use many samples:

###### Theorem 3 (Sample complexity lower bound for learning sums of weighted independent Bernoullis).

Let be a weighted sum of unknown independent Bernoullis (where the -th weight is simply ). Let be any learning algorithm which, given and access to independent draws from , outputs a hypothesis distribution such that with probability at least Then must use samples.

### 1.2 Related work.

At a high level, there has been a recent surge of interest in the theoretical computer science community on fundamental algorithmic problems involving basic types of probability distributions, see e.g., [KMV10, MV10, BS10, VV11] and other recent papers; our work may be considered as an extension of this theme. More specifically, there is a broad literature in probability theory studying various properties of PBDs; see [Wan93] for an accessible introduction to some of this work. In particular, many results study approximations to the Poisson Binomial distribution via simpler distributions. In a well-known result, Le Cam [Cam60] shows that for any PBD with , it holds that

 dTV(X,Poi(n∑i=1pi))≤2n∑i=1p2i,

where is the Poisson distribution with parameter . Subsequently many other proofs of this result and similar ones were given using a range of different techniques; [HC60, Che74, DP86, BHJ92] is a sampling of work along these lines, and Steele [Ste94] gives an extensive list of relevant references. Much work has also been done on approximating PBDs by normal distributions (see e.g., [Ber41, Ess42, Mik93, Vol95]) and by Binomial distributions (see e.g., [Ehm91, Soo96, Roo00]). These results provide structural information about PBDs that can be well-approximated via simpler distributions, but fall short of our goal of obtaining approximations of an unknown PBD up to arbitrary accuracy. Indeed, the approximations obtained in the probability literature (such as the Poisson, Normal and Binomial approximations) typically depend only on the first few moments of the target PBD. This is attractive from a learning perspective because it is possible to efficiently estimate such moments from random samples, but higher moments are crucial for arbitrary approximation [Roo00].

Taking a different perspective, it is easy to show (see Section 2 of [KG71]) that every PBD is a unimodal distribution over . (Recall that a distribution over is unimodal if there is a value such that for and for ) The learnability of general unimodal distributions over is well understood: Birgé [Bir87a, Bir97] has given a computationally efficient algorithm that can learn any unimodal distribution over to variation distance from samples, and has shown that any algorithm must use samples. (The [Bir87a, Bir97] upper and lower bounds are stated for continuous unimodal distributions, but the arguments are easily adapted to the discrete case.) Our main result, Theorem 1, shows that the additional PBD assumption can be leveraged to obtain sample complexity independent of with a computationally highly efficient algorithm.

So, how might one leverage the structure of PBDs to remove from the sample complexity? A first observation is that a PBD assigns of its mass to points. So one could draw samples to (approximately) identify these points and then try to estimate the probability assigned to each such point, but clearly such an approach, if followed naïvely, would give sample complexity. Alternatively, one could run Birgé’s algorithm on the restricted support of size , but that will not improve the asymptotic sample complexity. A different approach would be to construct a small -cover (under the total variation distance) of the space of all PBDs on variables. Indeed, if such a cover has size , it can be shown (see Lemma 10 in Section 3.1, or Chapter 7 of [DL01])) that a target PBD can be learned from samples. Still it is easy to argue that any cover needs to have size , so this approach too gives a dependence in the sample complexity.

Our approach, which removes completely from the sample complexity, requires a refined understanding of the structure of the set of all PBDs on variables, in fact one that is more refined than the understanding provided by the aforementioned results (approximating a PBD by a Poisson, Normal, or Binomial distribution). We give an outline of the approach in the next section.

### 1.3 Our approach.

The starting point of our algorithm for learning PBDs is a theorem of [DP11, Das08] that gives detailed information about the structure of a small -cover (under the total variation distance) of the space of all PBDs on variables (see Theorem 4). Roughly speaking, this result says that every PBD is either close to a PBD whose support is sparse, or is close to a translated “heavy” Binomial distribution. Our learning algorithm exploits this structure of the cover; it has two subroutines corresponding to these two different types of distributions that the cover contains. First, assuming that the target PBD is close to a sparsely supported distribution, it runs Birgé’s unimodal distribution learner over a carefully selected subinterval of to construct a hypothesis ; the (purported) sparsity of the distribution makes it possible for this algorithm to use samples independent of . Then, assuming that the target PBD is close to a translated “heavy” Binomial distribution, the algorithm constructs a hypothesis Translated Poisson Distribution [R0̈7] whose mean and variance match the estimated mean and variance of the target PBD; we show that is close to the target PBD if the target PBD is not close to any sparse distribution in the cover. At this point the algorithm has two hypothesis distributions, and , one of which should be good; it remains to select one as the final output hypothesis. This is achieved using a form of “hypothesis testing” for probability distributions.

The above sketch captures the main ingredients of Part (1) of Theorem 1, but additional work needs to be done to get the proper learning algorithm of Part (2). For the non-sparse case, first note that the Translated Poisson hypothesis is not a PBD. Via a sequence of transformations we are able to show that the Translated Poisson hypothesis can be converted to a Binomial distribution for some To handle the sparse case, we use an alternate learning approach: instead of using Birgé’s unimodal algorithm (which would incur a sample complexity of ), we first show that, in this case, there exists an efficiently constructible -cover of size and then apply a general learning result that we now describe.

The general learning result that we use (Lemma 10) is the following: We show that for any class of target distributions, if has an -cover of size then there is a generic algorithm for learning an unknown distribution from to accuracy that uses samples. Our approach is rather similar to the algorithm of [DL01] for choosing a density estimate (but different in some details); it works by carrying out a tournament that matches every pair of distributions in the cover against each other. Our analysis shows that with high probability some -accurate distribution in the cover will survive the tournament undefeated, and that any undefeated distribution will with high probability be -accurate.

Applying this general result to the -cover of size described above, we obtain a PBD that is -close to the target (this accounts for the increased running time in Part (2) versus Part (1)). We stress that for both the non-proper and proper learning algorithms sketched above, many technical subtleties and challenges arise in implementing the high-level plan given above, requiring a careful and detailed analysis.

We prove Theorem 2 using the general approach of Lemma 10 specialized to weighted sums of independent Bernoullis with constantly many distinct weights. We show how the tournament can be implemented efficiently for the class of weighted sums of independent Bernoullis with constantly many distinct weights, and thus obtain Theorem 2. Finally, the lower bound of Theorem 3 is proved by a direct information-theoretic argument.

### 1.4 Preliminaries.

#### Distributions.

For a distribution supported on we write to denote the value of the probability density function (pdf) at point , and to denote the value of the cumulative density function (cdf) at point . For , we write to denote and to denote the conditional distribution of restricted to Sometimes we write and for a subset , meaning and respectively.

#### Total Variation Distance.

Recall that the total variation distance between two distributions and over a finite domain is

 dTV(X,Y) := (1/2)⋅∑α∈D|X(α)−Y(α)|=maxS⊆D[X(S)−Y(S)].

Similarly, if and are two random variables ranging over a finite set, their total variation distance is defined as the total variation distance between their distributions. For convenience, we will often blur the distinction between a random variable and its distribution.

#### Covers.

Fix a finite domain , and let denote some set of distributions over Given , a subset is said to be a -cover of (w.r.t. the total variation distance) if for every distribution in there exists some distribution in such that We sometimes say that distributions are -neighbors if . If this holds, we also say that is -close to and vice versa.

#### Poisson Binomial Distribution.

A Poisson binomial distribution of order is a sum of mutually independent Bernoulli random variables . We denote the set of all Poisson binomial distributions of order by and, if is clear from context, just .

A Poisson binomial distribution can be represented uniquely as a vector satisfying . To go from to its corresponding vector, we find a collection of mutually independent Bernoullis such that is distributed according to and . (Such a collection exists by the definition of a Poisson binomial distribution.) Then we set for all . Lemma 1 of [DP13] shows that the resulting vector is unique.

We denote by the distribution of the sum of mutually independent indicators with expectations , for all . Given the above discussion is unique up to permutation of the ’s. We also sometimes write to denote the distribution of Note the difference between , which refers to the distribution of , and , which refers to the underlying collection of mutually independent Bernoulli random variables.

#### Translated Poisson Distribution.

We will make use of the translated Poisson distribution for approximating the Poisson Binomial distribution. We define the translated Poisson distribution, and state a known result on how well it approximates the Poisson Binomial distribution.

###### Definition 1 ([R0̈7]).

We say that an integer random variable is distributed according to the translated Poisson distribution with parameters and , denoted , iff can be written as

 Y=⌊μ−σ2⌋+Z,

where is a random variable distributed according to , where represents the fractional part of .

The following lemma gives a useful bound on the variation distance between a Poisson Binomial Distribution and a suitable translated Poisson distribution. Note that if the variance of the Poisson Binomial Distribution is large, then the lemma gives a strong bound.

###### Lemma 1 (see (3.4) of [R0̈7]).

Let be independent random indicators with . Then

 dTV(n∑i=1Ji,TP(μ,σ2))≤√∑ni=1p3i(1−pi)+2∑ni=1pi(1−pi),

where and .

The following bound on the total variation distance between translated Poisson distributions will be useful.

###### Lemma 2 (Lemma 2.1 of [Bl06]).

For and with , we have

 dTV(TP(μ1,σ21),TP(μ2,σ22))≤|μ1−μ2|σ1+|σ21−σ22|+1σ21.

#### Running Times, and Bit Complexity.

Throughout this paper, we measure the running times of our algorithms in numbers of bit operations. For a positive integer , we denote by its description complexity in binary, namely . Moreover, we represent a positive rational number as , where and are relatively prime positive integers. The description complexity of is defined to be . We will assume that all ’s and ’s input to our algorithms are rational numbers.

## 2 Learning a sum of Bernoulli random variables from poly(1/ϵ) samples

In this section, we prove Theorem 1 by providing a sample- and time-efficient algorithm for learning an unknown PBD . We start with an important ingredient in our analysis.

A cover for PBDs. We make use of the following theorem, which provides a cover of the set of all PBDs of order-. The theorem was given implicitly in [DP11] and explicitly as Theorem 1 in [DP13].

###### Theorem 4 (Cover for PBDs).

For all , there exists an -cover of such that

1. ; and

2. can be constructed in time linear in its representation size, i.e., .

Moreover, if , then the collection of Bernoulli random variables has one of the following forms, where is a positive integer, for some absolute constant :

• (-Sparse Form) There is some such that, for all , and, for all , .

• (-heavy Binomial Form) There is some and such that, for all , and, for all , ; moreover, satisfy and

Finally, for every for which there is no -neighbor in that is in sparse form, there exists some in -heavy Binomial form such that

• ; and

• if , , and , then and

We remark that the cover theorem as stated in [DP13] does not include the part of the above statement following “finally.” We provide a proof of this extension in Appendix A.

The Basic Learning Algorithm. The high-level structure of our learning algorithms which give Theorem 1 is provided in Algorithm Learn-PBD of Figure 1. We instantiate this high-level structure, with appropriate technical modifications, in Section 2.4, where we give more detailed descriptions of the non-proper and proper algorithms that give parts (1) and (2) of Theorem 1.

At a high level, the subroutine Learn-Sparse is given sample access to and is designed to find an -accurate hypothesis with probability at least , if the unknown PBD is -close to some sparse form PBD inside the cover . Similarly, Learn-Poisson is designed to find an -accurate hypothesis , if is not -close to a sparse form PBD (in this case, Theorem 4 implies that must be -close to some -heavy Binomial form PBD). Finally, Choose-Hypothesis is designed to choose one of the two hypotheses as being -close to The following subsections specify these subroutines, as well as how the algorithm can be used to establish Theorem 1. We note that Learn-Sparse and Learn-Poisson do not return the distributions and as a list of probabilities for every point in . They return instead a succinct description of these distributions in order to keep the running time of the algorithm logarithmic in . Similarly, Choose-Hypothesis operates with succinct descriptions of these distributions.

### 2.1 Learning when X is close to a sparse form PBD.

Our starting point here is the simple observation that any PBD is a unimodal distribution over the domain . (There is a simple inductive proof of this, or see Section 2 of [KG71].) This enables us to use the algorithm of Birgé [Bir97] for learning unimodal distributions. We recall Birgé’s result, and refer the reader to Appendix B for an explanation of how Theorem 5 as stated below follows from [Bir97].

###### Theorem 5 ([Bir97]).

For all , there is an algorithm that draws

 O(lognϵ3log1δ+1ϵ2log1δloglog1δ)

samples from an unknown unimodal distribution over , does

 ~O(log2nϵ3log21δ)

bit-operations, and outputs a (succinct description of a) hypothesis distribution over that has the following form: is uniform over subintervals , whose union , where In particular, the algorithm outputs the lists through and through , as well as the total probability mass that assigns to each subinterval , . Finally, with probability at least , .

The main result of this subsection is the following:

###### Lemma 3.

For all , there is an algorithm Learn-Sparse that draws

 O(1ϵ′3log1ϵ′log1δ′+1ϵ′2log1δ′loglog1δ′)

samples from a target PBD over , does

 logn⋅~O(1ϵ′3log21δ′)

bit operations, and outputs a (succinct description of a) hypothesis distribution over that has the following form: its support is contained in an explicitly specified interval , where , and for every point in the algorithm explicitly specifies the probability assigned to that point by 444In particular, our algorithm will output a list of pointers, mapping every point in to some memory location where the probability assigned to that point by is written. The algorithm has the following guarantee: if is -close to some sparse form PBD in the cover of Theorem 4, then with probability at least , , for some absolute constant , and the support of lies in the support of .

The high-level idea of Lemma 3 is quite simple. We truncate of the probability mass from each end of to obtain a conditional distribution ; since is unimodal so is . If is larger than then the algorithm outputs “fail” (and could not have been close to a sparse-form distribution in the cover). Otherwise, we use Birgé’s algorithm to learn the unimodal distribution . A detailed description of the algorithm is given in Figure 2 below.

Proof of Lemma 3: As described in Figure 2, algorithm Learn-Sparse first draws samples from and sorts them to obtain a list of values We claim the following about the values and defined in Step 2 of the algorithm:

###### Claim 4.

With probability at least , we have and .

###### Proof.

We only show that with probability at least , since the arguments for , and are identical. Given that each of these conditions is met with probability at least , the union bound establishes our claim.

To show that is satisfied with probability at least we argue as follows: Let . Clearly, while . Given this, if samples are drawn from then the expected number of them that are is at most . It follows then from the Chernoff bound that the probability that more than samples are is at most . Hence except with this failure probability, we have , which implies that . ∎

As specified in Steps 3 and 4, if , where is the constant in the statement of Theorem 4, the algorithm outputs “fail”, returning the trivial hypothesis which puts probability mass on the point . Otherwise, the algorithm runs Birgé’s unimodal distribution learner (Theorem 5) on the conditional distribution , and outputs the result of Birgé’s algorithm. Since is unimodal, it follows that is also unimodal, hence Birgé’s algorithm is appropriate for learning it. The way we apply Birgé’s algorithm to learn given samples from the original distribution is the obvious one: we draw samples from , ignoring all samples that fall outside of , until the right number of samples fall inside , as required by Birgé’s algorithm for learning a distribution of support of size with probability at least . Once we have the right number of samples in , we run Birgé’s algorithm to learn the conditional distribution . Note that the number of samples we need to draw from until the right number of samples fall inside is still , with probability at least . Indeed, since , it follows from the Chernoff bound that with probability at least , if samples are drawn from , at least fall inside .

Analysis: It is easy to see that the sample complexity of our algorithm is as promised. For the running time, notice that, if Birgé’s algorithm is invoked, it will return two lists of numbers through and through , as well as a list of probability masses assigned to each subinterval , , by the hypothesis distribution , where . In linear time, we can compute a list of probabilities , representing the probability assigned by to every point of subinterval , for . So we can represent our output hypothesis via a data structure that maintains pointers, having one pointer per point inside . The pointers map points to probabilities assigned by to these points. Thus turning the output of Birgé’s algorithm into an explicit distribution over incurs linear overhead in our running time, and hence the running time of our algorithm is also as promised. (See Appendix B for an explanation of the running time of Birgé’s algorithm.) Moreover, we also note that the output distribution has the promised structure, since in one case it has a single atom at and in the other case it is the output of Birgé’s algorithm on a distribution of support of size .

It only remains to justify the last part of the lemma. Let be the sparse-form PBD that is close to; say that is supported on where Since is -close to in total variation distance it must be the case that . Since by Claim 4, it must be the case that . Similar arguments give that . So the interval is contained in and has length at most . This means that Birgé’s algorithm is indeed used correctly by our algorithm to learn , with probability at least (that is, unless Claim 4 fails). Now it follows from the correctness of Birgé’s algorithm (Theorem 5) and the discussion above, that the hypothesis output when Birgé’s algorithm is invoked satisfies with probability at least , i.e., unless either Birgé’s algorithm fails, or we fail to get the right number of samples landing inside . To conclude the proof of the lemma we note that:

 2dTV(X,X[^a,^b]) = ∑i∈[^a,^b]|X[^a,^b](i)−X(i)|+∑i∉[^a,^b]|X[^a,^b](i)−X(i)| = ∑i∈[^a,^b]∣∣1X([^a,^b])X(i)−X(i)∣∣+∑i∉[^a,^b]X(i) = ∑i∈[^a,^b]∣∣11−O(ϵ′)X(i)−X(i)∣∣+O(ϵ′) = O(ϵ′)1−O(ϵ′)∑i∈[^a,^b]∣∣X(i)∣∣+O(ϵ′) = O(ϵ′).

So the triangle inequality gives: , and Lemma 3 is proved.

### 2.2 Learning when X is close to a k-heavy Binomial Form PBD.

###### Lemma 5.

For all , there is an algorithm Learn-Poisson that draws

 O(log(1/δ′)/ϵ′2)

samples from a target PBD over , does

 O(logn⋅log(1/δ′)/ϵ′2)

bit operations, and returns two parameters and . The algorithm has the following guarantee: Suppose is not -close to any sparse form PBD in the cover of Theorem 4. Let be the translated Poisson distribution with parameters and . Then with probability at least we have for some absolute constant .

Our proof plan is to exploit the structure of the cover of Theorem 4. In particular, if is not -close to any sparse form PBD in the cover, it must be -close to a PBD in heavy Binomial form with approximately the same mean and variance as , as specified by the final part of the cover theorem. Hence, a natural strategy is to obtain estimates and of the mean and variance of the unknown PBD , and output as a hypothesis a translated Poisson distribution with parameters and . We show that this strategy is a successful one. Before providing the details, we highlight two facts that we will establish in the subsequent analysis and that will be used later. The first is that, assuming is not -close to any sparse form PBD in the cover , its variance satisfies

 σ2=Ω(1/ϵ′2)≥θ2for some % universal constant~{}θ. (1)

The second is that under the same assumption, the estimates and of the mean and variance of that we obtain satisfy the following bounds with probability at least :

 |μ−^μ|≤ϵ′⋅σand|σ2−^σ2|≤ϵ′⋅σ2. (2)

See Figure 3 and the associated Figure 4 for a detailed description of the Learn-Poisson algorithm.

Proof of Lemma 5: We start by showing that we can estimate the mean and variance of the target PBD .

###### Lemma 6.

For all , there exists an algorithm with the following properties: given access to a PBD of order , it produces estimates and for and respectively such that with probability at least :

 |μ−^μ|≤ϵ⋅σ~{}~{}~{}~{}~{}~{}~{} and ~{}~{}~{% }~{}~{}~{}~{}|σ2−^σ2|≤ϵ⋅σ2√4+1σ2.

The algorithm uses

 O(log(1/δ)/ϵ2)

samples and runs in time

 O(lognlog(1/δ)/ϵ2).
###### Proof.

We treat the estimation of and separately. For both estimation problems we show how to use samples to obtain estimates and achieving the required guarantees with probability at least (we refer to these as “weak estimators”). Then a routine procedure allows us to boost the success probability to at the expense of a multiplicative factor on the number of samples. While we omit the details of the routine boosting argument, we remind the reader that it involves running the weak estimator times to obtain estimates and outputting the median of these estimates, and similarly for estimating .

We proceed to specify and analyze the weak estimators for and separately:

• Weak estimator for : Let be independent samples from , and let . Then

 E[^μ]=μ  and  Var[^μ]=1mVar[X]=1mσ2.

So Chebyshev’s inequality implies that

 Pr[|^μ−μ|≥tσ/√m]≤1t2.

Choosing


text:
# State-dependent diversification with the ClaSSE model

## Introduction

In the previous examples we have modeled all character state transitions as anagenetic changes. Anagenetic changes occur along the branches of a phylogeny, within a lineage. Cladogenetic changes, on the other hand, occur at speciation events. They represent changes in a character state that may be associated with speciation events due to increased reproductive isolation, for example colonizing a new geographic area or a shift in chromosome number. Note that it can be quite tricky to determine if a character state shift is a cause or a consequence of speciation, but we can at least test if state changes tend to occur in the same time window as speciation events.

A major challenge for all phylogenetic models of cladogenetic character change is accounting for unobserved speciation events due to lineages going extinct and not leaving any extant descendants (Bokma 2002), or due to incomplete sampling of lineages in the present. Teasing apart the phylogenetic signal for cladogenetic and anagenetic processes given unobserved speciation events is a major difficulty. Commonly used biographic models like the dispersal-extinction-cladogenesis (DEC; Ree and Smith (2008)) simply ignore unobserved speciation events and so result in biased estimates of cladogenetic versus anagenetic change.

This bias can be avoided by using the Cladogenetic State change Speciation and Extinction (ClaSSE) model (Goldberg and Igić 2012), which accounts for unobserved speciation events by jointly modeling both character evolution and the phylogenetic birth-death process. ClaSSE models extend other SSE models by incorporating both cladogenetic and anagenetic character evolution. This approach has been used to model biogeographic range evolution (Goldberg et al. 2011) and chromosome number evolution (missing reference).

Here we will use RevBayes to examine biogeographic range evolution in the primates. We will model biogeographic range evolution similar to a DEC model, however we will use ClaSSE to account for speciation events unobserved due to extinction or incomplete sampling.

## Setting up the analysis

### Reading in the data

Begin by reading in the observed tree.

observed_phylogeny <- readTrees("data/primates_biogeo.tre")[1]

Get the taxa in the tree. We’ll need this later on.

taxa = observed_phylogeny.taxa()

Now let’s read in the biogeographic range data. The areas are represented as the following character states:

• 0 = 00 = the null state with no range
• 1 = 01 = New World only
• 2 = 10 = Old World only
• 3 = 11 = both New and Old World

For consistency, we have chosen to use the same representation of biogeographic ranges used in the \RevBayes biogeography/DEC tutorial. Each range is represented as both a natural number (0, 1, 2, 3) and a corresponding bitset (00, 01, 10, 11). The null state (state 0) is used in DEC models to represent a lineage that has no biogeographic range and is therefore extinct. Our model will include this null state as well, however, we will explicitly model extinction as part of the birth-death process so our character will never enter state 0.

data_biogeo = readCharacterDataDelimited("data/primates_biogeo.tsv", stateLabels="0123", type="NaturalNumbers", delimiter="\t", headers=TRUE)

Also we need to create the move and monitor vectors.

moves = VectorMoves()
monitors = VectorMonitors()

### Set up the extinction rates

We are going to draw both anagenetic transition rates and diversification rates from a lognormal distribution. The mean of the prior distribution will be $\ln(\frac{\text{#Taxa}}{2}) / \text{tree-age}$ which is the expected net diversification rate, and the SD will be 1.0 so the 95\% prior interval ranges well over 2 orders of magnitude.

num_species <- 424 # approximate total number of primate species
rate_mean <- ln( ln(num_species/2.0) / observed_phylogeny.rootAge() )
rate_sd <- 1.0

The extinction rates will be stored in a vector where each element represents the extinction rate for the corresponding character state. We have chosen to allow a lineage to go extinct in both the New and Old World at the same time (like a global extinction event). As an alternative, you could restrict the model so that a lineage can only go extinct if it’s range is limited to one area.

extinction_rates[1] <- 0.0 # the null state (state 0)
extinction_rates[2] ~ dnLognormal(rate_mean, rate_sd) # extinction when the lineage is in New World (state 1)
extinction_rates[3] ~ dnLognormal(rate_mean, rate_sd) # extinction when the lineage is in Old World (state 2)
extinction_rates[4] ~ dnLognormal(rate_mean, rate_sd) # extinction when in both (state 3)

Note \Rev vectors are indexed starting with 1, yet our character states start at 0. So \texttt{extinction_rate[1]} will represent the extinction rate for character state 0.

Add MCMC moves for each extinction rate.

moves.append( mvSlide( extinction_rates[2], weight=4 ) )
moves.append( mvSlide( extinction_rates[3], weight=4 ) )
moves.append( mvSlide( extinction_rates[4], weight=4 ) )

Let’s also create a deterministic variable to monitor the overall extinction rate.

total_extinction := sum(extinction_rates)

### Set up the anagenetic transition rate matrix

First, let’s create the rates of anagenetic dispersal:

anagenetic_dispersal_13 ~ dnLognormal(rate_mean, rate_sd) # disperse from New to Old World 01 -> 11
anagenetic_dispersal_23 ~ dnLognormal(rate_mean, rate_sd) # disperse from Old to New World 10 -> 11

Now add MCMC moves for each anagenetic dispersal rate.

moves.append( mvSlide( anagenetic_dispersal_13, weight=4 ) )
moves.append( mvSlide( anagenetic_dispersal_23, weight=4 ) )

The anagenetic transitions will be stored in a 4 by 4 instantaneous rate matrix. We will construct this by first creating a vector of vectors. Let’s begin by initalizing all rates to 0.0:

for (i in 1:4) {
for (j in 1:4) {
r[i][j] <- 0.0
}
}

Now we can populate non-zero rates into the anagenetic transition rate matrix:

r[2][4] := anagenetic_dispersal_13
r[3][4] := anagenetic_dispersal_23
r[4][2] := extinction_rates[3]
r[4][3] := extinction_rates[2]

Note that we have modeled the rate of 11 $\rightarrow$ 01 (3 $\rightarrow$ 1) as being the rate of going extinct in area 2, and the rate of 11 $\rightarrow$ 10 (3 $\rightarrow$ 2) as being the rate of going extinct in area 1.

Now we pass our vector of vectors into the \cl{fnFreeK} function to create the instaneous rate matrix.

ana_rate_matrix := fnFreeK(r, rescaled=false)

### Set up the cladogenetic speciation rate matrix

Here we need to define each cladogenetic event type in the form [ancestor\_state, daughter1\_state, daughter2\_state] and assign each cladogenetic event type a corresponding speciation rate.

The first type of cladogenetic event we’ll specify is widespread sympatry. Widespread sympatric cladogenesis is where the biogeographic range does not change; that is the daughter lineages inherit the same range as the ancestor. In this example we are not going to allow the speciation events like 11 $\rightarrow$ 11, 11, as it seems biologically implausible. However if you wanted you could add this to your model.

Define the speciation rate for widespread sympatric cladogenesis events:

speciation_wide_sympatry ~ dnLognormal(rate_mean, rate_sd)
moves.append( mvSlide( speciation_wide_sympatry, weight=4 ) )

clado_events[1] = [1, 1, 1] # 01 -> 01, 01
clado_events[2] = [2, 2, 2] # 10 -> 10, 10

and assign each the same speciation rate:

speciation_rates[1] := speciation_wide_sympatry/2
speciation_rates[2] := speciation_wide_sympatry/2

Subset sympatry is where one daughter lineage inherits the full ancestral range but the other lineage inherits only a single region.

speciation_sub_sympatry ~ dnLognormal(rate_mean, rate_sd)
moves.append( mvSlide( speciation_sub_sympatry, weight=4 ) )

Define the subset sympatry events and assign each a speciation rate:

clado_events[3] = [3, 3, 1] # 11 -> 11, 01
clado_events[4] = [3, 1, 3] # 11 -> 01, 11
clado_events[5] = [3, 3, 2] # 11 -> 11, 10
clado_events[6] = [3, 2, 3] # 11 -> 10, 11
speciation_rates[3] := speciation_sub_sympatry/4
speciation_rates[4] := speciation_sub_sympatry/4
speciation_rates[5] := speciation_sub_sympatry/4
speciation_rates[6] := speciation_sub_sympatry/4

Allopatric cladogenesis is when the two daughter lineages split the ancestral range:

speciation_allopatry ~ dnLognormal(rate_mean, rate_sd)
moves.append( mvSlide( speciation_allopatry, weight=4 ) )

Define the allopatric events:

clado_events[7] = [3, 1, 2] # 11 -> 01, 10
clado_events[8] = [3, 2, 1] # 11 -> 10, 01
speciation_rates[7] := speciation_allopatry/2
speciation_rates[8] := speciation_allopatry/2

Now let’s create a deterministic variable to monitor the overall speciation rate:

total_speciation := sum(speciation_rates)

Finally, we construct the cladogenetic speciation rate matrix from the cladogenetic event types and the speciation rates.

clado_matrix := fnCladogeneticSpeciationRateMatrix(clado_events, speciation_rates, 4)

Let’s view the cladogenetic matrix to see if we have set it up correctly:

clado_matrix

### Set up the cladogenetic character state-dependent birth-death process

For simplicity we will fix the root frequencies to be equal except for the null state which has probability of 0.

root_frequencies <- simplex([0, 1, 1, 1])

rho is the probability of sampling species at the present:

rho <- observed_phylogeny.ntips()/num_species

Now we construct a stochastic variable drawn from the cladogenetic character state-dependent birth-death process.

classe ~ dnCDCladoBDP( rootAge         = observed_phylogeny.rootAge(),
extinctionRates = extinction_rates,
Q               = ana_rate_matrix,
delta           = 1.0,
pi              = root_frequencies,
rho             = rho,
condition       = "time",
taxa            = taxa )

Clamp the model with the observed data.

classe.clamp( observed_phylogeny )
classe.clampCharData( data_biogeo )

### Finalize the model

Just like before, we must create a workspace model object.

mymodel = model(classe)

\subsection{Set up and run the MCMC}

First, set up the monitors that will output parameter values to file and screen.

monitors.append( mnModel(filename="output/primates_ClaSSE.log", printgen=1) )
monitors.append( mnJointConditionalAncestralState(tree=observed_phylogeny, cdbdp=classe, type="NaturalNumbers", printgen=1, withTips=true, withStartStates=true, filename="output/anc_states_primates_ClaSSE.log") )
monitors.append( mnScreen(printgen=1, speciation_wide_sympatry, speciation_sub_sympatry, speciation_allopatry, extinction_rates) )

Now define our workspace MCMC object.

mymcmc = mcmc(mymodel, monitors, moves)

We will perform a pre-burnin to tune the proposals and then run the MCMC. Note that for a real analysis you would want to run the MCMC for many more iterations.

mymcmc.burnin(generations=200,tuningInterval=5)
mymcmc.run(generations=1000)

## Summarize ancestral states

When the analysis has completed you now summarize the ancestral states. The ancestral states are estimated both for the “beginning” and “end” state of each branch, so that the cladogenetic changes that occurred at speciation events are distinguished from the changes that occurred anagenetically along branches. Make sure the include_start_states argument is set to true.

anc_states = readAncestralStateTrace("output/anc_states_primates_ClaSSE.log")
anc_tree = ancestralStateTree(tree=observed_phylogeny, ancestral_state_trace_vector=anc_states, include_start_states=true, file="output/anc_states_primates_ClaSSE_results.tree", burnin=0, summary_statistic="MAP", site=0)

### Plotting ancestral states

Like before, we’ll plot the ancestral states using the RevGadgets R package. Execute the script plot_anc_states_ClaSSE.R in R. The results can be seen in Figure \ref{fig:results_ClaSSE}. The maximum a posteriori (MAP) estimate for each node is shown as well as the posterior probability of the states represented by the size of the dots.

library(RevGadgets)

tree_file = "output/anc_states_primates_ClaSSE_results.tree"

plot_ancestral_states(tree_file, summary_statistic="MAPRange",
tip_label_size=3,
tip_label_offset=1,
xlim_visible=c(0,100),
node_label_size=0,
shoulder_label_size=0,
include_start_states=TRUE,
show_posterior_legend=TRUE,
node_size_range=c(4, 7),
alpha=0.75)

output_file = "RevBayes_Anc_States_ClaSSE.pdf"
ggsave(output_file, width = 11, height = 9)

# Exercise

• Using either R or Tracer, visualize the posterior estimates for different types of cladogenetic events. What kind of speciation events are most common?
• As we have specified the model, we did not allow cladogenetic long distance (jump) dispersal, for example 01 $\rightarrow$ 01, 10. Modify this script to include cladogenetic long distance dispersal and calculate Bayes factors to see which model fits the data better. How does this affect the ancestral state estimate? }
1. Bokma F. 2002. Detection of punctuated equilibrium from molecular phylogenies. Journal of Evolutionary Biology. 15:1048–1056. 10.1046/j.1420-9101.2002.00458.x
2. Goldberg E.E., Lancaster L.T., Ree R.H. 2011. Phylogenetic Inference of Reciprocal Effects between Geographic Range Evolution and Diversification. Systematic Biology. 60:451–465. 10.1093/sysbio/syr046
3. Goldberg E.E., Igić B. 2012. Tempo and Mode in Plant Breeding System Evolution. Evolution. 66:3701–3709. 10.1111/j.1558-5646.2012.01730.x
4. Ree R.H., Smith S.A. 2008. Maximum Likelihood Inference of Geographic Range Evolution by Dispersal, Local Extinction, and Cladogenesis. Systematic Biology. 57:4–14. 10.1080/10635150701883881