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...]

Sexuality and Population Dynamics

So I was reading this response to a Daily Mail opinion piece about those 'gay cure' bus adverts, and I see this quote (from the Mail article)

“Since only about one percent of us are [gay], homosexuality is obviously a departure from the norm... Reversing that proportion would spell the end of the human race, which is clearly undesirable.”

There are several things wrong with this statement - some of which the response, linked above, points out.

First of all, that 1% values. The actual number is a tricky figure to pin down. The UK Office of National Statistics puts the figure at 1.5% gay or bisexual; the linked response lists 3 studies which put the value around 7% (for the UK), and various other studies suggest number in the modern West is somewhere between 2-13%.

Secondly, just because something is (statistically) "a departure from the norm", doesn't mean it's bad - for example, ~10% of people are left-handed, ~2% of people have red hair, etc. (Historical beliefs notwithstanding.)


But all that's tangential. What I'm interested in is the second part - would a majority gay population spell the end of humanity?

[NB/ for the remainder of the blog, I shall be using 'gay' colloquially, as shorthand for all of LGBTQUA, etc.]


The Simple Model

It's an interesting population dynamics question, and it's easy to model if we make some assumptions and simplifications:

1) We assume that homosexuals don't reproduce at all (either directly or indirectly). This isn't true - but it seems to be a common belief.

2) We assume that the rest of the population reproduces with some fixed rate (b=0.02059), which doesn't change from year to year. Similarly, we assumes a fixed death rate (d=0.00812), which affects all members of the population equally.

3) We assume that sexuality is binary, definite and static, with a fixed incidence of homosexuality, g = 0.07 - that is, of all children born 7% will be gay.


So the model we're using is preposterously simplified - effectively,

new population = old population + births - deaths

And we can easily program this model and run simulations.

Now, there are two ways to interpret the article's statement in the context of our model -


1) "if 99% of the (initial) population were gay" assuming the 7% incidence

The outcome looks like this

In this case, the population drops at first, but then starts to grow again, exponentially - certainly not the end of the human race. Also notable is that the proportion of the population that is gay stabilises at 7% (as expected).

2) "if 99% of people were born gay"

Now, the outcome looks like this

Exponential decay - so, in this case, humanity would die out. But if 99% of the population isn't reproducing, what do you expect?

Here are the phase portraits of the above two cases, for a slightly different view of what's happening (gay population along x, straight population along y)

Again, for g=0.07 the population tends to infinity; for g=0.99 the population tends to zero. Note that this behaviour is (ultimately) the same for all initial population ratios.


What I Learned in MAS271

We don't actually need numerical simulations to study the behaviour of this (or similar) systems - the model can be solved analytically.

First of all, we generalise the system. We start with a pair of coupled, first-order differential equations describing the behaviour of the two populations

Now, these two populations are mutually exclusive subgroups of the total population - that is, total population z = x + y. So with a little mathematical trickery, we can combine these two into a single equation, describing the behaviour of the total population

This is a second-order, linear differential equation, so has the general solution of the form

Where lambda values are given by

And A-coefficients are constants, which are determined by the system's initial conditions. If the initial total population is z0, and the inital size of the x population is x0, then we get

In other words, we now have an (analytical) equation which completely describes how the population will behave for any given set of variable.

Incidentally, the equations for x and y have the same form as z. The only difference is the A-coefficients.


Boom or Bust

What we want to do is be able to determine the long term behaviour of the population based on the given set of variables - that is, will the population grow, or die out?

For the coupled system above, there is only one stationary point, at (x,y) = (0,0) - i.e. when the total population is zero. What this means is, if the population is zero, it will stay zero. If it has any other value, then the population size will change over time.

In terms of the stability of the system, if (0,0) is stable, the population will tend towards zero - i.e. will die out. If (0,0) is unstable, the population will tend to infinity - i.e. grow indefinitely.

As it turns out, the lambda values above are the eigenvalues of the coupled system, meaning their values determine the stability of the system:

i) if both lambdas are negative, the system is stable and the population will die out;

ii) if one or both of the lambdas are positive, the system is unstable and the population will grow to infinity.

Now, if we take the limit of z(t) for t tends to infinity, we can simplify the population equation

What this tells us is, if lambda1 is less than one the population will die out. If lambda1 is greater than one, the population will grow.

And we can use this to determine the conditions for instability - the conditions under which the population will grow. In this case the condition is lambda1 > 0. But this can be simplified down to

