Population Growth Models using R/simecol, Part 1 : The world population

In this entry, as part of the series on dynamic systems modeling with R and simecol, we'll take a look at population growth models, our main focus being on human population growth models and how they tie into other theoretical frameworks such as demographic transition theory and carrying capacity. These different models are then fitted to data on world population spanning a period between 10000 BC up to 2015, using R and simecol. In this latter part we will finally introduce the FME R package.

Malthus model of population growth

In 1798 the Reverend Thomas Malthus published An Essay on the Principle of Population, sparking a debate about population size and its effects on the wealth of nations. His Essay is framed within the wider context of considerations about the perfectibility of the human condition and it was Malthus' response to the likes of Condorcet and others who argued, in the typical positivist style of the era, that human reason would eventually surmount all obstacles in the way of mankind's progress.

At the core of the Malthusian argument was the observation that populations grow exponentially but the means of physical sustenance can only increase in arithmetic progression at best, and since population growth is commensurate to the availability of means of sustenance, any momentary improvement in the living condition of the human masses that would lift the restrictions on human population growth would inevitably result in long term famine and aggravation of man's plight. In modern day system dynamics terminology, we would say that population growth has a negative feedback mechanism that prevents further population growth and this idea in itself is not without justification. I mean, the planet cannot sustain exponential human population growth indefinitely, can it?

Malthus was in fact aware of negative feedback mechanisms and cyclical behavior in population growth. To paraphrase Malthus: the increase in population meant a decrease in wages and at the same time, an increase in food prices (less food for more population). During this time, the difficulties in rearing a family are so great that population growth is at a stand. But at the same time, a decrease in the price of labor meant that farmers could employ more laborers to work the land and increase food production thus loosening population growth restraints. The same retrograde and progressive movements in human well-being ocurred again and again- a cycle of human misery and suffering¹. Being a clergyman, this situation was seen by Malthus as "divinely imposed to teach virtuous behavior"².

At any rate, what is known today as the Malthus model for population growth, given in (1), only considers the exponential aspect of growth, as given by the rate parameter \(r\). No doubt the Malthus model is the simplest population growth model. It characterizes the growth of a population as dependent only upon the reproductive rate \(r\) and the population level at any given time (\(N\)). Note that the negative feedback controls we spoke of earlier are conspicuously missing in this model. In this model, there is only positive feedback: higher population levels beget even higher population levels, and so on to infinity.

\[\frac{dN}{dt}=rN\qquad\qquad(1)\]

Verhulst or Logistic Model of Population Growth

Considering the ideas on population growth that Malthus expounded in his Essay, Adolphe Quetelet judged that Malthus failed to provide a mathematical explanation as to how population levels, whose growth was by nature exponential, would top off at some point and not continue to rise to infinity³. Quetelet himself and his pupil, Pierre Verhulst, would assume this task in their effort to apply the methods of physics, which had seen a huge success with Newton's work, to the social sciences, as was characteristic of Positivist thought in the nineteenth century⁴.

Making an analogy with physics, Quetelet argued that in addition to the principle of geometric growth of the population, there was another principle at work according to which "the resistance or the sum of the obstacles encountered in the growth of a population is proportional to the square of the speed at which the population level is increasing"⁵. This premise bears semblance to the mathematical description of an object falling through a dense fluid: the deeper the object descends into this fluid, the greater the density, and hence the resistance to its downward trajectory, until the object reaches a depth beyond which it can't descend further⁶. By applying this analogy to demographic growth, we can infer that there is a maximum population level. As a result, a population finds the cause of its eventual equilibrium in its own growth⁷. In modern literature, this model of population growth is given by the following differential equation:

\[\frac{dN}{dt}=r_{max}N\left(1-\frac{N}{K}\right)\qquad\qquad(2)\]

Let us examine this equation in more detail to understand its behavior. For small population levels, when \(N\) is much smaller than \(K\), \(N\) exhibits exponential like growth, as dictated by the \(r_{max}\) parameter, also known as the intrinsic growth rate. This is because the \(\left(1-N/K\right)\) term is still relatively close to 1. However, as \(N\) becomes larger, this last term starts to matter because it approaches zero. As \(N(t)\) attains this point, the growth given by the \(dN/dt\) differential peters out and the population levels top off asymptotically at \(K\). This assymptotical level \(K\) is known as the saturation point or the carrying capacity of the population. We can integrate equation (2) to obtain \(N(t)\), whose graph is as follows:

Functions whose curve has this characteristic S-shape like the Verhulst model of population growth are called sigmoid functions⁸. We encounter sigmoid functions in many contexts: as the cumulative probability distribution function of the normal and T-Student distributions, as the hyperbolic tangent function and as the activation function of neurons in mathematical neural network models. Then of course there's also the logistic family of functions, which occur in statistics and machine learning contexts as well as in chemistry, physics, linguistics and sociology (in the latter two contexts the logistic function is used to model language or social innovation/adoption processes). Hence the alternate denomination for the Verhulst model of population growth as the logistic model of population growth. Another feature of this model is that it has an inflection point - after the \(N\) reaches \(K/2\), the rate of population growth starts to decelerate until it eventually levels off to zero. At this halfway inflection point, the differential rate of growth is highest and curve has its "steepest" tangent line. Notice the symmetry of the curve around this inflection point.

