<<Statistics`DiscreteDistributions`

<<Statistics`NormalDistribution`

grow[nlist_, 
          T_] := 
Module[{newlist, bdist, births, 
            t}, 
newlist  =  
            nlist; 
Do[
bdist  =  
              BinomialDistribution[Last[newlist][[2]], 
                alpha]; 
births  = 
              Random[bdist]; 
AppendTo[
              newlist, {t, Last[newlist][[2]] + 
                  births}], 
{t, 1, 
              T}
]; 
newlist
];

alpha  =  E^(0.01) - 1;
Print[alpha];

T  =  50;

    n0  =  Round[6.6*10];

    nlist  =  {{0, n0}};


    nlist  =  grow[nlist, T];


    dataplot  = 
ListPlot[
nlist, 
PlotStyle -> 
            RGBColor[1, 0, 0], 
AxesOrigin  -> {0, 
              0}, 
AxesStyle  ->  {Thickness[
                0.0025]}, 
TextStyle -> {FontSize -> 18}, 
Prolog -> {PointSize[0.01], 
              Thickness[0.0025]}, 
DisplayFunction -> 
            Identity
];


    point  =  
        Graphics[{PointSize[0.0000001], Point[{0, 0}]}];

    Show[
point, 
        dataplot, 
PlotRange -> 
          All, 
Axes -> 
          True, 
AxesOrigin  -> {0, 
            0}, 
AxesStyle  ->  {Thickness[
              0.0025]}, 
TextStyle -> {FontSize -> 18}, 
Prolog -> {PointSize[0.01], 
            Thickness[
              0.0025]}, 
DisplayFunction -> $DisplayFunction
];

expectedplot  = 
Plot[
n0*(1 + 
                  alpha)^t, {t, 0, T}, 
PlotStyle -> 
            RGBColor[0, 0, 1], 
AxesOrigin  -> {0, 
              0}, 
AxesStyle  ->  {Thickness[
                0.0025]}, 
TextStyle -> {FontSize -> 18}, 
Prolog -> {PointSize[0.01], 
              Thickness[0.0025]}, 
DisplayFunction -> 
            Identity
];

point  =  
        Graphics[{PointSize[0.0000001], 
            Point[{0, 0}]}];
Show[
point, expectedplot, 
        dataplot, 
PlotRange -> 
          All, 
Axes -> 
          True, 
AxesOrigin  -> {0, 
            0}, 
AxesStyle  ->  {Thickness[
              0.0025]}, 
TextStyle -> {FontSize -> 18}, 
Prolog -> {PointSize[0.01], 
            Thickness[
              0.0025]}, 
DisplayFunction -> $DisplayFunction
];

\(\(\(
\)\(\(nlist\  = \ {{0, 
            n0}};\)


    \(nlist\  = \ grow[nlist, T];\)


    \(dataplot\  = 
ListPlot[
nlist, 
PlotStyle -> 
            RGBColor[1, 0, 0], 
AxesOrigin\  -> {0, 
              0}, 
AxesStyle\  -> \ {Thickness[
                0.0025]}, 
TextStyle -> {FontSize -> 18}, 
Prolog -> {PointSize[0.01], 
              Thickness[0.0025]}, 
DisplayFunction -> 
            Identity
];\)


    \(expectedplot\  = 
Plot[
n0*\((1 + 
                  alpha)\)^t, {t, 0, T}, 
PlotStyle -> 
            RGBColor[0, 0, 1], 
AxesOrigin\  -> {0, 
              0}, 
AxesStyle\  -> \ {Thickness[
                0.0025]}, 
TextStyle -> {FontSize -> 18}, 
Prolog -> {PointSize[0.01], 
              Thickness[0.0025]}, 
DisplayFunction -> 
            Identity
];\)


    point\  = \ Graphics[{PointSize[0.0000001], Point[{0, 0}]}]; 

    \(Show[
point, expectedplot, 
        dataplot, 
PlotRange -> 
          All, 
Axes -> 
          True, 
AxesOrigin\  -> {0, 
            0}, 
AxesStyle\  -> \ {Thickness[
              0.0025]}, 
TextStyle -> {FontSize -> 18}, 
Prolog -> {PointSize[0.01], 
            Thickness[
              0.0025]}, 
DisplayFunction -> $DisplayFunction
];\)\)\)\)

