#IceBucketChallenge - A Viral Campaign

I know this blog is a little late to the game. In fact, I started writing it while it was still relevant... but then I got distracted. Such is life.


Background

You probably already know (or vaguely remember) what the Ice Bucket Challenge is. Basically if you're nominated you have to dump a bucket of ice water over your head, or else you have to donate money to some charity - most commonly ALSA, the Amyotrophic Lateral Sclerosis Association.

You then nominate 3 more people, who have 24 hours to do the same thing. In some variants, you also make a small donation even if you dump the ice water over your head, or else you make a bigger donation if you don't (some specify $100).

Anyway, what I was interested in is how the challenge spread - this was a 'viral' campaign in a very true sense. So why not try to model how the campaign spread as we would model the spread of a virus?


The Viral Model

I've written about this sort of thing several times before. The basic idea is this - the population is divided into three groups:

Susceptible (S) - The population that hasn't been exposed to a disease/virus, but is susceptible to infection
Infectious (I) - The people who have been exposed, and can infect other people
Recovered (R) - The people who have been infected, but have recovered. It's generally assumed that these people can't be re-infected. (But there are variations).

For the ice bucket challenge, we can look at a direct analogy as

S - Those who haven't been nominated
I - Those who have been nominated and are taking the challenge (so can nominate other people)
R - Those who have completed or those who declined the challenge (and can't nominate anyone else)

We can draw the model as a flowchart, showing how people move between the different groups,

Here, we've divided nominees into those who accept the challenge (S->I) and those who don't accept (S->R). So, now we can describe the model as a system of differential equations,

Where the factors (α, β, γ) describe what proportions of each group move to where.

The most interesting part is the 'S*I' terms in the first two equations - this basically says that the number of newly infected people is proportional to the number of susceptible people AND the number of infectious people.

This makes the system self-limiting, meaning that the number of infected people can't just grow to infinity - since that's not what we observe. Instead, what we see is an increase to some maximum, followed by a steady drop-off. For example,

We don't have hard data on how many people did the challenge over time, but we can get an idea of what happened by looking at the YouTube search numbers for the phrase 'ice bucket challenge' (above, via Google Trends).

I mean, it seems reasonable to assume some loose correlation between the number of people doing the challenge, the number of videos of people doing the challenge, and the number of searches for those videos. Searches peaked on August 21st, in case you were wondering.


A More Discrete Model

The downside to this 'S*I' term is that it makes the equations non-linear, meaning they can't be solved analytically - that is, you can't come up with an 'exact' equation for the number of people doing the challenge on a given day, for example.

But we can do a numerical simulation, instead.

In fact, this is a more reasonable way of looking at the system, since we're interested in a discrete time step of one day - i.e. the 24 hours nominees have to complete the challenge. For the simulation, we're also going to rounded the numbers of people in each group (after each time step) to whole numbers, since you can't (or shouldn't) divide a person into fractions.

Now, we need to define some parameters. First of all we need to define our start populations. We'll call the initial susceptible population S0 - this could be, for example, the Earth's population (which was around 7.16bn when I ran my simulations). We'll assume that one person is infected as a starting point  - the originator of the challenge (I0 = 1). And we'll assume that no-one else has done the challenge at the start, so R0 = 0.

For the constants (α,β,γ), the easiest one to define is γ - the 'recovery' rate. We assume that after the 24 hour challenge period nominees are no longer 'infectious', therefore γ = 1.

For α and β, we have the total rate of infection/nomination defined as (α+β). We're given that each challenge completer gets to nominate three new people, so we can define (α+β) such that in the first step (when there's only one infectious person) we have (α+β)*S0*I0 = (α+β)*S0 = 3. Therefore (α+β) = 3/S0.

Now, α and β are related to the proportions of nominees that accept and decline the challenge, respectively. So we can redefine the constants as α = a/S0 and β = b/S0, such that (a+b) = 3, or alternatively α = a/S0 and β = (3-a)/S0. So now we can look at 'a' as the average number of nominees who accept the challenge.

So to tie it all together, we have the system of (difference) equations,

And it's pretty straightforward to write some code that'll run through those equations.


So we've got our model, what now?

Having a model is all well and good, but why bother? Well, now we can start asking questions. For example, how fast does the campaign spread? How long will it take for the challenge to die out, and how many people will have taken part by that point? And what happens when we change the number of people who accept the challenge?

First of all, if we run the simulation (with a = 2.5) and plot I(n) - the number of people doing the challenge on any given day - we get something like this

[where I(n) has been normalised so that the maximum is 100, as Google Trends does].

This is about what we expected - a rise and a fall. Though it's worth noting the shape is a little different from the YouTube searches above; the simulation drops off quickly, while the search numbers have a longer tail. This could be because, even after the challenges are done, there's a latent interest in (re)watching the videos. Or it could just be that this model is not entirely accurate...


So what happens when we vary 'a'?

We can start by assuming everyone accepts the challenge (a = 3). In this case, eventually everyone in the world does the challenge - specifically within 21 days of the first challenge. But that isn't very realistic.

