Randomized discrete exponential and logistic growth models 

What is the definition of a randomized population model? A function F plus a probability distribution X? 

> restart; -1
 

> with(plots); -1; with(Statistics); -1; with(LinearAlgebra); -1
 

> `assign`(expo, proc (r) option operator; return proc (P) options operator, arrow; `+`(P, `*`(r, `*`(P))) end proc end proc); -1
 

> `assign`(logi, proc (r, K) option operator; return proc (P) options operator, arrow; `+`(P, `*`(r, `*`(P, `*`(`+`(1, `-`(`/`(`*`(P), `*`(K)))))))) end proc end proc); -1
 

> `assign`(iterate, proc (F, P0, t) option remember; if `<=`(t, 0) then return P0 else return F(iterate(F, P0, `+`(t, `-`(1)))) end if end proc); -1
`assign`(iterate, proc (F, P0, t) option remember; if `<=`(t, 0) then return P0 else return F(iterate(F, P0, `+`(t, `-`(1)))) end if end proc); -1
 

> `assign`(orbit, proc (F, P0, T) local j, orb; `assign`(orb, [P0]); for j from 2 to T do `assign`(orb, [op(orb), F(orb[`+`(j, `-`(1))])]) end do; return op(orb) end proc); -1
`assign`(orbit, proc (F, P0, T) local j, orb; `assign`(orb, [P0]); for j from 2 to T do `assign`(orb, [op(orb), F(orb[`+`(j, `-`(1))])]) end do; return op(orb) end proc); -1
 

