# Population growth model using cohort data method: Is there a way to make the repetitive code simpler?

This is for a little simple project in which we try to model population of a certain species of fish and the impacts of overfishing

I do find it rather rudimentary to write it all out as such and was thinking if there is a way to code these equations in a shorter version

PopulationlistBFT[T_, γ_, Smax_, List0_, S111_, S220_,  S30_,
S40_, S50_, S60_, S70_,  S80_, S90_, S100_,  S110_, S120_, S130_,
S140_, S150_, S160_, S170_, S180_, S190_,
S200_, ρ1_, ρ2_, ρ3_ , ρ4_ , ρ5_, \
ρ15_, H1_, H2_, H3_, H4_, H5_, H6_, H7_, H8_, H9_,  H10_, H11_,
H12_, H13_, H14_, H15_, H16_, H17_, H18_, H19_, H20_, W1_, W2_, W3_,
W4_, W5_, W6_, W7_, W8_, W9_, W10_, W11_, W12_, W13_, W14_, W15_,
W16_, W17_, W18_, W19_, W20_] :=

For[t = 0;  {S1 = S111, S2 = S220, S3 = S30, S4 = S40, S5 = S50,
S6 = S60, S7 = S70, S8 = S80, S9 = S90, S10 = S100, S11 = S110,
S12 = S120, S13 = S130, S14 = S140, S15 = S150, S16 = S160,
S17 = S170, S18 = S180, S19 = S190, S20 = S200}, t < T + 1,
t++, {Stotal =

S1 + S2 + S3 + S4 + S5 + S6 + S7 + S8 + S9 + S10 + S11 + S12 +
S13 + S14 + S15 + S16 + S17 + S18 + S19 - S20,
Sjuvenile = S1 + S2 + S3 + S4 + S5,
Sm = S6 + S7 + S8 + S9 + S10 + S11 + S12 + S13 + S14 + S15 + S16 +
S17 + S18 + S19 - S20,
YTotal =
Y1 + Y2 + Y3 + Y4 + Y5 + Y6 + Y7 + Y8 + Y9 + Y10 + Y11 + Y12 +
Y13 + Y14 + Y15 + Y16 + Y17 + Y18 + Y19 + Y20,
HTotal =
H1 + H2 + H3 + H4 + H5 + H6 + H7 + H8 + H9 + H10 + H11 + H12 +
H13 + H14 + H15 + H16 + H17 + H18 + H19 + H20,

S20 = ρ15 ((S19) - If[H20 < S19, H20, H20/HTotal* If[S19 < 0, 0, S19]]),
Y20 = W20*If[H20 < S19, H20, H20/HTotal*If[S19 < 0, 0, S19]],
S19 = ρ15 ((S18) - If[H19 < S18, H19, H19/HTotal*If[S18 < 0, 0, S18]]),
Y19 = W19*If[H19 < S18, H19, H19/HTotal*If[S18 < 0, 0, S18]],
S18 = ρ15 ((S17) - If[H18 < S17, H18, H18/HTotal*If[S17 < 0, 0, S17]]),
Y18 = W18*If[H18 < S17, H18, H18/HTotal*If[S17 < 0, 0, S17]],
S17 = ρ15 ((S16) - If[H17 < S16, H17, H17/HTotal*If[S16 < 0, 0, S16]]),
Y17 = W17*If[H17 < S16, H17, H17/HTotal*If[S16 < 0, 0, S16]],
S16 = ρ15 ((S15) - If[H16 < S15, H16, H16/HTotal*If[S15 < 0, 0, S15]]),
Y16 = W16*If[H16 < S15, H16, H16/HTotal*If[S15 < 0, 0, S15]],
S15 = ρ15 ((S14) - If[H15 < S14, H15, H15/HTotal*If[S14 < 0, 0, S14]]),
Y15 = W15*If[H15 < S14, H15, H15/HTotal*If[S14 < 0, 0, S14]],
S14 = ρ5 ((S13) - If[H14 < S13, H14, H14/HTotal*If[S13 < 0, 0, S13]]),
Y14 = W14*If[H14 < S13, H14, H14/HTotal*If[S13 < 0, 0, S13]],


And this goes on until S1

S3 = ρ2 ((S2) - If[H3 < S2, H3, H3/HTotal*If[S2 < 0, 0, S2]]),
Y3 = W3*If[H3 < S2, H3, H3/HTotal*If[S2 < 0, 0, S2]],
S2 = ρ1 ((S1) - If[H2 < S1, H2, H2/HTotal*If[S1 < 0, 0, S1]]),
Y2 = W2*If[H2 < S1, H2, H2/HTotal*If[S1 < 0, 0, S1]],
S1 = γ ( 1 -  Stotal/Smax  ) (Sm) - If[H1 < γ ( 1 -  Stotal/Smax  ) (Sm), H1,
H1/HTotal*If[γ ( 1 -  Stotal/Smax  ) (Sm) < 0, 0, γ ( 1 -  Stotal/Smax  ) (Sm)]],
Y1 = W1*If[H1 < γ ( 1 -  Stotal/Smax  ) (Sm), H1, H1/HTotal* If[γ ( 1 -  Stotal/Smax  ) (Sm) < 0, 0, γ ( 1 -  Stotal/Smax  ) (Sm)]]};

List0[t] = {t, If[Sm < 0, 0, Sm], If[Sjuvenile < 0, 0, Sjuvenile],
If[Stotal < 0, 0, Stotal], YTotal}]


