EXPLOOT project
Predicting the capacity in a logistic model

http://we.vub.ac.be/exploot/website

I. Cnop
icnop@vub.ac.be

Dotcom crashes: can we predict the capacity in a system?

Introduction

At the end of the nineties , nobody would have believed that the portable phone (GSM)  exponential expansion could level off.  
Already in the Fall 2000, shares of Nokia and Ericsson saw their equities decline. 
After its 1Q 2001, Siemens reported a loss exceeding 300 M  DEM (about 150 M Euro) in its GSM production & sales division. 
Price competition between manufacturers is fierce and investments (including underwritten debt) in third generation content delivery technology is still not paying off.
The reason is that unlimited expansion has changed to a replacement economy in the original markets.

Is this maximum capacity predictable? The answer is:  YES  !

The market for portable (mobile cell phone or GSM) use is a typical example of the logistic equation. 
Its success can be explained by the fact that users tend to incite non-users to join in, thus creating a nonlinear multiplier effect. 
The logistic equation can be solved exactly (up to rescaling of the time-axis) if the capacity of the system is known.

The method below was originally developed to answer another problem: to recover the asymptotic capacity in population estimates for several countries and continents published in graphical form in the seventies by UNESCO, relying on census figures from 1950 to 1970. 
This problem was complicated by the fact that only the graphical output was available and original census figures where not.

Solution of the logistic model

General solution

We assume increase per time unit (say one month) is proportinal to the number of portable users multiplied by the number of non users:

sol  = 
      DSolve[ n'[t]  == r n[t] ( 1  -  n[t]/K), n[t], 
        t]

n[t] /. First[%] // TraditionalForm

Initial condition and parametrisation

Let us start with 1 k potential customers and 10 users to begin with. Select some constant r which indicates that over 1 month, the probability that an encounter will lead to 30% chance that the non-user becomes a user himself

K = 1000; r =  .3

Clear[n, sol, t];
sol  = 
      NDSolve[{ n'[t]  == r  n[t] ( 1  -  n[t]/K), 
          n[0] == 10}, n[t], {t, 0, 50}]

graph1 = Plot[n[t] /. First[%], {t, 0, 45}]

Classical approximations

It is not possible to approximate such graphic by a linear, quadratic, ... polynomial  function

Start with points selected on or near this graph

pts = {{0.407056, 17.4921}, {2.65786, 38.6217}, {4.55385, 
          106.795}, {7.03749, 281.014}, {8.81152, 472.377}, {10.5856, 
          672.111}, {12.5924, 863.474}}

Fit[%, {1, x, (x^2)}, x]

Show[{graph1, Plot[%, {x, 0, 45}]}, 
      PlotRange -> {0, 2500}]

Exercise: repeat this with a third or fourth degree polynomial.

What conclusions can you make?

None of these polynomial functions has an horizontal asymptote indicating the asymptotic limiting value, called capacity K.

This is precisely the value we are looking for. But the search is a delicate one:

Changing K gives a completely different solution

K = 2000

Clear[n, sol];
sol  = 
      NDSolve[{ n'[t]  == r  n[t] ( 1  -  n[t]/K), 
          n[0] == 10}, n[t], {t, 0, 50}]

graph2 = 
      Plot[n[t] /. First[%], {t, 0, 45}, 
        PlotRange -> {0, 2500}]

Show[{graph1, graph2}, PlotRange -> {0, 2500}]

The difference in the initial stage is small. 
It is only after some time that differences become important, and the top graph keeps growing while the bottom one is tapering off. 
This shows it is a delicate task to predict the asymptotic value.

Can we predict the capacity  K  of a system?

It is known that the capacity is twice the value at the inflection point. 
This inflection point is hard to locate experimentally for several reasons: the graph is steepest there, and taking succesive differences may lead to large errors if time intervals are uneven. 
Moreover we would like to predict capacity before attaining large values in the system.

We produce a blank graph in which we can put some actual data:

Plot[0, {t, 0, 45},  GridLines -> Automatic, 
      PlotRange -> {0, K}]

These data can be saved automatically (by Copy & Paste behind  "pts =" followed by enter), and plotted:

pts = {{0.72839, 15.5049}, {2.35802, 35.0368}, {3.91005, 
            71.3102}, {5.30688, 132.696}, {6.1605, 185.711}, {6.5485, 
            213.613}, {7.09172, 269.419}, {7.40212, 302.902}, {7.94533, 
            369.868}, {8.41094, 428.464}, {8.79895, 473.108}, {9.10935, 
            506.591}};