There's one last trick we can do - determine what proportion of the total population will be in the sub-population x (e.g. what what fraction of the population will be gay) as t tends to infinity.

For this, we determine the equation for population x (in a similar fashion to z), and take advantage of the limits to get the equation

Note that this equation is time-independant, i.e. the proportion will tend to some constant value.


Back to Basics

Just as a sort of check, we'll go back to the basic model. For this, we have

Which gives lambda and A values

Which gives the population equation

Which behaves like the graphs at the top. And this tends to

And gives

That is, the proportion of the population that will be gay stabilises at the pre-defined 'birth' ratio - which is what we saw from the numerical model. Note also that this value is independent of initial conditions and doesn't vary over time.


Simple Instability

Finally, for the simple model, the condition for instability (growth) is given by

For the numerical values chosen above, this gives the condition g < 0.606 - that is, if less than 61% of the population is 'born' gay, then the population won't die out. And 61% is a pretty high upper limit - a majority, in fact.

Similarly, we can use this inequality to find under what conditions g=0.99 is viable. In this case, it requires a birth rate 100 times the death rate.

So, if the death rate is 0.008 (i.e. 8 death for every 1,000 people), then it requires a birth rate of 0.8 (800 birth per 1,000 people). Now, for each 1,000 people, 99% is gay, meaning 10 straight people per 1,000. Of those, about half is female (5). In other words, each woman would have to give birth to 160 babies per year!

Clearly this isn't feasible. But for comparisons sake, consider bees or ants - these species are eusocial, meaning populations made of one fertile queen, and hundreds of (mostly) sterile drones. In this case, the high birth-to-death ratio is more reasonable, and the predominantly non-reproducing population is still viable.


A More Complex Model

In the simple model we assumed that only the straight population reproduces, and that both populations die at the same rate. So we make the following modifications:

1) First of all, lets tweak the model so that the gay population has children as well, with some separate birth rate, b1.

2) Of those born to gay couples, some proportion g1 will be gay (NB/ studies show that the majority of children raised by same-sex parents grow up to be heterosexual).

3) And for completeness, lets say that the gay population has a different death rate as well, d1.

Here's what the model looks like

And here are the coupled equations

Now, all the previously derived general equations and such from earlier still apply. But in this case, the variables can't really be simplified, and writing them out in full would just be a mess.

So if we take, for example, b and d as before, and b1 = 0.00103, d1 = 0.01218, and g1 = 0.03, here are the new phase diagrams for g = 0.07 and 0.99

For this model, the phase portraits are basically the same as for the simple model - i.e. the effect of using the modified model over the simple model is negligible.


Complex Conditions

So, for this model, the condition for growth looks like this

The important thing to note here is, if the condition for the simple model (g<1-d/b) is satified, then the condition for the complex model is necessarilty satisfied.

In fact, if we set g to its maximum possible value 1 (i.e. all babies born to straight couple are gay), then the human race will still persist so long as g1 > d1/b1.

And what this means is, so long as the birth rate (b1) is greater than the death rate (d1), then every baby born could be gay, the human race still wouldn't die out.

However, given the low birth rate for gay couples (compared to death rate), this might not actually be satisfied IRL. In which case, an entirely gay population is still not sustainable.

For the b and d values already chosen, we can use the condition to plot how the upper limit of g varies with g1

So for these values, the upper limit of g can increase by no more than 3.4%. So, again, the extra complexity doesn't give a significantly different result compared to the simple model.

In fact, the additional variables only become significant when (b1/d1)g1 > 0.15 (for the values above); or in general, when

So, if we keep b and d as before, g=0.07, and re-set b1=b, d1=d, and g1=0.06, then (b1/d1)g1=0.152, and the phase portrait looks like this

Here, the upper limit on g has risen to 0.67. (cf: 0.61 for the simple model). And, if we set the initial condition x0=0.07*z0, then the limit of x/z = 0.0697; just less than g.



tl;dr - So, would homosexuality ever mean the end of the human race? - Under some circumstances, yes. But in most (real world) cases, no. And even if a majority of the population were LGBT, the human race still wouldn't necessarily die out.

BUT, as a final note - remember that these are very simplified models of population dynamics. So take the results with a grain of salt. The science of sexuality is much more complex and much less abstract than this.


Oatzy.


[Putting my degree to good use.]