The beauty of the logistic model of population growth lies in its simplicity (only two parameters) and the interpretability of its parameters. The intrinsic growth rate (parameter \(r_{max}\)) is the rate of exponential growth when the population is small and the carrying capacity parameter \(K\) is simply the maximum population level attainable. One important disadvantage of the logistic model is the fact that its inflection point is exactly half the carrying capacity \(K\), which severly limits the applicability of this model. Nonetheless, Verhulst himself used this model to fit census data in France (1817-1831), Belgium (1815-1833), Essex County (1811-1831) and Russia (1796-1827), all with relative success⁹.

Variants of the logistic population growth model

In seeking to improve the applicability of the basic logistic model for population growth, many authors have since proposed models with more parameters that still retain the basic sigmoid features of the logistic model and include one inflection point. Given the inflexibility of the basic logistic growth model about the inflection point, Blumberg introduced what he called the hyperlogistic function whose derivative is given by (3), which we will in turn call the Bloomberg model of population growth:

\[\frac{dN}{dt}=r_{max}N^\alpha\left(1-\frac{N}{K}\right)^\gamma\qquad\qquad(3)\]

The Blumberg model introduced two additional parameters to surmount problems raised by treating the intrinsic growth rate as a time-dependent polynomial¹⁰, while at the same time contributing some measure of pliancy to the inflection point, given by (4). I would add that equation (3) cannot always be integrated - Blumberg cataloged the various conditions on \(\alpha\) and \(\gamma\) that result in analytic solutions to (3) when explicit integration can be carried out.

\[N_{inf}=\frac{\alpha}{\alpha+\gamma}K\qquad\qquad(4)\]

The Richards growth model, originally developed to fit empirical plant mass data, is given by:

\[\frac{dN}{dt}=r_{max}N\left(1-\left(\frac{N}{K}\right)^\beta\right)\qquad\qquad(5)\]

\[N_{inf}=\left(\frac{1}{1+\beta}\right)^{\tfrac{1}{\beta}}K\qquad\qquad(6)\]

The inflection point in the Richards model is given by (6). Note that for \(\beta=1\), the Richards model is the original logistic growth model with the same inflection point. For extreme cases of \(\beta\), we have \(N_{inf}=e^{-1}K\) when \(\beta\rightarrow 0\) and \(N_{inf}=K\) when \(\beta\rightarrow \infty\).

A further (and logical) generalization of the Blumberg and Richards models leads us to the so-called generalized logistic growth function with five parameters¹¹ (equation 7), where \(\alpha\), \(\beta\) and \(\gamma\) are all real positive numbers. The inflection point for this model is given by (8), which makes sense as long as \(N_0 < N_{inf}\). If \(N_0 > N_{inf}\), then \(N_{inf}\) will never be attained, because the intrinsic growth rate is assumed to be positive.

\[\frac{dN}{dt}=r_{max}N^\alpha\left(1-\left(\frac{N}{K}\right)^\beta\right)^\gamma\qquad\qquad(7)\]

\[N_{inf}=\left(1+\frac{\beta\gamma}{\alpha}\right)^{-\tfrac{1}{\beta}}K\qquad\qquad(8)\]

Yet another growth model similar to the logistic growth model is the Gompertz growth model¹². The Gompertz model has been used to describe exponential decay in mortality rates, the return on financial investments¹³ and tumor growth. It's equation takes the following form:

\[\frac{dN}{dt}=r_{max}N\,\log\left(\frac{N}{K}\right)\qquad\qquad(9)\]

An interesting variant for this model's representation is to use two stock variables, one representing the mass or population that's growing (\(N\)) and the other representing the diminish rate or growth (\(r\)). So, we would have two differential equations¹⁴:

\(\frac{dr}{dt}=-kr\)	(10)
\(\frac{dN}{dt}=rN\)

It must be noted that all the models above imply exponential growth from the first time instant, where in fact, the effective exponential rate of growth is highest, as the population or mass levels are still far from attaining the carrying capacity. In some situations, it would seem that some populations begin in a "dormant" stage, where they are not exhibiting this sort of exponential growth. Such may be the case, for example, in societies before the industrial revolution or even before the agricultural revolution, when scarcity was the order of the day and high birth rates were decompensated with high mortality rates, resulting in population numbers that grew ever so slowly.

While going over the literature in writing this post, I found an interesting two-stage population growth model in Petzoldt (2008). He describes the situation in which a microorganism culture is dormant while it is adapting to its environment, after which it begins to reproduce and increase in numbers according to the logistic growth mechanism outlined above. In other words, in this model, growth is delayed or lagged. Petzoldt's two stage model for population growth is in fact a two-compartment model (see my entry on epidemological models, where I discuss compartment models), with one compartment representing the part of the population in the "dormant" stage and the other representing the "active" population:

\[\frac{dN_d}{dt}=-r_1N_d\]	(11)
\[\frac{dN_a}{dt}=r_1N_d+r_2N_a\left(1-\frac{N_a}{K}\right)\]

Hyperbolic models of population growth

Some authors have pointed out that logistic population growth models, although fairly accurate for short-term growth of populations on a scale of centuries, fail to accurately describe the long-term population growth of the entire human race on a scale of thousands of years. Namely, the logistic models do not account for the fact that the world's population has grown from 2 billion to 7 billion in the last 50 years alone, in comparison to which the population growth curve seems almost flat for time periods in the remote past. In a seminal work by Von Foerster, Mora and Amiot¹⁵, they proposed that according to the available data by 1960, the world's population growth could be described by the the following model:

\[\frac{dN}{dt}=\frac{C}{(T_c-t)^2}\qquad\qquad(12)\]

This model has one essential flaw: as \(t\) approaches \(T_c\), the denominator approaches zero and the population shoots up towards infinity. Even when modeling the world's population, this model predicts an ever increasing population rate when in fact, in 1962, the world population growth rate \(N\cdot dN/dt\) peaked at 2.1% and has since been decreasing¹⁶. To avoid the singularity predicted by Von Foerster et al's model, the Russian physicist S.P. Kapitsa suggested the following approach:

\[\frac{dN}{dt}=\frac{C}{(T_c-t)^2 + \tau^2}\qquad\qquad(13)\]

The differential equation in (13) can be easily integrated, resulting in:

\[N(t)=\frac{C}{\tau}arccot\left(\frac{T_c-t}{\tau}\right)\qquad\qquad(14)\]

With (14), Kapitsa fitted his model to the world population data available at that time and found that \(C=(186\pm 1)\cdot 10^9\), \(T_c=2007\pm 1\), and \(\tau=42\pm1\). The parameter \(\tau\) represents, according to Kapitsa, the average lifespan of a human generation.

Models with variable carrying capacity

The logistic models covered so far implicitly assume that the carrying capacity is a constant, presumably dependent on the environment's capacity to sustain a species population. While this might be plausible with different populations of plants and animal species, the human species sets itself apart from the others in one important aspect: its unprecedented and unmatched capacity to significantly alter the environment (and indeed the planet) to suit its purposes.

As a first approach, we could consider that the carrying capacity of the earth is gradually expanded by a growing human population, an idea grounded on the observation that a greater number of people imply greater food productivity and eventually, a greater number of inventions to amplify the productivity of one individual. We would have the following model for population growth:

\[\frac{dN}{dt}=rN\left(1-\frac{N}{K}\right)\]	(15)
\[\frac{dK}{dt}=\gamma K N\]

In the first equation of (15) we can easily spot our old friend, the Verhulst logistic model. This time however, as the second equation of (15) implies a growing capacity always one step ahead of the population, it doesn't take much to convince ourselves that this model implies a runaway an unbounded population. In fact, it can be shown that prior to the critical time point \(t_c\), this model behaves very much like the hyperbolic model described by equation (12)¹⁷. Still, we intuitively know that there must be some physical limits for the earth's population and that human technological innovation cannot push the carrying capacity indefinitely towards infinity. For example, we have already pointed out that since 1962, the relative population growth rate \(1/N\cdot dN/dt\) has started to gradually to decrease. This occurred in the same time as the so-called Green Revolution, which brought about significant increases in food productivity and crop yields¹⁸. However, the Green Revolution did not bring about a surge in the relative population growth rate and quite the opposite occurred - the world population's growth seems to be slowing down.

At any rate, there is no reason to suppose that the carrying capacity remains constant throughout human history, and surely, there must be positive and negative influences on the carrying capacity. On the one hand, the idea that each person contributes to the growth in human carrying capacity is not to be discarded entirely. On the other hand, the contribution of each additional person depends on the resources available, and these resources must be shared among more people as the population increases. What is being described here can be summarized in the following population-growth model proposed by Cohen¹⁹, which he called the Condorcet-Mill model, after the British philosopher John Stuart Mill (1806-1873) who is cited as foreseeing that a stationary population is both inevitable and desirable.

\[\frac{dN}{dt}=rN\left(1-\frac{N}{K}\right)\]	(16)
\[\frac{dK}{dt}=\frac{L}{N}\frac{dN}{dt}\]

Introducing the FME package and some comments on the world population data we'll be using

We have finally come to the part of this blog post in which we will try to fit the world population data to the various population growth models covered. The data was obtained fro the ourworldindata.org site²⁰ and cosists of a tabular csv file with two columns, which i have renamed "time" (the year) and "pop" (population). This data covers a wide time-span, from 10,000 BC until 2015 and to be perfectly honest, I cannot vouch for the validity of the population counts in the remote past nor do I have any idea as to how they obtained census data for this time period. I'll just assume that the population data for the pre-historic and ancient time periods is fairly accurate. At any rate, this large time-scale population data allows us to compare the family of models based on logistic population growth to the singular approach by Kapitsa.

For the model fitting process, we will use an R-package called FME (FME stands for Flexible Modelling Environment). FME introduces added functionalities over the simecol package to aid in the process of system dynamics model building:

• FME's modFit incorporates additional non-linear optimization methods for parameter fitting over simecol's fitOdeModel: "Marq" for the Levenberg-Marquardt algorithm, which is a least squares method, "Pseudo" for the pseudo-random search algorithm, "Newton" for a Newton-type algorithm, and "bobyqua" for a derivative-free optimization using quadratic approximations.
• Global sensitivity analysis, which consists of assessing the change of certain model output variables according to changes in parameters over a large area of the parameter space, is done in FME via the sensRange function.
• The sensFunc function in the FME package provides a way to define sensitivity functions in order to perform a Local sensitivity analysis, which consists of studying the effects of one parameter when it varies over a very small neighboring region near its nominal value.
• The sensitivity functions are also used by FME to estimate the collinearity index of all possible subsets of the parameters. This has to do with the identifiability of the parameters in the model, a concept well worth delving into.