So let's say that on average 1.5 or 2 of those challenged accept. In these cases, the ice bucket phenomenon ends before everyone can be challenged. For a = 2, the challenge ends after 50 day, with ~13% of the population going unchallenged. On the other hand, for a = 1.5 the challenge takes significantly longer to end (over a year).

In fact it turns out that for any a < 1.6030165.. the challenge will take a significant amount of time to end - the number of 'infectious' people will eventually reach 1, and stay there until S <= S0/(2a) (remember we're rounding each group to the nearest whole number). For a <= 0.5 the challenge ends almost immediately.

Now, we know that for a = 3, the population S goes to zero (everyone is challenged), whereas for a = 2 the challenge stops before S can go to zero. So we have the question - what is the smallest value of 'a' for which S goes to zero? If you do a bit of interpolating, you can figure out that this critical value comes out at around 2.7405349.. At this value of 'a', the challenge ends after 23 days. In other words, if on average 2.74.. (~91%) of the people nominated accept the challenge, then after 23 days everyone in the world will have been challenged.

From what I've seen myself, it seems like nearer 1 in 3 people accept the challenge on average. So, as viral as the campaign was, it was never going to take over the world.

If you play around a bit, you'll find that the critical values of 'a' are dependent on the initial population (S0). But, as far as I can tell, it's not possible to derive these critical values analytically. (Answers in the comments if you can prove otherwise).


Why wait to be nominated..?

At this point, you're probably thinking this isn't a very realistic model. And you'd be right. So lets make it a bit more complicated.

In particular, lets add spontaneous participation - that is, people who aren't directly nominated, but who see all the other people doing the challenge, and decide they want to take part too.

We'll assume that this participation is proportional to those who have already done the challenge (I and R). So to start with, we need to separate the 'recovered' group into those who actually did the challenge (R), and those who declined (D).

The new model looks something like this,

With difference equations,

In the flowchart, we've introduced this new factor 'δ', the 'inspiration' rate. If we re-define it, like we did for the infection rate, as d/S0, then 'd' can be loosely interpreted as the average number of people inspired to take part by each person who's already done the challenge.

Now we can look at how this 'd' factor affects the things we looked at before - how long it takes for the challenge to end, etc.

Let's start by assuming that the average number of nominees accepting the challenge, 'a', is 1 (out of 3). What value of 'd' do we need for S to go to zero? Do a little interpolating again and you get d = 0.8668653 - that is, if each person who pours a bucket of water over their head inspires (on average) 0.867.. people to do the same, then everyone in the world will participate, within 26 days of the challenge starting. For a = 2, we need d = 0.4046021 for S goes to 0. And so on...

What's a plausible inspiration rate? For a normal person, probably zero, while for a celebrity... I don't know. But on average 'd' is probably very close to zero. I mean, we know for a fact that significantly less than the entire population of Earth has been nominated/taken part in the challenge.

If you plot I(n) for a = 1 and d = 0.1 you get something like this,

In this case, you have that same rise and fall - but this time, you have a longer tail, like we see in the YouTube searches. Is this proof that this iteration of the model is more accurate? Maybe. But as I pointed out before, the YouTube searches don't necessarily accurately represent the number of people taking the challenge over time.


Social Pressure

When a friend does a thing for charity, then publicly calls you out to take part, there's a certain amount of social pressure to comply. I mean, if a friend dumped water over their head (just because), then asked you to do the same thing, probably you'd look at them like they were a crazy person. Anyway, that kind of social pressure is implicitly included in the 'infection' rate - more social pressure, bigger 'a'.

Instead, the sort of social pressure I'm talking about here is the kind that goes: "I should accept the challenge because so many other people have already done it". Or alternatively, "it's okay for me to accept the challenge, since so many other people have already done it".

In other words, the more people accept and complete the challenge, the more likely a nominated person is to accept too. Mathematically, we can introduce this with a term in 'S*I*R'. Or alternatively, we can keep the term as 'a*S*I', but make the factor 'a' a function of R -> a(R).

Anyway, if you're interested, you can try investigating that yourself. Or try adapting the model in some other ways. But beyond a point you can end up complicating a model more than improving it.


A Network Theory Approach to Nominating

So I was eventually nominated for the challenge. But I'm a wimp, so I declined to dump ice water over my head, instead making a donation to the Motor Neuron Disease Association (the UK equivalent of ALSA).

For my nominations, I wanted to try and maximise spread. So I nominated the 3 of the people in my Twitter network who are the most active and well connected, and who I thought would be up for accepting the challenge. Plus, as a secondary effect, I figured they might nominate other people in my Twitter network - the network theory equivalent of wishing for more wishes.

In the end, one ignored the nomination, one acknowledged but didn't accept, and one accepted (in the form of a donation) but didn't nominate anyone else. So I guess that theory didn't quite pan out. But I did at least encourage more charitable giving.


So Yeah

If we had real world data we could maybe test the accuracy of these models. But even without, we can get a sense of how the challenge behaves - for example, we find that there's a critical 'challenge acceptance ratio' that determines whether the viral campaign will go 'pandemic'.

In theory, you could apply this sort of model to any viral campaign, or just anything that spreads 'virally'. The nice thing about the Ice Bucket Challenge in particular, though, is that it has well defined rules for how the challenge spreads from one person to the next.

So, yeah..


Oatzy.


[I need an editor, my pronouns are all over the place.]

How Long Will a Prime Number Be?

File this under things I think about when I'm trying to get to sleep.


The Question

The idea is this - pick a random integer between 0 and 9. If it's prime, stop. If it's not prime, pick another random integer and stick it on the end of the first number. If this new number is prime, stop. If it isn't prime, pick another random integer, and so on.

So, for example:

1 -> not prime
2 -> 12 -> not prime
7 -> 127 -> prime -> stop.

The question is this - on average, how long would we expect the number to be when we get a prime?


Sometime Programming Just Won't Cut It

Normally, this is where I'd bust out Python and write a little program to run through the procedure a few times and see what happens. The problem is, testing whether a number is prime tends to be really slooooow.

The 'easiest' algorithm, trial division - literally checking if the number is divisible by any of the integers less than it - has (worst case) speed O(sqrt(n)). What this means is, each time we add a digit to our number, it will take ~sqrt(10) = 3.16 times longer to check if it's prime. And that time soon adds up. There are quicker algorithms, for example AKS, but they're too complicated to bother with.

Instead, we can figure out the answer with some good old-fashioned maths.


Some Good Old-Fashioned Maths

Lets re-state the problem - what we want to know is the expected length of the prime number. Here we use the standard definition for the expectation value

So now the problem is find p(N), which in this case is the probability that a number of length N is prime. Now, as any romantic mathematician will tell you, prime numbers are mysterious things, which are scattered about the number line seemingly at random. However, what we do have is the 'Prime Number Theorem', which estimates the number of prime numbers smaller than n, as

where ln(n) is the natural logarithm. Here's a Numberphile video explaining the theorem.

So, for one digit numbers (numbers less than 10) there are ~10/log(10) = 4.2 prime numbers (well, technically 4, but this is an estimate). So the probability that a one digit number is prime is estimated as 0.42. The estimate gets better as the number of digits gets bigger.

Now, for numbers of two or more digits, the number of primes that are N digits long is (approximately) the number of primes less than 10^N minus the number of primes less than 10^(N-1). So then, the probability of an N-digit number being prime is the number of N-digit primes divided by the total number of N-digit numbers. It can be shown (homework) that this probability is given by

Plotting this function, we see that the probability goes to zero as N goes to infinity

In other words, the longer the number gets the less likely it is to be prime. In this case, you can probably guess what the expected length of our prime will be.

So we have

And as N->infinity, N p(N) -> 1/ln(10). What this means is the sum doesn't converge - it goes to infinity. So on average we expect our prime to be infinitely long.


So Basically

To cut a long story short, the procedure is most likely to stop when the number is only one digit long. Or else, the longer the number gets, the less likely the procedure is to stop. This is another reason why the programming solution was a non-starter.


Oatzy.


[Now stop doing maths and go to sleep.]

Modelling the Survival of Chocolates

So it's been about a year since I last blogged. I've been kinda busy with university stuff, and growing my facial hair into something approaching a beard. I might write a blog about it at some point.

Anyway, I've got a little free time now - at least until after Christmas - so here we are. Let's waste some time...


"This is a real paper, published in a proper journal"

So the other day, my dad hands me this print out of an academic paper, "this seems like your sort of thing". The paper is called

"The survival time of chocolates on hospital wards: covert observational study" [1]

The basic set-up is this: staff at various hospitals put out boxes of Roses and Quality Street chocolates on wards, and (covertly) recorded every time one was taken. From this, they analysed the 'survival' of the chocolates, according to Kaplan-Meier survival analysis.

They found that the rate of chocolate consumption fit an exponential decay, with the median survival time for a chocolate being 51 minutes, and the survival half life ("time taken for 50% of the chocolates to be eaten") being 99 minutes. They also found that Roses chocolates were preferred over Quality Street.


"Damn it, man, I'm a physicist, not a doctor!"
 
So I'm reading this paper, and I'm thinking "I reckon I could model this". So I had a crack at it.

Basically, what we've got here is an agent-based modelling problem. Individuals encounter a box of chocolates at random intervals. Each box has a variety of different types of chocolates. So each individual will look through the box for any chocolate(s) they like the look of, and help themselves.

We can make a few simplifications to this. First of all, we work in fixed, discrete time intervals. In each time interval, some number of people will encounter the box of chocolates - one person per interval in the simplest model, and some random number of people in a more complex model.

Now, different people have different preferences. We could generate for each person a random set of preferences, but for simplicity we'll assume some 'average preference probability'. In other words, imagine there are equal numbers of each type of chocolates - what is the probability that any given person will pick type 1, type 2, etc. ? 

So, in the model, each person will pick a random chocolate based on this preference probability distribution. This accounts for the fact that some types of chocolate (strawberry) are naturally more popular than others (toffee penny).


Here's Where it Gets Maths-y

We have a box of N chocolates total. There are m different types of chocolate, with each type i=1,..,m having a count of ni such that 

Each type also has a 'preference probability' pi with 0 <= pi <= 1 and

In each time step dt an individual randomly chooses a type of chocolate, i in {1,...,m}, based on the preference distribution, and one chocolate of that type is removed from the box according to 

This process is repeated until N(t) = 0

In the more complex model, some random number of individuals (following a Poisson distribution) each choose chocolates in each given time step. This is also functionally equivalent to a single individual choosing multiple chocolates in a single given time step.


What Differential Does It Make?

In the simplest model, one person chooses a single chocolate in each time step. This is nice, because it means we can express the model as a set of differential equations, which can be solved exactly.

For each type of chocolate we have the rate of consumption given by

which has standard (exponential decay) solution 

And since the total count N is just the sum of the individual counts ni we get the full solution

In other words, it's just the sum of several exponential decays (with different decay rates) - making it, roughly, an exponential decay itself. 

In particular, in the special case where all the chocolates are equally preferred such that pi = p for all i=1,..,m we get the solution

And all that is loosely what was observed on the hospital wards. 

But these 'exact' solutions don't capture some of the subtleties of the real life set-up. For example, the toffee pennies might be completely ignored until they are the only option left; whereas, the exact solution assumes that (fractions of) all of the types of chocolate are being taken all the time, just at different rates.


"We did not seek ethical approval for this study..."

As an improvement on the analytical solution, we can create a computer simulated model, that captures more of the randomness of real life (code below). For the smplest model, the results look something like this

Notice that in this version of the model, the number of chocolates decays more slowly, and has some of the granularity seen in the observed data plot. 

Plus, the randomness means the specific behaviour varies from run to run - for example, the orange line decays relatively quickly, while the blue has a long stretch when nothing is being taken.

Interestingly, this model seems to have a linear decrease for the first ~20-30 time steps in every sample run. This is probably because, for those first few steps, everyone can get what they want. But eventually, certain types run out, and the rate of consumption decays away.


Why Do Chocoholics Come in Threes?

I wrote about the Poisson distribution before, back when I was modelling cafes. Basically, while we might see an average of one person per time interval, we don't expect to see them so evenly spaced out. Rather, people will arrive in 'bursts' of varying size. This behaviour is described by the Poisson distribution.

So I wrote my simulation in Python, using the numpy library for the Poisson distributed random numbers. The results look something like this

In this case, the initial decrease is less linear than in the simple model, but the overall decay times are roughly the same (except for the blue line above). This is what we'd expect since, in the above, I chose the average number of people per interval to be 1.

So yeah. Try the code, have a play with the parameters, see what happens. The output is the total number of chocolates remaining after each time step.

One thing to look at would be how changing the probabilities affects the rate of decay - for example, does the number of chocolates decay faster when all the chocolates are equally popular, or when some are much more popular than others. I suspect the former is the case.

It might also be worth adding an option where a person can pick again if their prefered chocolate isn't avaiable.


And the IgNobel Prize for Medicine goes to...

So that's that. And I can definitely see the original paper being nominated for an IgNobel, at the very least.

As for me, I have some gift wrapping I really need to get to...

Merry Christmas!


Citation

[1] P R Gajendragadkar, et al. "The survival time of chocolates on hospital wards: covert observational study" BMJ 2013;347:f7198 (Published 14 December 2013)


Oatzy.


[Next job, modelling Her Majesty's nuts...]

Is There a Formula for The 'Perfect' Christmas Tree?

The other week I was waiting for my tram, reading Metro, and I see this article - Treegonometry: Formula for perfect Christmas tree discovered

In other word, a set of equations that will tell you, for example, how much tinsel, or how many baubles are 'required' for the 'perfect' Christmas tree.

Here is the press release, with the actual formulas.

Now, I have a few issues with this. For one thing, it's a bit silly. Surely the perfect tree is a matter of personal preference, and in that sense can't be defined by a set of mathematical rules.

Plus, it seems kind of arrogant to declare that their way is the 'perfect' way, and that (implicitly) any other way is wrong.

Also, the fact that the equations are preposterously simple (perfect thing = constant*tree height) makes me suspicious of them. But, worse than that, it really bothers me that I have no idea how they were derived.

But I'm starting to rant.

The equations in question were actually from a 'study' by the Sheffield University Maths Society (student from my university), and were commissioned by Debenhams (department store). So, I dunno, maybe I'm just bitter that no-one's ever commissioned me to do maths.

Anyway. If they'd said that these were equations for an 'ideal' Christmas tree, then I'd consider that more reasonable. You can say, for example, that the ideal Christmas tree should strike a balance between too many and too few decorations, etc.

So, to that end, I'm going to have a crack at deriving some, more general, 'Equations for an Ideal Christmas Tree' of my own.


Lights

Before we start, there's one simplification we need to make - we will assume that Christmas trees can be approximated to a smooth, regular, right-circular cone - height h, and base radius r.

There's going to be a certain degree of error as a result of this approximation, but I'm not going for perfect.

Okay. So, for lights, there's actually a justifiable basis for an ideal - the ideal tree should have its lights evenly distributed.

Or, to put it another way, you ideally want to wrap your lights such that you have a (roughly) constant light surface density (lights per unit surface area); because you don't want your tree to be covered in clumps and bald spots.

This seems like a vague point; how do you decide on the correct 'density'?

Well, we actually have two constraints to work from: firstly, we have to wrap a single line around our conical shape, and second, the distance (s) between each pair of adjacent lights on a string will be fixed.

From this, we can layout a (hypothetical) grid on the surface of the tree, so that all the lights are equally spaced out.

So, how do we do that?

1) Draw yourself a cone (tree)
2) Draw a grid
3) Draw lights at the points where the grid lines meet
4) Draw in the proper row lines - these represent the actual string of lights

