# Logistic growth model¶

In this example, we use the logistic growth model in order to show how to define a function which has a vector input and a field output. We use the OpenTURNSPythonPointToFieldFunction class to define the derived class and its methods.

## Introduction¶

The logistic growth model is the differential equation: for any , with the initial condition: where

• and are two real parameters,

• is the size of the population at time ,

• is the initial time,

• is the initial population at time ,

• is the final time.

The parameter sets the growth rate of the population. The parameter acts as a competition parameter which limits the size of the population by increasing the competition between its members.

In , the author uses this model to simulate the growth of the U.S. population. To do this, the author uses the U.S. census data from 1790 to 1910. For this time interval, R. Pearl and L. Reed  computed the following values of the parameters: Our goal is to use the logistic growth model in order to simulate the solution for a larger time interval, from 1790 to 2000: Then we can compare the predictions of this model with the real evolution of the U.S. population.

We can prove that, if , then the limit population is: In 1790, the U.S. population was 3.9 Millions inhabitants: We can prove that the exact solution of the ordinary differential equation is: for any .

We want to see the solution of the ordinary differential equation when uncertainties are taken into account in the parameters:

• the initial U.S. population ,

• the parameters and .

Indeed, Pearl and Reed  estimated the parameters and using the U.S. census data from 1790 to 1910 while we have the data up to 2000. Moreover, the method used by Pearl and Reed to estimate the parameters could be improved; they only used 3 dates to estimate the parameters instead of using least squares, for example. Finally, Pearl and Reed did not provide confidence intervals for the parameters and .

## Normalizing the data¶

The order of magnitude of the population is . This is why we consider the normalized population in millions: for any .

Let be the initial population: ## Uncertainties¶

Uncertainty can be accounted for by turning , and into independent random variables , and with Gaussian distributions. From this point onward, , and respectively denote , and .

Variable

Distribution gaussian, mean , coefficient of variation 10% gaussian, mean , coefficient of variation 30% gaussian, mean , coefficient of variation 30%

No particular probabilistic method was used to set these distributions. An improvement would be to use calibration methods to get a better quantification of these distributions. An improvement would be to use calibration methods to get a better quantification of these distributions.

## Define the model¶

:

from __future__ import print_function
import openturns as ot
from numpy import linspace, exp, maximum


The data is based on 22 dates from 1790 to 2000.

:

ustime=list(range(1790,2001,10))
uspop=[3.9,5.3,7.2,9.6,13.,17.,23.,31.,39., 50.,62.,76.,92.,106.,123.,132.,151.,179., 203.,221.,250.,281.]


Creation of the input distribution.

:

y0 = 3.9e6
a = 0.03134
b = 1.5887e-10
distY0 = ot.Normal(y0, 0.1 * y0)
distA  = ot.Normal(a, 0.3 * a)
distB  = ot.Normal(b, 0.3 * b)
distX = ot.ComposedDistribution([distY0, distA, distB])


The model.

:

class Popu(ot.OpenTURNSPythonPointToFieldFunction):

def __init__(self, t0 = 1790.0, tfinal = 2000.0, nt = 1000):
grid = ot.RegularGrid(t0, (tfinal - t0) / (nt - 1), nt)
super(Popu, self).__init__(3, grid, 1)
self.setInputDescription(['y0', 'a', 'b'])
self.setOutputDescription(['N'])
self.ticks_ = [t for t in grid.getVertices()]
self.phi_ = ot.SymbolicFunction(['t', 'y', 'a', 'b'], ['a*y - b*y^2'])

def _exec(self, X):
y0 = X
a  = X
b  = X
phi_ab = ot.ParametricFunction(self.phi_, [2, 3], [a, b])
phi_t = ot.ParametricFunction(phi_ab, , [0.0])
solver = ot.RungeKutta(phi_t)
initialState = [y0]
values = solver.solve(initialState, self.ticks_)
return values * [1.0e-6]

F = Popu(1790.0, 2000.0, 1000)
popu = ot.PointToFieldFunction(F)


## Generate a sample from the model¶

Sample from the model

:

size = 10
inputSample = distX.getSample(size)
outputSample = popu(inputSample)

:

ot.ResourceMap.SetAsUnsignedInteger('Drawable-DefaultPalettePhase', size)


Draw some curves

:

graph = outputSample.drawMarginal(0)
graph.setTitle('US population')
graph.setXTitle(r'$t$ (years)')
graph.setYTitle(r'$N$ (millions)')
cloud = ot.Cloud(ustime, uspop)
cloud.setPointStyle('circle')
cloud.setLegend('Data')

: 