The nonlinear optimization methods used to fit ODE models like the ones I have been dealing with in this series of posts are very different from methods like the ordinary least squares used in linear regression, for example. The ordinary least squares method, on account of being linear on the model parameters, will always come up with a unique estimation of these parameters. The least squares estimation in linear regression is unique because there's only one global minimum in the parameter search space. In the case of nonlinear optimization problems like the ones that we deal with in ODE models, however, the parameter search space is much more complex and there may be many many local minima. This situation substantially complicates ODE model fitting, particularly when the model has several parameters.

When one is fitting ODE models, different initial "guesses" for the parameter values will result in entirely different model parameters and indeed, in an entirely behavior of the model's variables. To determine just how much do small changes in the initial parameter values affect the model's variables is the reason why we perform sensitivity analysis on the models. It may occur that some parameters are linearly or almost linearly dependent on others, and just like in the linear regression case, this multicollinearity negatively affects the identifiability of the model's parameters. In particular, the colinearity index given by the FME collin function gives a measure of the parameter set's identifiability based on the data we use to fit the model.

The collin function provides the collinearity index for all possible subsets of a model's parameters. The larger the collinearity index for any given set of parameters, the less identifiable those parameters are. As a rule of thumb, a parameter set with a collinearity index less than 20 is identifiable²¹.

Like sands through the hourglass, so are the days of our lives

We now begin the process of fitting various population models to the world population data described above. We will begin with the classic (Verhulst) Logistic growth model with 2 parameters:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31

library("simecol") library("FME") #Read census data census <- read.table("world_pop.csv",header=TRUE,sep=",") #Verhulst logistic growth model verhulst_model <- new("odeModel", main=function(time,init,parms, ...) { x <- init p <- parms dpop <- x["pop"]*p["r"]*(1-(x["pop"]/p["K"])) list(dpop)}, init=c(pop=2431214), parms=c(r=1e-02,K=1e10), times=seq(from=-10000,to=2015,by=1), solver="lsoda") #cost function fCost <- function(p) { parms(verhulst_model) <- p out <- out(sim(verhulst_model)) return(modCost(out, census, weight="std")) } #model fitting result <- modFit(fCost, parms(verhulst_model), lower=c(r=0,K=2431214),upper=c(r=1,K=1e+12)) (parms(verhulst_model) <- result$par) cat("SC : ",ssqOdeModel(NULL,verhulst_model,census$time,census[2]),"\n") #identifiability analysis sF <- sensFun(func = fCost, parms = parms(verhulst_model)) (colin <- collin(sF))

The model's output - parameter estimations, sum of squares measure of goodness of fit as calculated in line 28 of the above code, and the collinearity index of the parameter set (produced by lines 30-31 of the code) - is given below, along with the plot for the real population curve (in red) and the fitted model (in cyan). I have not included the source code for generating these plots, but could do so if any reader is interested.

           r            K 
6.065198e-04 9.999754e+11 
SC :  92.33517
  r K N collinearity
1 1 1 2          7.9

Let's interpret the results above. First of all, with a collinearity index of 7.9, the model's parameters are identifiable, so that's not an issue. An inspection of the r and K parameters reveal that the Verhulst Logistic growth model predicts a maximum world population of almost \(1\cdot 10^{12}\) or one trillion inhabitants. Considering we're currently only at 8 billion inhabitans, we would have a long way to go before we reach that saturation point. One look at the graph reveals that this model is simply not adequate. The model's curve seems to bulge too much above the real population curve, indicating that this model predicts too much population growth at remote times in the past when we know that the real growth was barely noticeable. All the while, the model predicts a population of 3-4 billion in 2015 when we know the real population level to be over 7.3 bilion this year. So we have two major flaws with this model: 1) it predicts the modern, explosive, population growth to begin much sooner than it really did and 2) in modern times (1800-present), the model's predicted growth is simply not fast enough compared to the real data. These flaws are confirmed by the model's high sum of squares value of 92.34.

One of the limitations of the classic logistic growth model is the inability to "tweak" the inflection point. In this case, with the model's K value estimated at 1 trillion, the inflection point would be half of that, or 500 billion - too far into the future. We know for a fact that the increase in the world population growth rate started to slow down during the 1960's, so what we would need is a logistic growth model in which the inflection point can be closer to the maximum population levels. We will next attempt to fit the generalized logistic growth model, which would, in theory, afford us the flexibility in tweaking with the inflection point.

The problem with the generalized logistic growth model is the large number of parameters (5) with just one differential equation. This results in a very complex parameter space with potential problems in parameter identifiability. To explore this situation, I decided to sample one hundred initial parameter guess values (N=100) using the uniform probability distribution with a given range for each parameter and storing these initial values in a data-frame (see code snippet below), from where I proceeded to use them as initial values for the model fitting procedure.

1 2 3 4 5

parameters <- data.frame( r=runif(N,min=0,max=0.1), K=runif(N,min=7349472099,max=1e+11), alpha=runif(N,min=0.1,max=1.2), beta=runif(N,min=0.1,max=2), gamma=sample((1:10),size=N,replace=TRUE) )