I think part of knowing the answer is knowing how to phrase my questions better. Would appreciate any help in redirecting me to a more simplified version of the code.

Context

Thanks again for the comment! Really appreciate this.

To answer some of the comments, this is for a project for an Environmental Economics class. We were trying to model discrete-time (one year period) of the population growth of the Pacific Bluefin Tuna.

We started out at Year 1, in which have a γ variable. It acts as a growth factor for the increase in fish stocks. Secondly, we included the variable ( 1 - Stotal/Smax ) in which it's essentially penalising the growth factor if the stock approaches the maximum carrying capacity. Intuition of it is to model a Malthusian dynamic inside our data. Finally we find that the growth of fish would only depend on the matured stock as only they are able to reproduce. We gathered from the data that the matured stock would be >5 years of age.

What we wanted to capture with ρ is the non-fishing mortality rate. This is to capture the data that affects the fish such as climate change and other shocks in the ocean.

Finally for the harvesting. We thought that given the limitation of technology, it would be better to use an absolute amount in harvesting rate denoted by H. The H is different in each period and we gathered an average harvesting per year from the data that we used. But we also note that as population become sparse, fisherman would change methods and use a proportional harvesting method. The way we identified the proportions is just by taking the average harvest for each year and dividing it by the total harvest of all the fish.

Since this is an Economics course, we would be interested in changing fishing policies and modelling the yield. Thus, the Y. The yield is dependent on the catch and hence why it uses the catch numbers. The weights are different for each year the fish is in, hence why we have the W ranging from 1 - 20.

Thank you again for helping out. It's truly nice to have communities that are interested in helping to understand Mathematica.

Some Data

So the data that we pulled are from the International Scientific Committee for Tuna and Tuna-like Species. The data is as below

 PopulationlistBFT[T_, Y_, Smax_, List0_,
S111_, S220_,  S30_, S40_, S50_,
S60_, S70_, S80_, S90_, S100_,
S110_,S120_, S130_, S140_, S150_,
S160_, S170_, S180_, S190_, S200_,
ρ1_, ρ2_, ρ3_ , ρ4_ , ρ5_, ρ15_,
H1_, H2_, H3_, H4_, H5_,
H6_, H7_, H8_, H9_,  H10_,
H11_, H12_, H13_, H14_, H15_,
H16_, H17_, H18_, H19_, H20_,
W1_, W2_, W3_, W4_, W5_,
W6_, W7_, W8_, W9_, W10_,
W11_, W12_, W13_, W14_, W15_,
W16_, W17_, W18_, W19_, W20_]