\(\(\(
\)\(
\)\(\(counts[
          finalpop_]\  := 
Module[{distinct, 
            i}, 
distinct\  = \ 
            Sort[Union[
                finalpop]]; 
Table[
{distinct[\([\)\(i\)\(]\)], 
              Count[finalpop, distinct[\([\)\(i\)\(]\)]]}, {i, 1, 
              Length[distinct]}
]
];\)\(
\)
    \)\)\)

Nruns  =  10000;
T  =  50;
\(\(n0\  = \ Round[6.6*10];\)
\)
finalpop  =  {};
Do[
nlist  =  {{0, 
              n0}}; 
nlist  =  
          grow[nlist, T]; 
AppendTo[
          finalpop, Last[
              nlist][[2]]], 
{run, 1, 
          Nruns}
];

mean  =  Mean[finalpop];
var  =  Variance[finalpop];
Print["mean = ", N[mean], "; predicted = ", 
        n0*(1 + alpha)^T];
Print["variance = ", N[var]];
Print[];

\(\(popcounts\  = \ counts[finalpop];\)
\)
dataplot  = 
ListPlot[
popcounts, 
PlotStyle -> 
            RGBColor[1, 0, 
              0], 
AxesOrigin  -> {n0*(1 + alpha)^
                  T, 0}, 
AxesStyle  ->  {Thickness[
                0.0025]}, 
TextStyle -> {FontSize -> 18}, 
Prolog -> {PointSize[0.01], 
              Thickness[0.0025]}, 
DisplayFunction -> 
            Identity
];

fitplot  = 
Plot[
Nruns*
            PDF[NormalDistribution[mean, Sqrt[var]], x], {x, 
            mean - 4  Sqrt[var], 
            mean + 4*Sqrt[var]}, 
PlotStyle -> 
            RGBColor[0, 0, 
              1], 
AxesOrigin  -> {n0*(1 + alpha)^
                  T, 0}, 
AxesStyle  ->  {Thickness[
                0.0025]}, 
TextStyle -> {FontSize -> 18}, 
Prolog -> {PointSize[0.01], 
              Thickness[0.0025]}, 
DisplayFunction -> 
            Identity
];

point  =  
        Graphics[{PointSize[0.0000001], 
            Point[{n0*(1 + alpha)^T, 0}]}];
Show[
point, fitplot, 
        dataplot, 
PlotRange -> 
          All, 
Axes -> 
          True, 
AxesOrigin  -> {n0*(1 + alpha)^
                T, 
            0}, 
AxesStyle  ->  {Thickness[
              0.0025]}, 
TextStyle -> {FontSize -> 18}, 
Prolog -> {PointSize[0.01], 
            Thickness[
              0.0025]}, 
DisplayFunction -> $DisplayFunction
];

\(\(\(
\)\(\(Nruns\  = \ 10000;\)

    \(T\  = \ 50;\)

    \(n0\  = \ Round[6.6*10^2];\)


    \(finalpop\  = \ {};\)

    \(Do[
nlist\  = \ {{0, 
              n0}}; 
nlist\  = \ 
          grow[nlist, T]; 
AppendTo[
          finalpop, \(Last[
              nlist]\)[\([\)\(2\)\(]\)]], 
{run, 1, 
          Nruns}
];\)\)\)\)

mean  =  Mean[finalpop];
var  =  Variance[finalpop];
Print["mean = ", N[mean], "; predicted = ", 
        n0*(1 + alpha)^T];
Print["variance = ", N[var]];
Print[];