I must comment on the reason why I sampled only integer values for the gamma (\(\gamma\)) parameter. Taking a look at model equation (7), we see that if the expression within parentheses on the right hand side that's raised to the gamma exponent ever becomes negative, this would raise numerical error issues because we'd be raising a negative base to a fractional exponent. While there are possibly ways around this issue, it was simpler for me to just use integer values for gamma. Now when fitting the generalized logitic model using these 100 different initial parameter sets, I obtained different parameter estimations with different sum of squares goodness of fit values, which I wrote to a csv file with the following code:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48

library("simecol") library("FME") #Read census data census <- read.table("world_pop.csv",header=TRUE,sep=",") #Generalized logistic growth model generalized_model <- new("odeModel", main=function(time,init,parms, ...) { x <- init p <- parms dpop <- x^p["alpha"]*p["r"]*(1-(x/p["K"])^p["beta"])^p["gamma"] list(dpop)}, init=c(pop=2431214), parms=c(r=1e-08,K=7349472099,alpha=1,beta=1,gamma=4), times=seq(from=-10000.0,to=2015.0,by=1), solver="lsoda") #cost function fCost <- function(p) { parms(generalized_model) <- p out <- out(sim(generalized_model)) return(modCost(out, census, weight="std")) } parameters <- read.table("initial_parms.csv",sep="\t",header=TRUE) parameters <- as.matrix(parameters) for (i in 1:N) { #get initial parameter estimations and fit the model parms(generalized_model) <- parameters[i,] options(show.error.messages = FALSE) try( fitted_model <- modFit(fCost, parms(generalized_model) , lower=c(r=0,K=7349472099,alpha=0.1,beta=0.1,gamma=1), upper=c(r=0.1,K=1e+11,alpha=1.2,beta=2,gamma=20) ) ) options(show.error.messages = TRUE) parms(generalized_model) <- fitted_model$par r <- as.numeric(fitted_model$par["r"]) K <- as.numeric(fitted_model$par["K"]) a <- as.numeric(fitted_model$par["alpha"]) b <- as.numeric(fitted_model$par["beta"]) g <- as.numeric(fitted_model$par["gamma"]) Ninf <- as.numeric((1+b*g/a)^(-1/b)*K) output <- out(sim(generalized_model)) SC <- ssqOdeModel(NULL,generalized_model,census$time,census[2]) results <- rbind(results, data.frame( r=r,K=K,alpha=a,beta=b,gamma=g, Ninf=Ninf,SC=SC) ) cat(i," ",SC,"\n") }

Choosing the parameter set with the least sum-of-squares (SC) value, I ran the simulation of the model and obtained the following graph and the following collinearity indexes:

           r           K   alpha     beta    gamma       N_inf       SC
1.836653e-05 99966063867 1.19538 1.996857 1.009536 60943709007 76.60192

   r K alpha beta gamma N collinearity
1  1 1     0    0     0 2          5.8
2  1 0     1    0     0 2        214.0
3  1 0     0    1     0 2          6.2
4  1 0     0    0     1 2          5.8
5  0 1     1    0     0 2          6.0
6  0 1     0    1     0 2         81.3
7  0 1     0    0     1 2       6442.9
8  0 0     1    1     0 2          6.3
9  0 0     1    0     1 2          5.9
10 0 0     0    1     1 2         81.0
11 1 1     1    0     0 3        403.1
12 1 1     0    1     0 3        131.9
13 1 1     0    0     1 3       9122.3
14 1 0     1    1     0 3        425.7
15 1 0     1    0     1 3        403.4
16 1 0     0    1     1 3        133.0
17 0 1     1    1     0 3        133.7
18 0 1     1    0     1 3       9018.9
19 0 1     0    1     1 3       6751.7
20 0 0     1    1     1 3        134.9
21 1 1     1    1     0 4        608.5
22 1 1     1    0     1 4      14594.6
23 1 1     0    1     1 4      11876.4
24 1 0     1    1     1 4        604.7
25 0 1     1    1     1 4      11704.0
26 1 1     1    1     1 5      14777.0

R Puzzle No. 1 - Five Sailors, a Monkey and some Coconuts

In a brief diversion from the simecol dynamic system models, which I promise to resume in the future, I'll start another series of entertaining posts about solving mathematical/logic puzzles using R, of which this post you're reading is the first. This time around, we'll take a crack at the Problem of the Five Sailors, One Monkey and a Pile of Coconuts. Read on to find out what this is all about.

The Puzzle

This problem first appeared in the Saturday Evening Post in the 1920's and was later popularized by Martin Gardner. It goes as follows:

Five sailors arrive at a deserted island that has only coconuts and one monkey. After laboriously collecting all the coconuts into one big pile, the sailors all agree to divide up the coconuts into equal shares the next morning. However, in the middle of the night a sailor wakes up and, not trusting the others, decides to give the monkey one coconut (to buy his silence?) and takes away his 1/5 share of the coconuts, leaving the rest in the pile. The other 4 sailors wake up at different time intervals and each does the same thing: give one coconut to the monkey and stash away his 1/5 share of the coconuts that remain in the pile. In the morning they divide what is left of the pile into equal shares and there is still one coconut left for the monkey. How many coconuts were in the pile the day before?