PopulationlistBFT[50, 0.98, 50000, List0,
4000, 1000, 500, 500, 100,
60, 60, 20, 7, 1.1,
1.1, 1, 1, 1, 1,
1, 1, 0.2, 1, 1,
0.95, 0.99, 0.99, 0.99, 0.99, 0.95,
2400, 970, 375, 80, 20,
10, 10, 3, 1.2, 0.18,
0.18, 0.16, 0.16, 0.15, 0.15,
0.15, 0.15, 0, 0.15, 0.15,
5, 20, 30, 50, 65,
90, 100, 105, 125, 165,
185, 195, 200, 205, 210,
213, 216, 220, 222, 224]


From this data, I would actually developed a graphical representation to show the decline in stocks in about 30 years. Data seems bleak but it is what it is.

PopulationRunMatured = Table[{List0[j][], List0[j][]}, {j, 1, 50}];
PopulationRunJuvenile = Table[{List0[j][], List0[j][]}, {j, 1, 50}];
PopulationRunTotal = Table[{List0[j][], List0[j][]}, {j, 1, 50}];
PopulationYield = Table[{List0[j][], List0[j][]}, {j, 1, 50}];

ListPlot[PopulationRunMatured, Frame -> True, PlotLabel -> "Stock of Matured as a function of time", FrameLabel -> {"t","S[t] in '000"}, PlotRange -> Full, Joined -> True]
ListPlot[PopulationRunJuvenile, Frame -> True, PlotLabel -> "Stock of Juvenile as a function of time",  FrameLabel -> {"t","S[t] in '000"}, PlotRange -> Full, Joined -> True]
ListPlot[PopulationRunTotal, Frame -> True, PlotLabel -> "Stock of Total as a function of time", FrameLabel -> {"t","S[t] in '000"}, PlotRange -> Full, Joined -> True]
ListPlot[PopulationYield, Frame -> True, PlotLabel -> "Total Yield as a function of time", FrameLabel -> {"t","S[t] in '000"}, PlotRange -> Full, Joined -> True]