\(\(popcounts\  = \ counts[finalpop];\)
\)
dataplot  = 
ListPlot[
popcounts, 
PlotStyle -> 
            RGBColor[1, 0, 
              0], 
AxesOrigin  -> {n0*(1 + alpha)^
                  T, 0}, 
AxesStyle  ->  {Thickness[
                0.0025]}, 
TextStyle -> {FontSize -> 18}, 
Prolog -> {PointSize[0.01], 
              Thickness[0.0025]}, 
DisplayFunction -> 
            Identity
];

fitplot  = 
Plot[
Nruns*
            PDF[NormalDistribution[mean, Sqrt[var]], x], {x, 
            mean - 4  Sqrt[var], 
            mean + 4*Sqrt[var]}, 
PlotStyle -> 
            RGBColor[0, 0, 
              1], 
AxesOrigin  -> {n0*(1 + alpha)^
                  T, 0}, 
AxesStyle  ->  {Thickness[
                0.0025]}, 
TextStyle -> {FontSize -> 18}, 
Prolog -> {PointSize[0.01], 
              Thickness[0.0025]}, 
DisplayFunction -> 
            Identity
];

point  =  
        Graphics[{PointSize[0.0000001], 
            Point[{n0*(1 + alpha)^T, 0}]}];
Show[
point, fitplot, 
        dataplot, 
PlotRange -> 
          All, 
Axes -> 
          True, 
AxesOrigin  -> {n0*(1 + alpha)^
                T, 
            0}, 
AxesStyle  ->  {Thickness[
              0.0025]}, 
TextStyle -> {FontSize -> 18}, 
Prolog -> {PointSize[0.01], 
            Thickness[
              0.0025]}, 
DisplayFunction -> $DisplayFunction
];

\(\(\(
\)\(\(Nruns\  = \ 10000;\)

    \(T\  = \ 50;\)

    \(n0\  = \ Round[6.6*10^3];\)


    \(finalpop\  = \ {};\)

    \(Do[
nlist\  = \ {{0, 
              n0}}; 
nlist\  = \ 
          grow[nlist, T]; 
AppendTo[
          finalpop, \(Last[
              nlist]\)[\([\)\(2\)\(]\)]], 
{run, 1, 
          Nruns}
];\)\)\)\)

mean  =  Mean[finalpop];
var  =  Variance[finalpop];
Print["mean = ", N[mean], "; predicted = ", 
        n0*(1 + alpha)^T];
Print["variance = ", N[var]];
Print[];

\(\(popcounts\  = \ counts[finalpop];\)
\)
dataplot  = 
ListPlot[
popcounts, 
PlotStyle -> 
            RGBColor[1, 0, 
              0], 
AxesOrigin  -> {n0*(1 + alpha)^
                  T, 0}, 
AxesStyle  ->  {Thickness[
                0.0025]}, 
TextStyle -> {FontSize -> 18}, 
Prolog -> {PointSize[0.01], 
              Thickness[0.0025]}, 
DisplayFunction -> 
            Identity
];

fitplot  = 
Plot[
Nruns*
            PDF[NormalDistribution[mean, Sqrt[var]], x], {x, 
            mean - 4  Sqrt[var], 
            mean + 4*Sqrt[var]}, 
PlotStyle -> 
            RGBColor[0, 0, 
              1], 
AxesOrigin  -> {n0*(1 + alpha)^
                  T, 0}, 
AxesStyle  ->  {Thickness[
                0.0025]}, 
TextStyle -> {FontSize -> 18}, 
Prolog -> {PointSize[0.01], 
              Thickness[0.0025]}, 
DisplayFunction -> 
            Identity
];

point  =  
        Graphics[{PointSize[0.0000001], 
            Point[{n0*(1 + alpha)^T, 0}]}];
Show[
point, fitplot, 
        dataplot, 
PlotRange -> 
          All, 
Axes -> 
          True, 
AxesOrigin  -> {n0*(1 + alpha)^
                T, 
            0}, 
AxesStyle  ->  {Thickness[
              0.0025]}, 
TextStyle -> {FontSize -> 18}, 
Prolog -> {PointSize[0.01], 
            Thickness[
              0.0025]}, 
DisplayFunction -> $DisplayFunction
];