Readers in the 1920's did not have the benefit of Google or the Internet so after racking their brains for a while, many gave up in frustration and soon, the premises of the newspaper were flooded by an angry mob bearing pitchforks and torches, demanding to know the answer¹. Now, before you fire up the Google search engine for the answer to this conundrum, let's first try a...

Brute-force solution

By brute-force, I really do mean brute as opposed to smart. It would seem that most people think mathematicians arrive at sophisticated answers directly in one step. Not me - I'm a dumb one so I usually try to explore the problem situation and try a little brute-force searching first. I don't actually have to do the computations myself, that's what R is for. Only after I understand the problem more thoroughly do I attempt to tackle it using more formal mathematical artifices.²

Our brute-force strategy will be to try out successive integer values for the size of the coconut pile to search for one where we will be able to deal out the coconut pile to the 5 sailors and the monkey in the conditions described by the problem. For each candidate integer value, we will simulate each of the five sailors taking their shares and giving one to the monkey, and in the end, dividing the pile into five equal shares plus one for the monkey. The key here is that the divisions in each step must be exact, since we can't deal with fractions of coconuts, otherwise the problem would be too easy, wouldn't it? Here is my brute-force solution script:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37

#------------------------------------------------------------------ # Five Sailors, a monkey and a pile of coconuts # José Romero - 05/09/2016 #------------------------------------------------------------------ n_sailors <- 5 n_iterations <- 6 max_bound <- 1000000 verify <- function(pile, verbose=FALSE) { solution_DF <- data.frame() for (i in 1:n_iterations) { sailors_share <- pile%/%n_sailors for_monkey <- pile%%n_sailors if (for_monkey == 1) { solution_DF <- rbind(solution_DF, data.frame("Pile"=pile, "Sailor.s.share"=sailors_share, "For.monkey"=for_monkey)) pile <- pile - (sailors_share+for_monkey) } } return(solution_DF) } solve_monkey_coconuts <- function(max_bound) { for (i in 1:max_bound) { if (nrow(verify(i))==6) return(i) } return(0) } solution <- solve_monkey_coconuts(max_bound) if (solution>0) { cat("Number of coconuts in pile: ",solution,"\n") print(verify(solution)) } else sprintf("No solution in [1,%d]",max_bound)

     Number of coconuts in pile:  15621 
        Pile Sailor.s.share For.monkey
     1 15621           3124          1
     2 12496           2499          1
     3  9996           1999          1
     4  7996           1599          1
     5  6396           1279          1
     6  5116           1023          1

There it is - the pile originally had 15621 coconuts and there were 1023 coconuts for each sailor (and one for the monkey) at daybreak. The one thing that strikes us when we run the script above is how agonizingly slow it is. Since we initially have no idea how large the solution is, we brute-search in the first one million positive integers using the for loop in lines 27-29³.

To check if an integer candidate is a solution to the problem, we simulate 6 partitionings (one for each of the five sailors and the final division of the coconuts done in the morning). Each partitioning consists of dividing the coconut pile into five equal parts and verifying that the remainder (integer modulo division) is one, which is the coconut for the monkey. If at some point in the process we do not get a remainder of 1, we will be exiting the verfication process, knowing that we could not perform all 6 iterations of it. This is what the code in lines 10-24 does.

Finally, note the output below the R script, indicating how many coconuts remain in the pile, the share for each sailor and the share for the monkey at each iteration. Could we have somehow guessed this would be the answer? Well, we could have ventured that since we were going to divide the pile by 5 six times, the answer would need to be divisible by \(5^6=15626\), this being our approximate answer. Of course, we would have known this could not be the answer to the puzzle because we had not taken into consideration the part about giving one coconut to the monkey at each step, which complicates things a great deal.

Diophantine equations to the rescue

Now that we know the answer to the problem, let's see if we can deduce it mathematically. Let \(C_0\) be the number of coconuts the five sailors originally collected. When the first sailor comes along in the middle of the night, he divides this pile into five equal parts and gives the remaining coconut to the monkey. So if \(S_1\) is his one fifth share of the coconuts, then the following equation holds: \(C_0=5S_1+1\). So when the next sailor comes along, there will be \(C_1=4S_1\) coconuts in the pile, the previous sailor having taken his share, plus one for the monkey. But if he also divides that pile into 5 equal parts and throws the remaining coconut to the monkey, then \(C_1=5S_2+1\) also holds, and so on. Writing down all six of these equations describing the division of coconuts at each stage:

\[\begin{align*} C_0&=1+5S_1\\ C_1=4S_1&=5S_2+1\\ C_2=4S_2&=5S_3+1\qquad\qquad(1)\\ C_3=4S_3&=5S_4+1\\ C_4=4S_4&=5S_5+1\\ C_5=4S_5&=5S_6+1 \end{align*} \]

Oh great, six equations and twelve unknows?! Please bear with me a little longer. With the equations in (1), we can try to express five of the \(S_i\)'s in function of their sucessors \(S_{i+1}\):

\[\begin{align*} S_1&=\frac{1+5S_2}{4}\\ S_2&=\frac{1+5S_3}{4}\\ S_3&=\frac{1+5S_4}{4}\qquad\qquad(2)\\ S_4&=\frac{1+5S_5}{4}\\ S_5&=\frac{1+5S_6}{4} \end{align*}\]