The expected output would be a graph similar to this Update on using Jagra's method

 PopulationlistBFT[T_, γ_, smax_, List0_, sInputs1_, sInputs2_, sInputs3_,
sInputs4_, sInputs5_, sInputs6_, sInputs7_, sInputs8_, sInputs9_,
sInputs10_, sInputs11_, sInputs12_, sInputs13_, sInputs14_, sInputs15_,
sInputs16_, sInputs17_, sInputs18_, sInputs19_, sInputs20_ , ρ1_, ρ2_, ρ3_,
ρ4_, ρ5_, ρ6_, ρ7_, ρ8_, ρ9_, ρ10_, ρ11_, ρ12_, ρ13_, ρ14_, ρ15_, ρ16_,
ρ17_,ρ18_, ρ19_, ρ20_, h1_, h2_, h3_, h4_, h5_, h7_, h8_, h9_, h10_, h11_,
h12_, h13_, h14_, h15_, h16_, h17_, h18_, h19_, h20_, w1_, w2_, w3_, w4_,
w5_, w6_, w7_, w8_, w9_, w10_, w11_, w12_, w13_, w14_, w15_, w16_, w17_,
w18_, w19_, w20_] :=
For[t = 0 ;
{s1 = sInputs1, s2 = sInputs2, s3 = sInputs4, s4 = sInputs4,
s5 = sInputs5, s6 = sInputs6, s7 = sInputs8, s9 = sInputs9,
s10 = sInputs10, s11 = sInputs11, s12 = sInputs12, s13 = sInputs13,
s14 = sInputs14, s15 = sInputs15, s16 = sInputs16,
s17 = sInputs17, s18 = sInputs18, s19 = sInputs19,
s20 = sInputs20} , t < T + 1,
t++,
{sTotal = Total[s[[1 ;; 19]]] - s[],
sJuvenile = Total[s[[1 ;; 5]]],
sm = Total[s[[6 ;; 19]]] - s[],
yTotal = Total[y[[1;;19]], hTotal = Total[h[[1;;20]]],

s[[20 - i + 1]] = Table[ρ[[20 - i + 1]] ((s[[20 - i]]) - If[h[[20 - i + 1]] < s[[20 - i]], h[[20 - i + 1]], h[[20 - i + 1]]/hTotal* If[s[[20 - i]] < 0, 0, s[[20 - i]]]]), {i, 1, 19}],
y[[20 - i + 1]] = Table[w[[20 - i + 1]]* If[h[[20 - i + 1]] < s[[20 - i]], h[[20 - i + 1]], h[[20 - i + 1]]/hTotal*If[s[[20 - i]] < 0, 0, s[[20 - i]]]], {i, 1, 19}],
s[] = γ (1 - sTotal/smax) (sm)-If[h[] < γ (1 - sTotal/smax) (sm), h[], h[]/hTotal*If[γ (1 - sTotal/smax) (sm) < 0, 0, γ (1 - sTotal/smax) (sm)]],
y[] = w[]*If[h[] < γ (1 - sTotal/smax) (sm), h[],h[]/hTotal* If[γ(1 - sTotal/smax) (sm) < 0, 0, γ (1 - sTotal/smax) (sm)]]};
List0[t] = {t, If[sm < 0, 0, sm], If[sJuvenile < 0, 0, sJuvenile], If[sTotal < 0, 0, sTotal], yTotal}]


I'm pretty sure I got the syntax quite wrong. Because when I input the data accordingly, it wasn't able to detect the output.

PopulationlistBFT[50, 0.98, 50000, List0,
4000, 1000, 500, 500, 100,
60, 60, 20, 7, 1.1,
1.1, 1, 1, 1, 1,
1, 1, 0.2, 1, 1,
0.95, 0.99, 0.99, 0.99, 0.99,
0.99, 0.99, 0.99, 0.99, 0.99,
0.99, 0.99, 0.99, 0.99, 0.95,
0.95, 0.95, 0.95, 0.95, 0.95,
2400, 970, 375, 80, 20,
10, 10, 3, 1.2, 0.18,
0.18, 0.16, 0.16, 0.15, 0.15,
0.15, 0.15, 0, 0.15, 0.15,
5, 20, 30, 50, 65,
90, 100, 105, 125, 165,
185, 195, 200, 205, 210,
213, 216, 220, 222, 224]

PopulationRunMatured = Table[{List0[j][], List0[j][]}, {j, 1, 50}];
PopulationRunJuvenile = Table[{List0[j][], List0[j][]}, {j, 1, 50}];
PopulationRunTotal = Table[{List0[j][], List0[j][]}, {j, 1, 50}];
PopulationYield = Table[{List0[j][], List0[j][]}, {j, 1, 50}];

ListPlot[PopulationRunMatured, Frame -> True, PlotLabel -> "Stock of Matured as a function of time", FrameLabel -> {"t","S[t] in '000"}, PlotRange -> Full, Joined -> True]
ListPlot[PopulationRunJuvenile, Frame -> True, PlotLabel -> "Stock of Juvenile as a function of time",  FrameLabel -> {"t","S[t] in '000"}, PlotRange -> Full, Joined -> True]
ListPlot[PopulationRunTotal, Frame -> True, PlotLabel -> "Stock of Total as a function of time", FrameLabel -> {"t","S[t] in '000"}, PlotRange -> Full, Joined -> True]
ListPlot[PopulationYield, Frame -> True, PlotLabel -> "Total Yield as a function of time", FrameLabel -> {"t","S[t] in '000"}, PlotRange -> Full, Joined -> True]


I suspect that it's not able to detect the information for the next period. But I could be wrong. Many thanks to Jagra though.

• I hung on the title: if the code needs to be of some length, then if you can make it shorter, it doesn't need to be that long - that's a contradiction that I can't turn my head around :/ – corey979 Dec 21 '16 at 13:49
• Why not have Y as a list instead of all those Yn. Then Ytotal is just Total @ Y. – Kuba Dec 21 '16 at 13:51
• It will make this more fun if you can put your problem in context. What do you want to accomplish? What real world problem? – Jagra Dec 21 '16 at 14:27
• Also, you'll find it much easier and faster to get help on this sort of thing if you present the minimum amount of code that illustrates the issue. Reduce your Ss & Hs & Ws & Ys down to a few rather than so many. AND then supply some actual values to plug into the function. – Jagra Dec 21 '16 at 14:29
• Why not make the apparent lists of inputs to your function actual lists? Why does the For[] statement need to assign the inputs to other variables? It seems that the repetitive nature of the code would suit itself to the use of Table[]. You could additionally likely benefit from dispensing with the use of For[] in favor of something more functional. – Jagra Dec 21 '16 at 14:48

## 1 Answer

Just a stab at this, because I don't really think I appreciate the problem in its entirety, but consider something like the code below as a step towards simplification. I think you can replace much of our receptive code with the use of Table[]:

PopulationlistBFT[
T_, γ_, smax_, List0_,
sInputs_ (* 20 items *),
ρ_ (* 20 items originally designated 1-5, 15 *),
h_ (* 20 items *), w_ (* 20 items *)] :=

For[t = 0 (*start *);
{s = sInputs} (* test *),
t < T + 1, t++,

{
sTotal = Total[s[[1 ;; 19]]] - s[],
sJuvenile = Total[s[[1 ;; 5]]],
sm = Total[s[[6 ;; 19]]] - s[],

yTotal = Total[y],
hTotal = Total[h],

s = Table[ρ[] ((s[[20 - i]]) -If[h[[20 - i + 1]] < s[[20 - i]], h[[20 - i + 1]], h[[20 - i + 1]]/hTotal* If[s[[20 - i]] < 0, 0, s[[20 - i]]]]), {i, 1, 20}],
y = Table[w[[20 - i + 1]]* If[h[[20 - i + 1]] < s[[20 - i]], h[[20 - i + 1]], h[[20 - i + 1]]/hTotal * If[s[[20 - i]] < 0, 0, s[[20 - i]]]], {i, 1, 20}],

(* The above should get you to the following where you appear to have started to change how you use ρ ... You could use the same idea as above and just introduce something to iterate the values associtated with ρ  *)

s[] = ρ[] ((s[]) - If[h[] < s[], h[], h[]/hTotal*If[s[] < 0, 0, s[]]]),
y[] = w[]*If[h[] < s[], h[],h[]/hTotal*If[s[] < 0, 0, s[]]],

s[] = ρ[] ((s[]) - If[h[] < s[], h[], h[]/hTotal*If[s[] < 0, 0, s[]]]),
y[] = w[]*If[h[] < s[], h[], h[]/hTotal*If[s[] < 0, 0, s[]]],

s[] = γ (1 - sTotal/smax) (sm) - If[h[] < γ (1 - sTotal/smax) (sm), h[], h[]/hTotal* If[γ (1 - sTotal/smax) (sm) < 0, 0, γ (1 - sTotal/smax) (sm)]],
y[] = w[]*If[h[] < γ (1 - sTotal/smax) (sm), h[], h[]/hTotal* If[γ (1 - sTotal/smax) (sm) < 0,0, γ (1 - sTotal/smax) (sm)]]

};

List0[t] = {t,
If[sm < 0, 0, sm],
If[sJuvenile < 0, 0, sJuvenile],
If[sTotal < 0, 0, sTotal],
yTotal}
]


Note, I have changed your variable names so they do not begin with capital letters.

Also, you can likely do away with the For[] statement and replace it with something more functional. doing so will likely simplify my simplistic suggestions more and far more elegantly.

Also, please supply some data inputs and expected outputs.

• I've added some more information in the question. Thanks again for helping out. – Adam Azraei Dec 27 '16 at 15:57