ListPlot[pts, PlotRange -> {{0, 45}, {0, 1000}}]

Show[%, graph1]

t = First[Transpose[pts]];

times1[times_] := Take[times, Length[pts] - 1]

times2[times_] := Take[RotateLeft[times], Length[pts] - 1]

d = Last[Transpose[pts]];

data1[data_] :=  Take[data, Length[pts] - 1]

data2[data_] :=  Take[RotateLeft[data], Length[pts] - 1]

ListPlot[
      Transpose[{(data2[d]  +  data1[d])/2, 
          2 (data2[d]  -  data1[d])/(data2[d]  +  
                    data1[d])/(times2[t] - times1[t])}], 
      PlotRange -> {{0, 1000}, {0, 1}}]

Linear regression

Fit [
      Transpose[{(data2[d]  +  data1[d])/2, 
          2 (data2[d]  -  data1[d])/(data2[d]  +  
                    data1[d])/(times2[t] - times1[t])}], {1, x}, 
      x]

sol = 
      Solve[{a   +  b == % /. x -> 1, 
          b == % /. x -> 0}, {a, b}]

Show[{Plot[%%, {x, 0, 1000}], %%%}, 
      PlotRange -> {{0, 1000}, {0, 1}}]

If we obtain coefficients a and b for this linear regression a x + b, we obtain  r  and   K  by
  multiplying with x and equating with the right hand side of the differential equation.

a = -0.0004589037117878614`
b = 0.48302807562678013`

r =. ; K =. ; Expand[ r x (1  -  x/K)]

Solve[(Coefficient[%, (x^2)] /. r -> b) == a, 
      K]

a =. ; b =. ;

Discarding the first data

Manipulating inputs

Since division by small numbers induces bigger errors, it is sometimes more accurate to discard the first few measurements:

pts = {{3.91005, 71.3102}, {5.30688, 132.696}, {6.1605, 
          185.711}, {6.5485, 213.613}, {7.09172, 269.419}, {7.40212, 
          302.902}, {7.94533, 369.868}, {8.41094, 428.464}, {8.79895, 
          473.108}, {9.10935, 506.591}}

inputs = 
      ListPlot[pts, PlotRange -> {{0, 45}, {0, 1000}}]

t = First[Transpose[pts]];

times1[times_] := Take[times, Length[pts] - 1]

times2[times_] := Take[RotateLeft[times], Length[pts] - 1]

d = Last[Transpose[pts]];

data1[data_] :=  Take[data, Length[pts] - 1]

data2[data_] :=  Take[RotateLeft[data], Length[pts] - 1]

ListPlot[
      Transpose[{(data2[d]  +  data1[d])/2, 
          2 (data2[d]  -  data1[d])/(data2[d]  +  
                    data1[d])/(times2[t] - times1[t])}], 
      PlotRange -> {{0, 1000}, {0, 1}}]

Linear regression

Fit [
      Transpose[{(data2[d]  +  data1[d])/2, 
          2 (data2[d]  -  data1[d])/(data2[d]  +  
                    data1[d])/(times2[t] - times1[t])}], {1, x}, 
      x]

sol = 
      Solve[{a   +  b == % /. x -> 1, 
          b == % /. x -> 0}, {a, b}]

Show[{Plot[%%, {x, 0, 1000}], %%%}, 
      PlotRange -> {{0, 1000}, {0, 1}}]

If we obtain coefficients a and b for this linear regression a x + b, we obtain  r  and   K  by
  multiplying with x and equating with the right hand side of the differential equation.

a = -0.0004860318758329596`
b = 0.49286493038009366`

r = b; K =. ; solK =. ; Expand[ r x (1  -  x/K)]

solK = 
      Solve[(Coefficient[%, (x^2)]) == a, K]

a
b
r
K = K /. First[solK]

Plotting output

Clear[n, sol];
sol  = 
      NDSolve[{ n'[t]  == r  n[t] ( 1  -  n[t]/K), 
          n[0] == 10}, n[t], {t, 0, 30}]

graph1 = Plot[n[t] /. First[%], {t, 0, 25}]

Show[%, inputs]

The fit is remarkably good.

Statistical techniques

Statistical techniques can be used to discard outliers and to predict error intervals when error margins on measurements are known.
Moreover it is interesting to introduce measurement data one by one to see how the capacity is approximated as more measurements become available.

Simulations

Simulations can be produced to illustrate the use of different estimators and to refine the statistics involved.