Some backward substitution, starting with the last two equations in (2) and working on up to the top equation there allows us to express \(S_1\) in terms of \(S_6\) and substituting this in the first equation of (1), we get:

\[ C_0=\frac{5}{4}\left(\frac{5}{4}\left(\frac{5}{4}\left(\frac{5}{4}\left(\frac{5}{4}\left(5S_6+1\right)+1\right)+1\right)+1\right)+1\right)+1\qquad(3) \]

Working with this last equation:

\[ C_0=\frac{5^6}{4^5}S_6+\sum_{i=0}^5 \left(\frac{5}{4}\right)^i\qquad(4) \]

In the right hand side of (4) we have a geometric series so we apply the pertinent formula:

\[ C_0=\frac{5^6}{4^5}S_6+\frac{\left(\tfrac{5}{4}\right)^6-1}{\tfrac{5}{4}-1}\qquad(5) \]

This means that:

\[\begin{align*} C_0&=\frac{5^6}{4^5}-4+\frac{5^6}{4^5}S_6\quad\rightarrow\\ 4^5C_0&=5^6-4^6+5^6S_6\quad\rightarrow\quad&(5a)\\ 1024C_0&-15625S_6=11529\qquad\quad&(5b) \end{align*}\]

Equation (5b), which looks just like an innocent linear equation on two variables, is called a linear diophantine equation. If \(C_0\) and \(S_6\) were allowed to assume real values, we simply would have infinite solutions, but there's a catch here: we're looking for integer values of \(C_0\) and \(S_6\) as solutions to (5b). This is typical of diophantine problems, which usually involve less equations than unknowns but in which integer solutions are required.

What comes next is a little sleight of hand. Due to the symmetry in equation (5a), we can see that if \(C_0\) and \(S_6\) are one solution to the equation, so are \(C_0'=C_0+5^6\) and \(S_6'=S_6+4^5\). If we substite \(C_0'\) for \(C_0\) and \(S_6'\) for \(S_6\) in equation (5a) and notice that the \(4^55^6\) terms cancel out on both sides of the equality, the resulting equation is essentially the same:

\[ 4^5(C_0+5^6)=5^6(S_6+4^5)+5^6-4^6 \]

What this all amounts to is that we can find other solutions by adding some constants to each of the variables in a known solution pair. But, how do we find one solution pair? The trick is to add \(4^6\) to both sides of (5a) and to factor out both sides:

\[ 4^5(C_0+4)=5^6(S_6+1)\qquad(6) \]

Now this last equation is easy to solve. We just need to make both sides of (6) equal to zero and obtain \(C_0=-4\) and \(S_6=-1\). But the problem is, how can we have negative coconuts?⁴ This is not really a problem because if we add \(5^6\) to \(C_0\) and \(4^5\) to \(S_6\), we know we can get another solution pair just as mathematically valid, but more sensible, given the problem. We finally get \(C_0=15625-4=15621\) and \(S_6=1024-1=1023\), which concur with the solutions we obtained by brute-force.

Taking a crack at diophantine equations with R

We're not about to give up on trying to find a more efficient and possibly vectorized way to solve this problem using R. I looked for an R package to solve diophantine equations and found that the "elliptic" package has a function called "congruence" which deals with linear diophantine equations of the type \(mx+by=l\).

To be honest, I could never get this function to give me the solution to equation (5b). Poring through the package's vignette, I encountered some frightening number theoretic stuff, so given my ignorance and aversion to number theory⁵, I gave up on using this package but before doing so, I inspected the R code for the congruence function:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32

function (a, l = 1) { l <- as.integer(l) m <- as.integer(a[1]) n <- as.integer(a[2]) zero <- as.integer(0) one <- as.integer(1) if (m == zero & n == one) { return(NULL) } if (m == one & n == zero) { return(c(NA, 1)) } if (m == 1) { return(rbind(c(a), c(1, n + l), c(0, l))) } if (n == 1) { return(rbind(c(a), c(m - l, 1), c(l, 0))) } q1 <- which((+l + (1:m) * n)%%m == zero) if (!any(q1)) { q1 <- NA } q2 <- which((-l + (1:n) * m)%%n == zero) if (!any(q2)) { q2 <- NA } out <- rbind(a, cbind(q1, q2)) rownames(out) <- NULL colnames(out) <- NULL return(out) }

Hmmm, the code above doesn't look so inscrutable. I mean, it uses the expected integer modulo operator in R (the %%'s in the code). Maybe if I knew a little bit more about diophantine equations I could reimplement this function from scratch and get it to work. So I looked up the Wikipedia article on diophantine equations and found the following theorem:

The linear diophantine equation \(ax+by=c\), where \(a\),\(b\),\(c\),\(x\) and \(y\) are integers, has a solution if and only if \(c\) is a multiple of the greatest common divisor of \(a\) and \(b\). Moreover, if \((x, y)\) is a solution, then the other solutions have the form \((x + kv,y − ku)\), where \(k\) is an arbitrary integer, and \(u\) and \(v\) are the quotients of \(a\) and \(b\) (respectively) by the greatest common divisor of \(a\) and \(b\).

This theorem is really handy for two reasons: (1) it gives us the conditions for the existence of solutions to the diophantine equation and (2) it gives us a mechanism for generating infinitely more solutions if we are able to find one solution. This last point, finding one solution, is what I would have to tackle in my implementation of the code.