> `assign`(exporand, proc (r, P0, T) local j, orb, rd; description
`assign`(exporand, proc (r, P0, T) local j, orb, rd; description
`assign`(exporand, proc (r, P0, T) local j, orb, rd; description
`assign`(exporand, proc (r, P0, T) local j, orb, rd; description
`assign`(exporand, proc (r, P0, T) local j, orb, rd; description
`assign`(exporand, proc (r, P0, T) local j, orb, rd; description
`assign`(exporand, proc (r, P0, T) local j, orb, rd; description
`assign`(exporand, proc (r, P0, T) local j, orb, rd; description
`assign`(exporand, proc (r, P0, T) local j, orb, rd; description
`assign`(exporand, proc (r, P0, T) local j, orb, rd; description
 

> `assign`(logirand, proc (r, K, P0, T, X) local j, orb, rdr, rdK; description
`assign`(logirand, proc (r, K, P0, T, X) local j, orb, rdr, rdK; description
`assign`(logirand, proc (r, K, P0, T, X) local j, orb, rdr, rdK; description
`assign`(logirand, proc (r, K, P0, T, X) local j, orb, rdr, rdK; description
`assign`(logirand, proc (r, K, P0, T, X) local j, orb, rdr, rdK; description
`assign`(logirand, proc (r, K, P0, T, X) local j, orb, rdr, rdK; description
`assign`(logirand, proc (r, K, P0, T, X) local j, orb, rdr, rdK; description
`assign`(logirand, proc (r, K, P0, T, X) local j, orb, rdr, rdK; description
`assign`(logirand, proc (r, K, P0, T, X) local j, orb, rdr, rdK; description
 

> `assign`(poprand, proc (F, X, P0, T) local j, orb, rd, G; description
`assign`(poprand, proc (F, X, P0, T) local j, orb, rd, G; description
`assign`(poprand, proc (F, X, P0, T) local j, orb, rd, G; description
`assign`(poprand, proc (F, X, P0, T) local j, orb, rd, G; description
`assign`(poprand, proc (F, X, P0, T) local j, orb, rd, G; description
`assign`(poprand, proc (F, X, P0, T) local j, orb, rd, G; description
`assign`(poprand, proc (F, X, P0, T) local j, orb, rd, G; description
`assign`(poprand, proc (F, X, P0, T) local j, orb, rd, G; description
`assign`(poprand, proc (F, X, P0, T) local j, orb, rd, G; description
`assign`(poprand, proc (F, X, P0, T) local j, orb, rd, G; description
 

> `assign`(tempplot, proc (S) local N; description
`assign`(tempplot, proc (S) local N; description
`assign`(tempplot, proc (S) local N; description
`assign`(tempplot, proc (S) local N; description
`assign`(tempplot, proc (S) local N; description
`assign`(tempplot, proc (S) local N; description
`assign`(tempplot, proc (S) local N; description
`assign`(tempplot, proc (S) local N; description
 

Exponential growth 

Suppose we have counted the size of a population over T days and obtained the population sizes shown in plot1. 

> `assign`(P0, 1); -1; `assign`(T, 14); -1
 

> `assign`(r, NormalDistribution(.4, .2)); -1
 

> `assign`(popobs, exporand(r, P0, T, X)); -1
 

> popobs; 1
 

[1, 1.218296384, 1.787718988, 2.597335872, 3.445250543, 5.260351783, 7.503081341, 11.71192258, 15.14186888, 18.79296858, 26.19041587, 38.07821000, 62.38167275, 80.74863187]
[1, 1.218296384, 1.787718988, 2.597335872, 3.445250543, 5.260351783, 7.503081341, 11.71192258, 15.14186888, 18.79296858, 26.19041587, 38.07821000, 62.38167275, 80.74863187]
[1, 1.218296384, 1.787718988, 2.597335872, 3.445250543, 5.260351783, 7.503081341, 11.71192258, 15.14186888, 18.79296858, 26.19041587, 38.07821000, 62.38167275, 80.74863187]		 (1)
 

> `assign`(P[max], ceil(max(popobs[T]))); -1
 

> `assign`(plot1, pointplot([seq(`<,>`(j, popobs[j]), j = 1 .. T)], symbol = solidcircle, color = red, symbolsize = 15)); -1
`assign`(plot1, pointplot([seq(`<,>`(j, popobs[j]), j = 1 .. T)], symbol = solidcircle, color = red, symbolsize = 15)); -1
 

> display(plot1, caption =
 

Plot_2d
 

Can these data be described reasonably by an exponential growth model? 

We compute the growth and the per capita growth rate during each of the T-1 time intervals and try to fit it with a horizontal line (the average of the PCGRs). 

> `assign`(GR, seq(`+`(popobs[`+`(j, 1)], `-`(popobs[j])), j = 1 .. `+`(T, `-`(1)))); -1
 

> `assign`(PCGR, seq(`/`(`*`(GR[j]), `*`(popobs[j])), j = 1 .. `+`(T, `-`(1)))); -1
 

> `assign`(L, LinearFit([1], popobs[1 .. `+`(T, `-`(1))], [PCGR])); 1
 

Vector[column](%id = 137675756) (2)
 

> `assign`(r, L[1]); 1
 

.407180759038461515 (3)
 

> display([plot(r, P = 0 .. P[max], color = blue, thickness = 2), pointplot([seq(`<,>`(popobs[j], PCGR[j]), j = 1 .. `+`(T, `-`(1)))], symbol = solidcircle, color = red, symbolsize = 15)], caption =
display([plot(r, P = 0 .. P[max], color = blue, thickness = 2), pointplot([seq(`<,>`(popobs[j], PCGR[j]), j = 1 .. `+`(T, `-`(1)))], symbol = solidcircle, color = red, symbolsize = 15)], caption =
display([plot(r, P = 0 .. P[max], color = blue, thickness = 2), pointplot([seq(`<,>`(popobs[j], PCGR[j]), j = 1 .. `+`(T, `-`(1)))], symbol = solidcircle, color = red, symbolsize = 15)], caption =
 

Plot_2d
 

The estimate is that the observed data can be descibed as discrete exponential growth with per capita growth rate r. 

We finally compare the estimated discrete exponential growth model (blue) with the actual observations (red). 

> `assign`(popmodel, orbit(expo(r), `+`(`*`(.9, `*`(popobs[1]))), T)); -1
 

> `assign`(plot2, pointplot([seq(`<,>`(j, popmodel[j]), j = 1 .. T)], symbol = solidcircle, color = blue, symbolsize = 15)); -1
`assign`(plot2, pointplot([seq(`<,>`(j, popmodel[j]), j = 1 .. T)], symbol = solidcircle, color = blue, symbolsize = 15)); -1
 

> display([plot1, plot2], caption =
 

Plot_2d
 

We didn't estimate the P0 in the model. Try to vary P0 for a better fit. 

Logistic growth 

Suppose we have counted the size of a population over T days and obtained the population sizes shown in plot1. 

> `assign`(P0, 1); -1; `assign`(T, 30); -1
 

> `assign`(r, NormalDistribution(.4, .2)); -1; `assign`(K, NormalDistribution(1200, 100)); -1
 

> `assign`(popobs, logirand(r, K, P0, T)); -1
 

> `assign`(P[max], ceil(max(popobs[T]))); -1
 

> `assign`(plot1, pointplot([seq(`<,>`(j, popobs[j]), j = 1 .. T)], symbol = solidcircle, color = red, symbolsize = 15)); -1
`assign`(plot1, pointplot([seq(`<,>`(j, popobs[j]), j = 1 .. T)], symbol = solidcircle, color = red, symbolsize = 15)); -1
 

> display(plot1, caption =
 

Plot_2d
 

Can these data be described reasonably by an logistic growth model? 

We compute the growth  and the per capita growth rate during each of the T-1 time intervals and try to fit it with a line (not necessarily horizontal). 

> `assign`(GR, seq(`+`(popobs[`+`(j, 1)], `-`(popobs[j])), j = 1 .. `+`(T, `-`(1)))); -1
 

> `assign`(PCGR, seq(`/`(`*`(GR[j]), `*`(popobs[j])), j = 1 .. `+`(T, `-`(1)))); -1
 

> `assign`(L, LinearFit([1, P], popobs[1 .. `+`(T, `-`(1))], [PCGR], P)); 1
 

`+`(.411168718258540733, `-`(`*`(0.329662314421960028e-3, `*`(P)))) (4)
 

> `assign`(L, unapply(L, P)); -1
 

> `assign`(K, solve(L(P) = 0, P)); 1
 

1247.242103 (5)
 

> `assign`(r, L(0)); 1
 

.4111687183 (6)
 

> display([plot(L(P), P = 0 .. max(P[max], K), color = blue), pointplot([seq(`<,>`(popobs[j], PCGR[j]), j = 1 .. `+`(T, `-`(1)))], symbol = solidcircle, color = red, symbolsize = 15)], caption =
display([plot(L(P), P = 0 .. max(P[max], K), color = blue), pointplot([seq(`<,>`(popobs[j], PCGR[j]), j = 1 .. `+`(T, `-`(1)))], symbol = solidcircle, color = red, symbolsize = 15)], caption =
display([plot(L(P), P = 0 .. max(P[max], K), color = blue), pointplot([seq(`<,>`(popobs[j], PCGR[j]), j = 1 .. `+`(T, `-`(1)))], symbol = solidcircle, color = red, symbolsize = 15)], caption =
 

Plot_2d
 

The estimate is that the observed data can be descibed as discrete logistic growth with per capita growth rate r and carrying capacity K. 

We finally compare the logistic growth model with parameters r and K  (blue) with the actual observations (red). 

> `assign`(popmodel, orbit(logi(r, K), `+`(`*`(1.0, `*`(popobs[1]))), T)); -1
 

> `assign`(plot2, pointplot([seq(`<,>`(j, popmodel[j]), j = 1 .. T)], symbol = solidcircle, color = blue, symbolsize = 15)); -1
`assign`(plot2, pointplot([seq(`<,>`(j, popmodel[j]), j = 1 .. T)], symbol = solidcircle, color = blue, symbolsize = 15)); -1
 

> display([plot1, plot2], caption =
 

Plot_2d
 

The choice of P0  in the model can be important for the fit of these two graphs. We didn't estimate the initial population. Try to modify P0 in the model by varying the factor in front pop_obs[1] in the model. 

>
 

Stochastic single population models 

Exponential growth 

> `assign`(T, 14); -1
 

> `assign`(r, NormalDistribution(.4, .2)); -1
 

We replicate the experiment 5000 times. and record the final population values. 

> `assign`(runs, [seq(exporand(r, 1, T), j = 1 .. 5000)]); -1
 

> `assign`(tops, [seq(runs[j][T], j = 1 .. nops(runs))]); -1
 

> DataSummary(tops); 1
 

[mean = 79.43496841, standarddeviation = 43.73085118, skewness = 1.658934330, kurtosis = 8.055955291, minimum = 7.914769298, maximum = 481.1628505, cumulativeweight = 5000.]
[mean = 79.43496841, standarddeviation = 43.73085118, skewness = 1.658934330, kurtosis = 8.055955291, minimum = 7.914769298, maximum = 481.1628505, cumulativeweight = 5000.]		 (7)
 

> Histogram(tops, frequencyscale = relative, averageshifted = 4); 1
 

Plot_2d
 

Guess: The logs of the final population sizes is a normal distribution 

> `assign`(logtops, [seq(ln(tops[j]), j = 1 .. nops(runs))]); -1
 

> DataSummary(logtops); 1
 

[mean = 4.238546730, standarddeviation = .5261152434, skewness = -0.7975874542e-1, kurtosis = 2.983301519, minimum = 2.068730545, maximum = 6.176205779, cumulativeweight = 5000.]
[mean = 4.238546730, standarddeviation = .5261152434, skewness = -0.7975874542e-1, kurtosis = 2.983301519, minimum = 2.068730545, maximum = 6.176205779, cumulativeweight = 5000.]		 (8)
 

> `assign`(pl1, plot(PDF(Normal(Mean(logtops), StandardDeviation(logtops)), t), t = min(logtops) .. max(logtops), color = red, thickness = 2)); -1
`assign`(pl1, plot(PDF(Normal(Mean(logtops), StandardDeviation(logtops)), t), t = min(logtops) .. max(logtops), color = red, thickness = 2)); -1
 

> `assign`(Hlog, Histogram(logtops, frequencyscale = relative, averageshifted = 4)); -1
 

> display([Hlog, pl1]); 1
 

Plot_2d
 

Here is a graph of the first runs of the model 

> `assign`(S, [seq([seq(runs[k][t], k = 1 .. `+`(`*`(`/`(1, 50), `*`(nops(runs)))))], t = 1 .. T)]); -1
 

> tempplot(S); 1
 

Plot_2d
 

Logistic growth 

> `assign`(P0, 1); -1; `assign`(T, 30); -1
 

> `assign`(r, NormalDistribution(.4, .2)); -1; `assign`(K, NormalDistribution(1200, 100)); -1
 

We replicate the experiment 5000 times. and record the final population values. 

> `assign`(runs, [seq(logirand(r, K, 1, T), j = 1 .. 5000)]); -1
 

> `assign`(tops, [seq(runs[j][T], j = 1 .. nops(runs))]); -1
 

> DataSummary(tops); 1
 

[mean = 1140.958515, standarddeviation = 87.66163470, skewness = -1.358680079, kurtosis = 7.292290075, minimum = 437.2256243, maximum = 1393.714210, cumulativeweight = 5000.]
[mean = 1140.958515, standarddeviation = 87.66163470, skewness = -1.358680079, kurtosis = 7.292290075, minimum = 437.2256243, maximum = 1393.714210, cumulativeweight = 5000.]		 (9)
 

Guess 1: The final population sizes follow a normal distribution (as the final population sizes are close to the carrying capacity when r is not too big) 

> `assign`(pl1, plot(PDF(Normal(Mean(tops), StandardDeviation(tops)), t), t = min(tops) .. max(tops), color = red, thickness = 2)); -1
`assign`(pl1, plot(PDF(Normal(Mean(tops), StandardDeviation(tops)), t), t = min(tops) .. max(tops), color = red, thickness = 2)); -1
 

> `assign`(H, Histogram(tops, frequencyscale = relative, averageshifted = 4)); -1
 

> display([H, pl1]); 1
 

Plot_2d
 

Guess 2:The log(final popsizes) follow normal distribution 

> `assign`(logtops, [seq(ln(tops[j]), j = 1 .. nops(runs))]); -1
 

> DataSummary(logtops); 1
 

[mean = 7.036379195, standarddeviation = 0.8277494682e-1, skewness = -2.095428191, kurtosis = 13.00661230, minimum = 6.080449365, maximum = 7.239727556, cumulativeweight = 5000.]
[mean = 7.036379195, standarddeviation = 0.8277494682e-1, skewness = -2.095428191, kurtosis = 13.00661230, minimum = 6.080449365, maximum = 7.239727556, cumulativeweight = 5000.]		 (10)
 

> `assign`(pl1, plot(PDF(Normal(Mean(logtops), StandardDeviation(logtops)), t), t = min(logtops) .. max(logtops), color = red, thickness = 2)); -1
`assign`(pl1, plot(PDF(Normal(Mean(logtops), StandardDeviation(logtops)), t), t = min(logtops) .. max(logtops), color = red, thickness = 2)); -1
 

> `assign`(Hlog, Histogram(logtops, frequencyscale = relative, averageshifted = 4)); -1
 

> display([Hlog, pl1]); 1
 

Plot_2d
 

None of the guesses look good 

Here is a graph of the first runs of the model 

> `assign`(S, [seq([seq(runs[k][t], k = 1 .. `+`(`*`(`/`(1, 50), `*`(nops(runs)))))], t = 1 .. T)]); -1
 

> tempplot(S); 1
 

Plot_2d
 

>