In this cases, I've drawn the grid a little tight, so you'd probably want two strings of lights to get this layout (hence the two colours).

And if you like, you can play with the distances between rows, and the angles, to get some slightly different layouts

The best grid to try for would be a triangular grid, where the distances between each light and the six nearest are the same

So, let's say we've picked a grid. How many light will we need?

There's a nice short-cut to working this out, without having to worry about lengths of spirals on conic surfaces. To do this, we take advantage of the regular grid layout to work out the lights density (lights per unit area).

Each section of grid is roughly a diamond, with all sides about the same length

So with a bit of trigonometry, you can get an approximate equation for the area of this diamond

Then, the density is one over that area (since there's one light per grid diamond).

We then multiply this density by the total conic surface area (excluding the base) to get the total number of lights needed

And, finally, multiplying that by the separation between lights (s), we get the total length needed:

Easy. Though, somewhat more complicated than the L = pi*h equation, derived by SUMS.

So, then, for the triangular grid (theta=60), you get

Cool.

The only problem then is getting the lights to line up on a grid on the actual tree. So... good luck with that.

Realistically, all of this is mostly irrelevant anyway - you can't go out and ask for, say, exactly 5.83m of lights; you buy your lights in pre-cut lengths. Can't get hold of the 'perfect' length of light? Then your tree is imperfect, and you should feel bad.

Anyway. The point is, if you can get your lights more or less even - so that there aren't any clumps or bald spots -, then as far as I'm concerned, you're on to a winner.


Tinsel

What's interesting about their equation for tinsel is this factor of 13/8. Now, again, I don't know how these equations were derived, so I don't know if this was intentional; But, 8 and 13 are consecutive Fibonacci numbers. Why is this important? Well, the Fibonacci sequence is closely related to spirals.

In particular, there's this thing you see in nature; for example, if you look at the spirals on a pineapple, the spirals going in one direction might be 13 and the number in the opposite direction might be 8. Or the numbers might be 21 left and 13 right... The point is, the numbers of spirals on a pineapple are always consecutive Fibonacci numbers (or sometimes Lucas numbers).

And you get the same effect on other things, like the spirals of seeds in the head of a sunflower, or the seeds on a strawberry, or the spirals on a pinecone, or a cauliflower, or all sorts of things. Hell, maybe the branches on a Christmas tree form Fibonacci spirals.

Vi Hart explains it better than me.

So maybe that has something to do with that pre-factor. Or maybe not. It's an interesting tidbit, though.

Anyway.


The thing with tinsel is different people like to do tinsel differently - some like to elegantly drape it across the outer branches, others like to wrap up their tree light they're restraining a hostage. It's a matter of preference. But it's going to affect the amount of tinsel you'll need.

Where tinsel differs from lights is, you're not trying to set up a grid, or get an even surface density. Rather, in this case, you'd probably want to wrap it such that the rows are more horizontal, with a roughly constant vertical separation.

Here's where things get messy; the equation for the length of a spiral on the surface of a cone is given by

I know, right? Maybe you would be better using the SUMS equation for this one.


Extras

The star/angel is going to be some fraction of the height of the tree.

I don't know what the ideal value of the fraction (alpha) would be, but the 10th they came up with seems reasonable.

Baubles, I haven't a clue how they came up with those numbers. The factor of sqrt(17) makes me think some geometry was probably involved, but I dunno.

You would probably want to figure it out as some ideal ornament surface density, Db (like with the lights). In this case, the number of baubles needed would be something like

In fact, if your baubles are all, more or less, the same, you can lay them out on a grid, like the lights. Though, this time, you'd want to have them a little more spaced out, since baubles are much bigger than fairy-lights. But at least this time you don't have the separation constraint.


Thoughts

Here are various things Ben Goldacre, of Bad Science, has said on the subject of commissioned, 'perfect' formulas. Here is an article by mathematician, Simon Singh. Here is an article on BBC News. And here is a collection of such formulas on Apathy Sketchpad.


To be honest, I wouldn't bother with any of this; their equations, or mine. I mean, would you really want a tree that was 'perfectly' decorated? Cold and artificial are the words that come to mind.

And, frankly, I'm not sure my equations would actually work in practice.

I'd say, use your best judgement on how much of everything you'll need, and just do your own thing. Have fun with it!


Our Christmas tree is imperfect. In fact, it's gloriously imperfect. No, seriously, it's a mess.

My ex's family used to construct these massive, elaborate, works-of-art trees; with yearly colour schemes, and matching baubles, and everything. By comparison, she described our tree as kitsch.

But it's adorned with all the baubles, and tinsel, and decorations we accumulated over the last 20-odd years; at least, the ones that haven't been lost or broken. And, in a sentimental sort of way, it is perfect.

Well, okay, not perfect. But, damn it, it's ours.





Oatzy.


[I want you to know, I had no part in decorating that tree.]

[...And, yes, that's a weeping angel on top.]

Picking a Seat on the Tram

When I'm at university, I have to get the tram twice a day. Now I quite like the trams in Sheffield (at least, more than the trains and buses). But if there's anything I can do to get a little more leg room or personal space, then I'm going to try it.

On the way in to university, I get on the tram at the start of the line - i.e. the tram is more or less empty when I get on, so I have pretty much free choice of where I sit. In the main, middle section of the tram, seats are grouped into fours - two forwards, two backwards, face to face - and there are there are 5 rows and 2 columns of these groups.

(The are more seats in the front and rear carriages, but with difference layouts. And I almost never sit in those sections, so I'm ignoring them.)

When the tram's busy, you have to accept that you're probably going to be squashed, and bumping knees with the person across from you, and desperately trying to avoid awkward eye-contact. But when it's less busy, there's the chance for a bit of personal space.

The question is, where's the best place to sit (assuming you have free choice) so you're more likely to get some of that personal space?

Truth be told, I don't actually know. I figured, probably the best place is the forward-facing window seat in either of the 3rd (middle) row groups.

The logic was that a person would either have to sit next to you, or so they're facing backwards - generally less preferable options. And since people are more likely to grab seats nearest doors, the 3rd row is 'best' because it's equidistant from both doors.

Of course, if the front and rear doors aren't used equally, then the choice of row would have to be tweaked.


From that, I started thinking about how people chose which seat to take in a given four-group. And since I have nothing better to do than stare out the window for the 20-odd minute journeys, I figured I'd try to model it.

From empty to full there are 16 possible seat use configurations ('states'). People chose seats based on personal preferences, as well as which seats are already taken. So, for the model, we look at the different states, and consider how they might change with the addition or removal of various numbers of people.

Here's an example, that I sketched, of possible transitions resulting from the addition of one person at a time

The red lines are possible changes of state where the number of people stays the same (usually after the removal of others). More on that later.


Now, the key thing here is that different states, and different transitions are more likely to occur than others.

For eample, a person with free choice is more likely to chose a forward-facing seat (some people tend to feel unwell traveling backwards), and a window-side seat so they've got something to stare at; unless they're planning on making a hasty exit.

So that would theoretically make the backward-facing aisle seat the least popular.

HOWEVER, if someone's already sat in FFWS, then BFAS becomes the more preferable, since it avoids having to sit next to, or directly in front of, some random stranger whose just staring blankly out the window in an unsettling sort of way. Also, personal space.

And from there, if you introduce another person, they're most likely to sit FFAS, since BFWS is a bit awkward to get to and leaves you feeling kind of boxed-in. Plus you have the extra arm room on the aisle-side, and, if necessary, you can turn slightly to give yourself a little more leg room.

And if you add a four person, they have a hell of a time getting to the free seat, but it's their only option. Unless BFAS moves across to the window.

Now, back to those red lines in the diagram. From the 3-config described above, say BFAS leaves. You now have two people sat next to each other. The person in FFAS may then chose to move to BFAS for the sake of more personal space. This is kind of like how physical systems tend towards their lowest energy (least socially awkward) state.


The model itself is actually pretty straight-forward. What you do is construct a stochastic matrix (sometimes called a transition matrix) - basically, a 16x16 'table' that tells you the probability of the seating changing from one state to another. For example

And this works for all possible transitions between all possible states.

The tricky part if determining the probabilities. You could just make random guesses at it. Or, if you were really determined, you could spend a load of time on trams, recording what transitions happen at each stop and how many times they happen. I don't plan on doing either.

But say you have your data and you've constructed your matrix, M. You could do a Monte Carlo type simulation using those probabilities. But with a stochastic matrix you can be more precise.

See, if you calculate M*M, then the entries are the probabilities of moving from state i to j after two 'stops'. And if you calculate M^n, they you get the transition probabilities for after N stops.

And if you add up the entries down each column, j, then you have the probability of being in state j after N stops - that is, you can find the most likely seating configuration after some arbitrary number of stops. You could also use this method to work out the popularity of each seat over all possible states.

Which I think is pretty cool.


And that's all well and good for the groups. But what if you want to model the whole tram? That's where things get tricky, since you have to model how people chose seat groups, which itself depends on what seats are already taken in groups.

Ultimately, it can be modeled the same way, with a transition matrix. Only now, you're working with a 40x40 matrix.

Not difficult, per se, but certainly requires a lot of data and number crunching...

Come to think of it, if you were going to go to the trouble of observing and counting state transitions, you could just count how many times each seat is sat in over some arbitrarily long period to figure out each seat's popularity. Then you just need to try to get the seat next to the least popular.

But counting isn't as fun as modeling.


Oatzy.


[I wonder what normal people think about on public transport..]

Sharing a Burger Between Three

There are several ways to divide a burger evenly between three people.

Probably the best is to cut it radially (like a pizza). It can be tricky working out exactly where to make the cuts, but if you can pull it off, then all the pieces will be roughly identical (topping distribution notwithstanding).

But for the sake of arguing, lets say you want to divide the burger by making two parallel cuts: Where, then, do you make the cuts so that all three people get the same amount of burger?

NB/ This gets quite maths-heavy, so if you're not interested in that sort of thing, feel free to skip right to the end for the solution.



Geometry

For simplicity, we're going to consider the burger as a circle, and make the cuts so that each chunk has the same area. The two cuts are going to be the same distance from, and parallel to, the central axis of the burger, so we only need to consider the position of one of the cuts.

Here's the set-up

We work out the area of the cut-off as the area of the circular segment, minus the area of the triangle.


- Aside: Radians

Radians are basically an alternative way of measuring angles. For maths and physics they're generally more useful than degrees.

They're relatively easy - there are 2pi radians in a full circle, so 2pi radians = 360 degrees

1 radian = 180/pi = 57.3 degrees
1 degree = pi/180 = 0.017 radians, etc.

*    *    *

Back to the circle; with angle x in radians, the area of the circular segment is

The area of a triangle is half base times height..


- Aside: Area of the Triangle

We start by splitting the triangle down the middle, so that we have two identical right angle triangles

The height, l = r*cos(x/2)

The base, b = 2*(r*sin(x/2))

So the area of the triangle is (l*b)/2 = r^2 sin(x/2)cos(x/2)

Finally, use the identity sin(2x) = 2sin(x)cos(x)
to get

*    *    *

So the area of the cut-off is

and it needs to equal a third the area of the circle = 1/3 pi r^2.

So first, we need to find x satisfying

or

Once we have a value for x, we find where to make the cut from

Intermission

The thing about this equation is it doesn't have an exact, analytical solution - to find the solution you have to use numerical methods. Well, I say you have to use numerical methods; these days you can just type the equation into WolframAlpha, and you'll get a solution like *snaps fingers*

Which is nice. I even have the WolframAlpha app on my phone. But when I thought up this question I was on holiday in Sherwood forest, where there was literally no mobile singal.

So that was out of the question. And since I'm not in the habit of carrying a scientific calculator around with me, I was stuck with the basic calculator on my phone. It looks like this:

No trig functions, no square roots, no pi button. It doesn't even do brackets, or have a memory function. Luckily, I am in the habit of carrying around a notepad and pen.

Anyway, there are two ways of working this out with only a basic calculator. The first is 'easier', but only if you know some stuff, and the numbers happen to be nice (in this case, they kind of are). The second is harder, in that it requires more number crunching, but it'll work with any numbers, and can be more precise.

Again, feel free to skip to the solution if you're not interested in the gritty details.



Method One

First of all, here's a graph of the two sides of the equation

hand-drawn with Skitch

We want to find the point at which the two graphs cross. We can see that that happens somewhere between 2pi/3 and pi (120 and 180 degrees). So, lets make a guess that it's exactly halfway between these two values: 5pi/6 (150 degrees).

For the right hand side of the equation: 5pi/6 - 2pi/3 = 0.52

NB/ I'm using pi=3.1416 (rounded to 4 decimal places). If you prefer, you could use the approximation 22/7. The result should be roughly the same.

For the left hand side of the equation, we need to work out sin(5pi/6)

At A-level, we were expected to memorise sin() and cos() of angles 0, 30, 45, 60, 90, and 180 (degrees). We were also expected to know the formulas for sin() and cos() of sums of angles. For sin(), it works like

Why is this important? Well 150 degrees = 180 - 30 (5pi/6 rads = pi - pi/6)

So

And since 0.5 is pretty close to 0.52 - less than 5% error - we can accept the convenience of that answer and say it's close enough.

So our approximate value of x is 5pi/6 = 2.618

[Incidentally, the identity for sin(2x) is just a special case of the above, with a=b=x; i.e. sin(2x) = sin(x+x) = 2sin(x)cos(x)]


Now we just need to work out l/r = cos(x/2) = cos(5pi/12)

For this one, 5pi/12 rads = 75 degrees = 45 + 30, so we can use

So

And we just have to evaluate that. But we don't have a square root button. Now, I just happen to know that sqrt(3) ~ 1.73 and sqrt(2) ~ 1.41.

But I'm just weird like that. Let's say you don't. How do you work it out?


- Aside: Square Roots

There are several ways of working out square roots with just basic operators. For two easy examples:

The first is 'Trial and Improvement' - pick a number, square it, does that give the right answer? If not, pick another number based on whether the last guess was too big or too small.

For example: sqrt(3)

1.5 -> 2.25 -> too small
1.7 -> 2.89 -> too small
1.8 -> 3.24 -> too big
1.75 -> 3.0625 -> too big
1.73 -> 2.9929 -> too small
1.74 -> 3.0276 -> too big
1.735 -> 3.010225 -> too big
1.7325 -> 3.00155625 -> too big
1.732 -> 2.999824 -> too small
etc.

The second method is the "Babylonian Method". It's more systematic, and can converge to the correct answer quicker than guessing. But it can be irritating if your calculator doesn't have a memory function.

It uses the recurrence relation

Basically, you make a guess xn. Divide the number you want to square root (S) by xn. If xn is lower than the actual square root, then S/xn will be greater than it. That means the actual root will be between xn and S/xn, so we make the next guess xn+1 the average of these two values. Repeat until x is sufficiently accurate.

For example: sqrt(2)

x0 = 1.5 -> 2/1.5 = 1.33
x1 = (1.5 + 1.33)/2 = 1.4166.. -> 2/1.4167 = 1.41176..
x2 = (1.4167 + 1.41176..)/2 = 1.41421.. -> 2/1.41421 = 1.41421..

*    *    *

Whatever way you do it, you repeat the process until you get the degree of accuracy you're happy with.

You should get the answer around l/r = 0.259



Method Two

We go back to the equation sin(x) = x - 2pi/3

We still have to find the solution numerically, we still don't have a calculator with a sin() function, and this time the numbers don't work out nicely.

So, the question is, how do we calculate sin(x)?


- Aside: Taylor Expansion

The Taylor Expansion of a function is a way of fitting a polynomial (sums of powers) to a more complicated function. It works like this

Basically, it gives a way of converting a function we can't calculate into an infinite sum of powers of x, which we can calculate.

It's usually expanded around the origin (x0=0), since the equations work out neater. But you can do it around any point, x0=a. This is useful if the value you are trying to calculate is far from x=0. The closer x is to x0=a, the quicker the sum converges.

Even though the expansion is an infinite sum, it's usually sufficient to just take the first few terms, since each additional term makes a smaller and smaller contribution to the sum.

So the trick is working out how many terms you need to include to get some desired level of accuracy.

*    *    *

In this case, I'm going to use the Taylor Expansion of sin(x) around x0=pi, since the approximate value (2.6) is nearer to pi than 0.

Here's what the expansion looks like

So, how many terms do we need to include?

Here's what the graph looks like for different numbers of terms

plotted with WolframAlpha

For a value around 2.6, it can be shown that including the first two terms is correct to ~3 decimal places; the first three terms is correct to ~5 decimal places; the first four terms to ~7 decimal places, etc.

So I would probably go to the third term (for 5dp), but only take the result to 3dp.

NB/ We shouldn't get too hung up on getting an extremely accurate value for x, since we're already getting rounding errors from the factors of pi in the expansion. Also, since we're calculating x to 5dp, we should use pi=3.14159

That means we want to solve

which can't be solved exactly.

So, for finding the correct value (without WolframAlpha), we can use any root-finding method. For what it's worth, I used Trial and Improvement; the other methods are easier with a computer.

But, note that the function is decreasing

So if the guess gives a value greater than zero, you need to increase the value of x (and vice versa).

If you run through all that (I won't go into detail), it gives a value around x=2.605


Alternatively, you could expand around x0=5pi/6 (if you know/can work-out sin and cos of 150 deg without a calculator).

In this case you'd only need up to the term in x^2 (correct to ~4dp). Using this expansion would mean solving a quadratic equation, which is easy. But using this expansion can introduce more rounding errors from the factors of sqrt(3). It's a matter of preference, I guess. The answer should be about the same.


Finally, we need to calculate cos(x/2)

Again, we use the Taylor Expansion to calculate cos(). In this case, we're doing the expansion around x0=0; the expansion is

In this case, you just keep adding terms until the result remains approximately constant to some desired degree of accuracy (3pd).

This gives a value around l/r = 0.265



So What is the Real Answer?

Once I got to somewhere where I could get at WolframAlpha, I checked the real numbers; here are the results:

The approximation of x from Method One (2.618) is an over estimate by ~0.5%, which is relatively acceptable. The approximation from Method Two (2.605) is correct to 3 decimal places, which is definitely acceptable.

And for the value of l/r

From Method One (0.259), the approximation is an under estimate by 2%, and correct to 2 decimal places, so is probably acceptable. The approximation from Method Two (0.265) is, again, correct to 3 decimal places. So that is also acceptable.

So, if the numbers happen to be convenient and you know some trigonometry, you're probably as well using Method One. If not, or if you just want more accuracy, then go for Method Two.



Applying the Results

The results are actually quite nice, in terms of practical application (dividing up a burger). The ratio of the radius (0.265) being close to one quarter, you find the cuts like this

That is, find the central axis, then find the (imaginary) line halfway between the centre and the edge - make the cut halfway between the centre and this imaginary line (maybe cut an extra hair's breadth towards the edge). Repeat on the other side.

Easy.


So, now you know. Obviously, all this applies to dividing any circular thing evenly between three people. You could probably even adapt the methods for sharing between even more people.

And in theory, You could do all this with just pen and paper (no calculator). Though you probably wouldn't want to. I know I wouldn't..


Oatzy.


[Wow, I really managed to stretch that one out.]