I decided that I would look for a first solution by separating each variable to one side of the \(ax+by=c\) equation and getting \(ax=c-by\) or \(by=c-ax\). Then, taking \(ax=c-by\) for example, I would have my code try out a range of values for \(y\) to verify when the expression \(c-by\) is divisible by \(a\). Intuitively, I realized that the range of values to try for \(y\) would be all integers between 0 and \(a\) (when using \(by=c-ax\), I would have to try all integers between 0 and \(b\) for \(x\)), so it's not like I would have to search through an infinitely large solution space.

Having found values of \(x\) or \(y\), which are guaranteed to exist if the solution existence conditions in the Wikipedia theorem are met, I would substitute this value in the equation to obtain the other value. Finally, I would have to implement a function for finding the greatest common divisor. All in all, I came up with the following code:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30

gcd <- function(m,n) { #greatest common divisor function m <- as.integer(abs(m)) n <- as.integer(abs(n)) divm <- (1:m)[which(m%%(1:m)==0)] divn <- (1:n)[which(n%%(1:n)==0)] return(max(divm[which(divm%in%divn)])) } #gcd diophantine <- function (m,n,l) { #solves the linear diophantine equation: mx+ny=l l <- as.integer(l) m <- as.integer(m) n <- as.integer(n) gcd <- gcd(m,n) if (l%%gcd!=0) return(NA) #no solution exists #otherwise obtain one solution (solx,soly) if (abs(m)>abs(n)) { #use multiples of the lower magnitude n mult <- (0:n)[which((l-m*(0:n))%%n==0)] solx <- mult[1]; soly <- (l-m*solx)%/%n } else { #use multiples of the lower magnitude m mult <- (0:m)[which((l-n*(0:m))%%m==0)] soly <- mult[1]; solx <- (l-n*soly)%/%m } return(rbind(c(solx,soly),c(n%/%gcd,-m%/%gcd))) } #diophantine diophantine(1024,-15625,11529)

##        [,1]  [,2]
## [1,]  15621  1023
## [2,] -15625 -1024

As you can see, it works like a charm. This is a code I can be proud of (look ma, no for/while loops!). My function not only gives me the solution in the first row of the results matrix, but gives me the coefficients whose multiples I would have to add to \(x\) and \(y\) to get other solution pairs. I haven't really proven that it works for all possible values of m, n and l, but I tried a few example problems which I solved by hand and it seems to be ok. Needless to say, I invite my readers to put this code under close scrutiny (you know what they say about many pairs of eyes seeing more than one). Perhaps some of you will find errors in the code, or a test-case which "breaks" the code. Also, I invite my readers to come up with an R function for solving general 3 or more variable linear diophantine equations. For further discussion or questions of this, I invite you to participate in the comments for this post.

Notes

1. I guess I got carried away with that pitchfork and torch bearing angry mob bit, but I do wish we Venezuelans would organize into such a mob to evict Maduro from power. Ironically, it seems like in the days when there were no social networks or smartphones, people somehow organized into action more effectively - for example during the 1958 ousting of the Perez-Jiménez dictatorship. Now, we somehow expect the Internet to "do things for us", while we stay at home.
2. When teaching math, it is never a good idea to just go ahead and show students the "clean" mathematical solution. Teaching math is all about modeling the problem-solving process, which in real life, is usually a messy, trial-and-error roundabout affair. When we have become experts in making mistakes and following dead-end leads, we will have gained sufficient familiarity with the problem and only then can we hope to come up with a clean solution.
3. R is not particularly fast when running for/while loops. What R programmers always strive to do is to vectorize their code. Vectorization consists of avoiding for/while loops all together and working with vectors and functions of vectors. In the next script, we will be using vectorized code. So you could say that in R, using for or while loops in brute-search is the most brutish way to go about it.
4. I wonder what negative coconuts taste like...
5. My apologies to number-theorists (they would cringe at a blog like this anyways), but number theory seems to me like the most useless branch of mathematics. I mean, of what practical use is Fermat's Last Theorem or Goldbach's Conjecture? Sure, probably the most difficult and even unresolved or unproven math problems are buried deep within number theory. Sure, solving these problems would be an instant tiket to mathematical fame and fortune - the stuff of legends. As for me, I prefer to steer clear of things for which I just don't have the talent, so that I won't suffer the same fate as Uncle Petros. Uncle Petros, by the way, is the main character in an extremely entertaining novel about mathematics by Apostolos Doxiadis, called Uncle Petros and Goldbach's Conjecture

Bibliography

• Diophantine equation. (2016, July 26). In Wikipedia, The Free Encyclopedia. Retrieved September 4, 2016, from https://en.wikipedia.org/wiki/Diophantine_equation
• HANKIN, R. K. S. (February 2006). Introducing Elliptic, an R package for Elliptic and Modular Functions. Journal of Statistical Software, Volume 15, Issue 7.
• PADILLA, T. and McPARTLAN, P. [Numberphile] (2015, April 29). Monkeys and Coconuts - Numberphile [Video file]. Retrieved September 4, 2016, from https://www.youtube.com/watch?v=U9qU20VmvaU
  .
• R DEVELOPMENT CORE TEAM (2009). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. ISBN 3-900051-07-0. http://www.R-project.org.

If you found this post interesting or useful, please share it on Google+, Facebook or Twitter so that others may find it too.