######## ######## ########
######## SIBMPOP4 ########
######## ######## ########
# SIBMPOP4 Version 2.1. beta
# (Simple: with all plots, not as function, no output table)
# (For academic purposes only. SIBMPOP4 comes with absolutely no guarantee)
# (R. Schwamborn, 2017)
#####################
# SIBMPOP4 #
####
##
#
# Age-structured, Stochastic Individual-Based Population Model (SIBM-POP4),
# with well-defined cohorts and non-seasonal von Bertalanffy growth,
# with fully stochastic mortality and growth.
# With error terms (% C.V. of Linf, K, Z and cohort strength)
# up to 11 cohorts
# seasonal recruitment
# adjustable duration of the recruitment season (days)
# Analysis of the Catch Curve
# with a "modified" Powell-Wetherall plot
# and output of bias (%) for Z/K and Linf estimates
############
## INPUTS ##
############
#Growth
Linf_input = 10
cvLinf_input = 5 # C.V,in % (e.g. 10 = 10%)
K_input = 1.0
cvK_input = 5 # C.V,in % (e.g. 3 = 3%)
t0_input = -0.1
# Mortality
Z_input = 2 # TOTAL Mortality Z
cvZ_input = 0 # C.V,in % (e.g. 10 = 10%)
# N Recruits per cohort
Nr_input = 1000 # N of recruits per year
cvNr_input = 0 # C.V. of recruitment in % (e.g. 30 = 30%)
# Duration (i.e., width) of the Recruitment season
cvt0_input = 180 # C.V. of t0, in % (e.g. 0 = 0%) -> used for Recruitment variab.
sdt0 <- abs(t0_input * (cvt0_input/100)) # standard devation of t0
(Days.season = round((2* (sdt0 <- abs(t0_input * (cvt0_input/100)))*365), 0)) # Recruitment season (days)
## [1] 131
# Days.season: Recruitment season (days) = 2 * standard deviation of the recruitment curve, in days.
###
#calculate SDs and random vectors
sdLinf <- Linf_input* (cvLinf_input/100)
Linfe_vec <- rnorm(10000, mean = Linf_input, sd = sdLinf)
sdK <- K_input * (cvK_input/100)
Ke_vec <- rnorm(10000, mean = K_input, sd = sdK)
sdt0 <- abs(t0_input * (cvt0_input/100))
t0e_vec <- rnorm(10000, mean = t0_input, sd = sdt0)
sdZ <- Z_input * (cvZ_input/100)
Ze_vec <- rnorm(10000, mean = Z_input, sd = sdZ)
sdNr <- Z_input * (cvNr_input/100)
Nre_vec <- rnorm(10000, mean = Nr_input, sd = sdNr)
####
## Mortality # N(t) as afunction of Nzero, Z, age "t"
M.fun <- function(N_0,Z_in,t) {
N <- N_0 * exp(-t*(Z_in))
;
N
}
M.fun(1000, 1.8,1)
## [1] 165.2989
###
# Apply Mortality to a population
# 1. make a t-vector
t_vec_4 <- seq(0,4, by = 1/12) # 0, 1, etc = jan
t_vec_10 <-seq(0,10, by =1/12)
#2. calculate mortality vector (N(t))
Nt_vec1 <- rep(0, 121)
N0_in <- 100
Z_sd <- (cvZ_input/100) * Z_input
Z_in_vec <- rnorm(121, mean= Z_input, sd = Z_sd)
for (w in 0:121) {
Nt_vec1[w] <- M.fun(N0_in, Z_in_vec[w], t_vec_10[w])
}
plot(Nt_vec1 ~ t_vec_10)
length(t_vec_10)
## [1] 121
## (NOT now: without increasing N -> step by step mortality!)
#3. generate a group of 1000 L(t) values (1 cohort) that grow with t
L0 <- rep(0,1000)
## L(t): VBGF
# L(t) length.at.age
VBG.fun.full <- function(t, Linf, K, t0) {
L = (Linf * (1 - exp(-K*(t-t0)))); L
}
VBG.fun.full(0.7, 10, 1.9, -0.1)
## [1] 7.812881
# apply VBGF to vector
Lt_vec1 <- rep(0, 121)
for (w in 0:121) {
Lt_vec1[w] <- VBG.fun.full(t_vec_10[w], 10, 0.9, 0)
}
plot(Lt_vec1 ~ t_vec_10)
length(t_vec_10)
## [1] 121
# L matrix: matrix of 120 columns (months) and 1000 lines.
# define new matrix
Lmat3 <- matrix( 0, nrow = 1000, ncol = 121)
# fill with growth trajectories
Lmat3 <- matrix( 0, nrow = 1000, ncol = 121)
for (i in 1:1000) {
for (w in 0:121) {
Lmat3[i,w] <- VBG.fun.full(t_vec_10[w], sample(Linfe_vec,1),
sample(Ke_vec,1), sample(t0e_vec,1))
}
}
# Plots the cohort
#plot( Lmat3[1,1:121] ~ t_vec_10,
# col=rgb(0, 0.1, 1, 0.12), pch=1, cex = 0.4,
# ylim = c(0, 15))
#INPUT: mortality vector Nt_vec1
plot(Nt_vec1 ~ t_vec_10)
Nt_vec1
## [1] 1.000000e+02 8.464817e+01 7.165313e+01 6.065307e+01 5.134171e+01
## [6] 4.345982e+01 3.678794e+01 3.114032e+01 2.635971e+01 2.231302e+01
## [11] 1.888756e+01 1.598797e+01 1.353353e+01 1.145588e+01 9.697197e+00
## [16] 8.208500e+00 6.948345e+00 5.881647e+00 4.978707e+00 4.214384e+00
## [21] 3.567399e+00 3.019738e+00 2.556153e+00 2.163737e+00 1.831564e+00
## [26] 1.550385e+00 1.312373e+00 1.110900e+00 9.403563e-01 7.959944e-01
## [31] 6.737947e-01 5.703549e-01 4.827950e-01 4.086771e-01 3.459377e-01
## [36] 2.928300e-01 2.478752e-01 2.098218e-01 1.776104e-01 1.503439e-01
## [41] 1.272634e-01 1.077261e-01 9.118820e-02 7.718914e-02 6.533920e-02
## [46] 5.530844e-02 4.681758e-02 3.963023e-02 3.354626e-02 2.839630e-02
## [51] 2.403695e-02 2.034684e-02 1.722323e-02 1.457915e-02 1.234098e-02
## [56] 1.044641e-02 8.842699e-03 7.485183e-03 6.336071e-03 5.363368e-03
## [61] 4.539993e-03 3.843021e-03 3.253047e-03 2.753645e-03 2.330910e-03
## [66] 1.973073e-03 1.670170e-03 1.413768e-03 1.196729e-03 1.013009e-03
## [71] 8.574939e-04 7.258529e-04 6.144212e-04 5.200963e-04 4.402521e-04
## [76] 3.726653e-04 3.154544e-04 2.670264e-04 2.260329e-04 1.913328e-04
## [81] 1.619597e-04 1.370959e-04 1.160492e-04 9.823351e-05 8.315287e-05
## [86] 7.038739e-05 5.958164e-05 5.043477e-05 4.269211e-05 3.613809e-05
## [91] 3.059023e-05 2.589407e-05 2.191886e-05 1.855391e-05 1.570555e-05
## [96] 1.329446e-05 1.125352e-05 9.525897e-06 8.063498e-06 6.825603e-06
## [101] 5.777749e-06 4.890759e-06 4.139938e-06 3.504382e-06 2.966395e-06
## [106] 2.510999e-06 2.125515e-06 1.799210e-06 1.522998e-06 1.289190e-06
## [111] 1.091276e-06 9.237450e-07 7.819332e-07 6.618922e-07 5.602796e-07
## [116] 4.742665e-07 4.014579e-07 3.398268e-07 2.876572e-07 2.434965e-07
## [121] 2.061154e-07
Mort_as_Factor_vec <- (Nt_vec1/Nt_vec1[1])
Abs_vector <- length(Lmat3[,1]) - (Nt_vec1/Nt_vec1[1] * length(Lmat3[,1]))
Abs_vector <- round(Abs_vector,0) # vector with absentees in each month
length(Abs_vector)
## [1] 121
Abs_vector[3] # number of absentees in month 3
## [1] 283
N_Mort <- 50 # number of absentees in a month
for (j in 1:121) {
N_Mort <- Abs_vector[j]
N_Mort_lines_vec <-sample((seq(1:length(Lmat3[,1]))), N_Mort, replace = FALSE)
# replace absentees with NaN for column j
for (y in 1: N_Mort) {
Lmat3[(N_Mort_lines_vec[y]), j] <- NA
}
}
# View(Lmat3)
# Plots the cohort
# plot( Lmat3[1,1:121] ~ t_vec_10, col=rgb(0, 0.1, 1, 0.12),
# pch=1, cex = 0.4, ylim = c(0,15))
#
#for (i in 1:1000) {
# points( Lmat3[i,1:121] ~ t_vec_10, col=rgb(0, 0.1, 1, 0.17), pch=1, cex = 0.5)
#
#}
# more plots
plot(Lmat3[12,] ~ t_vec_10 )
# plot N(t) of the matrix # plot N_alive_vec
sum(is.na(Lmat3[,14])) # number ofNAs in column 14
## [1] 885
N_alive <- length(Lmat3[,14]) - sum(is.na(Lmat3[,14])) # number ofNAs in column 14
N_alive_vec <- rep(0, ncol(Lmat3))
for (j in 1: ncol(Lmat3)) {
N_alive_vec[j] <- length(Lmat3[,j]) - sum(is.na(Lmat3[,j])) # number ofNAs in column 14
}
plot(N_alive_vec ~ t_vec_10 , main = "N ind. in cohort")
text(8, 800, paste("Z_input=", Z_input))
N_alive_vec[(9*12)+1] # N alive at the onset of year 10
## [1] 0
# ( n= 27 ind. when Z = 0.4)
# 5. make 12 histograms, make a catch curve for the last year
# vector for "mids"
mids_OK_05 <- seq(0, 15, by = 0.5)
breaks_OK_05 <- seq(-0.2, 15.2, by = 0.5)
hist_mat3 <- matrix(0, nrow = (12), ncol = (length(breaks_OK_05)-1))
for (j in 110:121) {
hist(Lmat3[,j] , xlim = c(0, 15), ylim = c(0, 500), breaks = breaks_OK_05 )
obj <- hist(Lmat3[,j] , xlim = c(0, 15), ylim = c(0, 500), breaks = breaks_OK_05 )
hist_mat3[(j-109),(1:(length(obj$counts)))] <- obj$counts
}
# View(hist_mat3)
### OK!
# 6. repeat for n= 10 cohorts
# Fill ALL with NA columns, restart at t = 1 y, with t0 = t0+0 - OK!
# Fill ALL with NA columns, restart at t = 1 y, with t0 = t0+1
# Fill ALL with NA columns, restart at t = 1 y, with t0 = t0+2
# etc. ...
# Fill ALL with NA columns, restart at t = 1 y, with t0 = t0+9
# write a loop
# 6.2. Cohort 2
### replace absentees with NaN for ALL columns
# define new matrix
Lmat32 <- matrix( NA , nrow = 1000, ncol = 121)
# 62a. fill with growth trajectories
Lmat32 <- matrix( NA , nrow = 1000, ncol = 121)
for (i in 1:1000) {
for (w in 0:121) {
Lmat32[i,w] <- VBG.fun.full(t_vec_10[w], sample(Linfe_vec,1),
sample(Ke_vec,1), (1+ sample(t0e_vec,1)))
}
}
# 62b. clean negative L values (first columns)
Lmat32 <- ifelse(Lmat32<0, NA,Lmat32)
# View(Lmat32)
plot(Lmat32[12,] ~ t_vec_10 )
## 62c. apply Mortality
#INPUT: mortality vector Nt_vec1
plot(Nt_vec1 ~ t_vec_10)
Nt_vec1
## [1] 1.000000e+02 8.464817e+01 7.165313e+01 6.065307e+01 5.134171e+01
## [6] 4.345982e+01 3.678794e+01 3.114032e+01 2.635971e+01 2.231302e+01
## [11] 1.888756e+01 1.598797e+01 1.353353e+01 1.145588e+01 9.697197e+00
## [16] 8.208500e+00 6.948345e+00 5.881647e+00 4.978707e+00 4.214384e+00
## [21] 3.567399e+00 3.019738e+00 2.556153e+00 2.163737e+00 1.831564e+00
## [26] 1.550385e+00 1.312373e+00 1.110900e+00 9.403563e-01 7.959944e-01
## [31] 6.737947e-01 5.703549e-01 4.827950e-01 4.086771e-01 3.459377e-01
## [36] 2.928300e-01 2.478752e-01 2.098218e-01 1.776104e-01 1.503439e-01
## [41] 1.272634e-01 1.077261e-01 9.118820e-02 7.718914e-02 6.533920e-02
## [46] 5.530844e-02 4.681758e-02 3.963023e-02 3.354626e-02 2.839630e-02
## [51] 2.403695e-02 2.034684e-02 1.722323e-02 1.457915e-02 1.234098e-02
## [56] 1.044641e-02 8.842699e-03 7.485183e-03 6.336071e-03 5.363368e-03
## [61] 4.539993e-03 3.843021e-03 3.253047e-03 2.753645e-03 2.330910e-03
## [66] 1.973073e-03 1.670170e-03 1.413768e-03 1.196729e-03 1.013009e-03
## [71] 8.574939e-04 7.258529e-04 6.144212e-04 5.200963e-04 4.402521e-04
## [76] 3.726653e-04 3.154544e-04 2.670264e-04 2.260329e-04 1.913328e-04
## [81] 1.619597e-04 1.370959e-04 1.160492e-04 9.823351e-05 8.315287e-05
## [86] 7.038739e-05 5.958164e-05 5.043477e-05 4.269211e-05 3.613809e-05
## [91] 3.059023e-05 2.589407e-05 2.191886e-05 1.855391e-05 1.570555e-05
## [96] 1.329446e-05 1.125352e-05 9.525897e-06 8.063498e-06 6.825603e-06
## [101] 5.777749e-06 4.890759e-06 4.139938e-06 3.504382e-06 2.966395e-06
## [106] 2.510999e-06 2.125515e-06 1.799210e-06 1.522998e-06 1.289190e-06
## [111] 1.091276e-06 9.237450e-07 7.819332e-07 6.618922e-07 5.602796e-07
## [116] 4.742665e-07 4.014579e-07 3.398268e-07 2.876572e-07 2.434965e-07
## [121] 2.061154e-07
Mort_as_Factor_vec <- (Nt_vec1/Nt_vec1[1])
Abs_vector <- length(Lmat32[,1]) - (Nt_vec1/Nt_vec1[1] * length(Lmat32[,1]))
Abs_vector <- round(Abs_vector,0) # vector with absentees in each month
length(Abs_vector)
## [1] 121
Abs_vector[3] # number of absentees in month 3
## [1] 283
N_Mort <- 50 # number of absentees in a month
#adjust start column or Mort.
# replace "(j in 1:121)" by "(j in (1+12):121)", # "Abs_vector[j]" by "Abs_vector[j-12]"
for (j in (1+(1*12)):121) {
N_Mort <- Abs_vector[j-(1*12)]
N_Mort_lines_vec <-sample((seq(1:length(Lmat32[,1]))), N_Mort, replace = FALSE)
# replace absentees with NaN for column j
for (y in 1: N_Mort) {
Lmat32[(N_Mort_lines_vec[y]), j] <- NA
}
}
# View(Lmat32)
#plots!
plot(Lmat32[12,] ~ t_vec_10 )
# plot N(t) of the matrix # plot N_alive_vec
sum(is.na(Lmat32[,14])) # number ofNAs in column 14
## [1] 277
N_alive <- length(Lmat32[,14]) - sum(is.na(Lmat32[,14])) # number ofNAs in column 14
N_alive_vec <- rep(0, ncol(Lmat32))
for (j in 1: ncol(Lmat32)) {
N_alive_vec[j] <- length(Lmat32[,j]) - sum(is.na(Lmat32[,j])) # number ofNAs in column 14
}
plot(N_alive_vec ~ t_vec_10 , main = "N ind. in cohort")
text(8, 800, paste("Z_input=", Z_input))
N_alive_vec[(9*12)+1] # N alive at the onset of year 10
## [1] 0
# coh 1( n= 27 ind. when Z = 0.4)
# coh 2( n= 41 ind. when Z = 0.4)
# 6.2d. make 12 histograms, make a catch curve for the last year
# vector for "mids"
mids_OK_05 <- seq(0, 15, by = 0.5)
breaks_OK_05 <- seq(-0.2, 15.2, by = 0.5)
hist_mat32 <- matrix(0, nrow = (12), ncol = (length(breaks_OK_05)-1))
for (j in 110:121) {
hist(Lmat32[,j] , xlim = c(0, 15), ylim = c(0, 500), breaks = breaks_OK_05 )
obj <- hist(Lmat32[,j] , xlim = c(0, 15), ylim = c(0, 500), breaks = breaks_OK_05 )
hist_mat32[(j-109),(1:(length(obj$counts)))] <- obj$counts
}
# View(hist_mat32)
### COHORT 2 # 6.2 OK!
###
### Now COHORT 3 # 6.3
### replace absentees with NaN for ALL columns
# define new matrix
Lmat33 <- matrix( NA , nrow = 1000, ncol = 121)
# 6.3a. fill with growth trajectories
Lmat33 <- matrix( NA , nrow = 1000, ncol = 121)
for (i in 1:1000) {
for (w in 0:121) {
Lmat33[i,w] <- VBG.fun.full(t_vec_10[w], sample(Linfe_vec,1),
sample(Ke_vec,1), (2+ sample(t0e_vec,1)))
}
}
# 6.3.b. clean negative L values (first columns)
Lmat33 <- ifelse(Lmat33<0, NA,Lmat33)
# View(Lmat33)
plot(Lmat33[12,] ~ t_vec_10 )
## 6.3.c. apply Mortality
#INPUT: mortality vector Nt_vec1
plot(Nt_vec1 ~ t_vec_10)
Nt_vec1
## [1] 1.000000e+02 8.464817e+01 7.165313e+01 6.065307e+01 5.134171e+01
## [6] 4.345982e+01 3.678794e+01 3.114032e+01 2.635971e+01 2.231302e+01
## [11] 1.888756e+01 1.598797e+01 1.353353e+01 1.145588e+01 9.697197e+00
## [16] 8.208500e+00 6.948345e+00 5.881647e+00 4.978707e+00 4.214384e+00
## [21] 3.567399e+00 3.019738e+00 2.556153e+00 2.163737e+00 1.831564e+00
## [26] 1.550385e+00 1.312373e+00 1.110900e+00 9.403563e-01 7.959944e-01
## [31] 6.737947e-01 5.703549e-01 4.827950e-01 4.086771e-01 3.459377e-01
## [36] 2.928300e-01 2.478752e-01 2.098218e-01 1.776104e-01 1.503439e-01
## [41] 1.272634e-01 1.077261e-01 9.118820e-02 7.718914e-02 6.533920e-02
## [46] 5.530844e-02 4.681758e-02 3.963023e-02 3.354626e-02 2.839630e-02
## [51] 2.403695e-02 2.034684e-02 1.722323e-02 1.457915e-02 1.234098e-02
## [56] 1.044641e-02 8.842699e-03 7.485183e-03 6.336071e-03 5.363368e-03
## [61] 4.539993e-03 3.843021e-03 3.253047e-03 2.753645e-03 2.330910e-03
## [66] 1.973073e-03 1.670170e-03 1.413768e-03 1.196729e-03 1.013009e-03
## [71] 8.574939e-04 7.258529e-04 6.144212e-04 5.200963e-04 4.402521e-04
## [76] 3.726653e-04 3.154544e-04 2.670264e-04 2.260329e-04 1.913328e-04
## [81] 1.619597e-04 1.370959e-04 1.160492e-04 9.823351e-05 8.315287e-05
## [86] 7.038739e-05 5.958164e-05 5.043477e-05 4.269211e-05 3.613809e-05
## [91] 3.059023e-05 2.589407e-05 2.191886e-05 1.855391e-05 1.570555e-05
## [96] 1.329446e-05 1.125352e-05 9.525897e-06 8.063498e-06 6.825603e-06
## [101] 5.777749e-06 4.890759e-06 4.139938e-06 3.504382e-06 2.966395e-06
## [106] 2.510999e-06 2.125515e-06 1.799210e-06 1.522998e-06 1.289190e-06
## [111] 1.091276e-06 9.237450e-07 7.819332e-07 6.618922e-07 5.602796e-07
## [116] 4.742665e-07 4.014579e-07 3.398268e-07 2.876572e-07 2.434965e-07
## [121] 2.061154e-07
Mort_as_Factor_vec <- (Nt_vec1/Nt_vec1[1])
Abs_vector <- length(Lmat33[,1]) - (Nt_vec1/Nt_vec1[1] * length(Lmat33[,1]))
Abs_vector <- round(Abs_vector,0) # vector with absentees in each month
length(Abs_vector)
## [1] 121
Abs_vector[3] # number of absentees in month 3
## [1] 283
N_Mort <- 50 # number of absentees in a month
#adjust start column or Mort.
# replace "(j in 1:121)" by "(j in (1+12):121)", # "Abs_vector[j]" by "Abs_vector[j-12]"
for (j in (1+(2*12)):121) {
N_Mort <- Abs_vector[j-(2*12)]
N_Mort_lines_vec <-sample((seq(1:length(Lmat33[,1]))), N_Mort, replace = FALSE)
# replace absentees with NaN for column j
for (y in 1: N_Mort) {
Lmat33[(N_Mort_lines_vec[y]), j] <- NA
}
}
# View(Lmat33)
#plots!
plot(Lmat33[12,] ~ t_vec_10 )
# plot N(t) of the matrix # plot N_alive_vec
sum(is.na(Lmat33[,14])) # number ofNAs in column 14
## [1] 1000
N_alive <- length(Lmat33[,14]) - sum(is.na(Lmat33[,14])) # number ofNAs in column 14
N_alive_vec <- rep(0, ncol(Lmat33))
for (j in 1: ncol(Lmat33)) {
N_alive_vec[j] <- length(Lmat33[,j]) - sum(is.na(Lmat33[,j])) # number ofNAs in column 14
}
plot(N_alive_vec ~ t_vec_10 , main = "N ind. in cohort")
text(8, 800, paste("Z_input=", Z_input))
N_alive_vec[(9*12)+1] # N alive at the onset of year 10
## [1] 0
# coh 1( n= 27 ind. when Z = 0.4)
# coh 2( n= 41 ind. when Z = 0.4)
# coh 3( n= 61 ind. when Z = 0.4)
# 6..3.d. make 12 histograms, make a catch curve for the last year
# vector for "mids"
mids_OK_05 <- seq(0, 15, by = 0.5)
breaks_OK_05 <- seq(-0.2, 15.2, by = 0.5)
hist_mat33 <- matrix(0, nrow = (12), ncol = (length(breaks_OK_05)-1))
for (j in 110:121) {
hist(Lmat33[,j] , xlim = c(0, 15), ylim = c(0, 500), breaks = breaks_OK_05 )
obj <- hist(Lmat33[,j] , xlim = c(0, 15), ylim = c(0, 500), breaks = breaks_OK_05 )
hist_mat33[(j-109),(1:(length(obj$counts)))] <- obj$counts
}
# View(hist_mat33)
### COHORT 3 # 6.3 OK!
###
### Now COHORT 4 # 6.4
### replace absentees with NaN for ALL columns
# define new matrix
Lmat34 <- matrix( NA , nrow = 1000, ncol = 121)
# 6.4.a. fill with growth trajectories
Lmat34 <- matrix( NA , nrow = 1000, ncol = 121)
for (i in 1:1000) {
for (w in 0:121) {
Lmat34[i,w] <- VBG.fun.full(t_vec_10[w], sample(Linfe_vec,1),
sample(Ke_vec,1), (3+ sample(t0e_vec,1)))
}
}
# 6.4.b. clean negative L values (first columns)
Lmat34 <- ifelse(Lmat34 < 0, NA,Lmat34)
# View(Lmat34)
plot(Lmat34[12,] ~ t_vec_10 )
## 6.4.c. apply Mortality
#INPUT: mortality vector Nt_vec1
plot(Nt_vec1 ~ t_vec_10)
Nt_vec1
## [1] 1.000000e+02 8.464817e+01 7.165313e+01 6.065307e+01 5.134171e+01
## [6] 4.345982e+01 3.678794e+01 3.114032e+01 2.635971e+01 2.231302e+01
## [11] 1.888756e+01 1.598797e+01 1.353353e+01 1.145588e+01 9.697197e+00
## [16] 8.208500e+00 6.948345e+00 5.881647e+00 4.978707e+00 4.214384e+00
## [21] 3.567399e+00 3.019738e+00 2.556153e+00 2.163737e+00 1.831564e+00
## [26] 1.550385e+00 1.312373e+00 1.110900e+00 9.403563e-01 7.959944e-01
## [31] 6.737947e-01 5.703549e-01 4.827950e-01 4.086771e-01 3.459377e-01
## [36] 2.928300e-01 2.478752e-01 2.098218e-01 1.776104e-01 1.503439e-01
## [41] 1.272634e-01 1.077261e-01 9.118820e-02 7.718914e-02 6.533920e-02
## [46] 5.530844e-02 4.681758e-02 3.963023e-02 3.354626e-02 2.839630e-02
## [51] 2.403695e-02 2.034684e-02 1.722323e-02 1.457915e-02 1.234098e-02
## [56] 1.044641e-02 8.842699e-03 7.485183e-03 6.336071e-03 5.363368e-03
## [61] 4.539993e-03 3.843021e-03 3.253047e-03 2.753645e-03 2.330910e-03
## [66] 1.973073e-03 1.670170e-03 1.413768e-03 1.196729e-03 1.013009e-03
## [71] 8.574939e-04 7.258529e-04 6.144212e-04 5.200963e-04 4.402521e-04
## [76] 3.726653e-04 3.154544e-04 2.670264e-04 2.260329e-04 1.913328e-04
## [81] 1.619597e-04 1.370959e-04 1.160492e-04 9.823351e-05 8.315287e-05
## [86] 7.038739e-05 5.958164e-05 5.043477e-05 4.269211e-05 3.613809e-05
## [91] 3.059023e-05 2.589407e-05 2.191886e-05 1.855391e-05 1.570555e-05
## [96] 1.329446e-05 1.125352e-05 9.525897e-06 8.063498e-06 6.825603e-06
## [101] 5.777749e-06 4.890759e-06 4.139938e-06 3.504382e-06 2.966395e-06
## [106] 2.510999e-06 2.125515e-06 1.799210e-06 1.522998e-06 1.289190e-06
## [111] 1.091276e-06 9.237450e-07 7.819332e-07 6.618922e-07 5.602796e-07
## [116] 4.742665e-07 4.014579e-07 3.398268e-07 2.876572e-07 2.434965e-07
## [121] 2.061154e-07
Mort_as_Factor_vec <- (Nt_vec1/Nt_vec1[1])
Abs_vector <- length(Lmat34[,1]) - (Nt_vec1/Nt_vec1[1] * length(Lmat34[,1]))
Abs_vector <- round(Abs_vector,0) # vector with absentees in each month
length(Abs_vector)
## [1] 121
Abs_vector[3] # number of absentees in month 3
## [1] 283
N_Mort <- 50 # number of absentees in a month
#adjust start column or Mort.
# replace "(j in 1:121)" by "(j in (1+12):121)", # "Abs_vector[j]" by "Abs_vector[j-12]"
for (j in (1+(3*12)):121) {
N_Mort <- Abs_vector[j-(3*12)]
N_Mort_lines_vec <-sample((seq(1:length(Lmat34[,1]))), N_Mort, replace = FALSE)
# replace absentees with NaN for column j
for (y in 1: N_Mort) {
Lmat34[(N_Mort_lines_vec[y]), j] <- NA
}
}
# View(Lmat34)
#plots!
plot(Lmat34[12,] ~ t_vec_10 )
# plot N(t) of the matrix # plot N_alive_vec
sum(is.na(Lmat34[,14])) # number ofNAs in column 14
## [1] 1000
N_alive <- length(Lmat34[,14]) - sum(is.na(Lmat34[,14])) # number ofNAs in column 14
N_alive_vec <- rep(0, ncol(Lmat34))
for (j in 1: ncol(Lmat34)) {
N_alive_vec[j] <- length(Lmat34[,j]) - sum(is.na(Lmat34[,j])) # number ofNAs in column 14
}
plot(N_alive_vec ~ t_vec_10 , main = "N ind. in cohort")
text(8, 800, paste("Z_input=", Z_input))
N_alive_vec[(9*12)+1] # N alive at the onset of year 10
## [1] 0
# coh 1( n= 27 ind. when Z = 0.4)
# coh 2( n= 41 ind. when Z = 0.4)
# coh 3( n= 61 ind. when Z = 0.4)
# coh 4( n= 91 ind. when Z = 0.4)
# percent Coh 1 at the onset of year 10
Ncoh1_ini10 <- 27
NTOTAL_ini10 <-sum (27,41,61,91)
NTOTAL_ini10
## [1] 220
Ncoh1_ini10 / NTOTAL_ini10
## [1] 0.1227273
Ncoh1_ini10 / NTOTAL_ini10 * 100 # COH1 in percent
## [1] 12.27273
# 6.4.d. make 12 histograms, make a catch curve for the last year
# vector for "mids"
mids_OK_05 <- seq(0, 15, by = 0.5)
breaks_OK_05 <- seq(-0.2, 15.2, by = 0.5)
hist_mat34 <- matrix(0, nrow = (12), ncol = (length(breaks_OK_05)-1))
for (j in 110:121) {
hist(Lmat34[,j] , xlim = c(0, 15), ylim = c(0, 500), breaks = breaks_OK_05 )
obj <- hist(Lmat34[,j] , xlim = c(0, 15), ylim = c(0, 500), breaks = breaks_OK_05 )
hist_mat34[(j-109),(1:(length(obj$counts)))] <- obj$counts
}
# View(hist_mat34)
### COHORT 4 # 6.4 OK!
########
### NOW COHORT 5 # 6.5
## check:
# 1.) increase t0 by 1
# 2.)increase by 1
#for (j in (1+(4*12)):121) {
#
# N_Mort <- Abs_vector[j-(4*12)]
### replace absentees with NaN for ALL columns
# define new matrix
Lmat35 <- matrix( NA , nrow = 1000, ncol = 121)
# 6.5.a. fill with growth trajectories
Lmat35 <- matrix( NA , nrow = 1000, ncol = 121)
for (i in 1:1000) {
for (w in 0:121) {
Lmat35[i,w] <- VBG.fun.full(t_vec_10[w], sample(Linfe_vec,1),
sample(Ke_vec,1), (4+ sample(t0e_vec,1)))
}
}
# 6.5.b. clean negative L values (first columns)
Lmat35 <- ifelse(Lmat35 < 0, NA,Lmat35)
# View(Lmat35)
plot(Lmat35[12,] ~ t_vec_10 )
## 6.5.c. apply Mortality
#INPUT: mortality vector Nt_vec1
plot(Nt_vec1 ~ t_vec_10)
Nt_vec1
## [1] 1.000000e+02 8.464817e+01 7.165313e+01 6.065307e+01 5.134171e+01
## [6] 4.345982e+01 3.678794e+01 3.114032e+01 2.635971e+01 2.231302e+01
## [11] 1.888756e+01 1.598797e+01 1.353353e+01 1.145588e+01 9.697197e+00
## [16] 8.208500e+00 6.948345e+00 5.881647e+00 4.978707e+00 4.214384e+00
## [21] 3.567399e+00 3.019738e+00 2.556153e+00 2.163737e+00 1.831564e+00
## [26] 1.550385e+00 1.312373e+00 1.110900e+00 9.403563e-01 7.959944e-01
## [31] 6.737947e-01 5.703549e-01 4.827950e-01 4.086771e-01 3.459377e-01
## [36] 2.928300e-01 2.478752e-01 2.098218e-01 1.776104e-01 1.503439e-01
## [41] 1.272634e-01 1.077261e-01 9.118820e-02 7.718914e-02 6.533920e-02
## [46] 5.530844e-02 4.681758e-02 3.963023e-02 3.354626e-02 2.839630e-02
## [51] 2.403695e-02 2.034684e-02 1.722323e-02 1.457915e-02 1.234098e-02
## [56] 1.044641e-02 8.842699e-03 7.485183e-03 6.336071e-03 5.363368e-03
## [61] 4.539993e-03 3.843021e-03 3.253047e-03 2.753645e-03 2.330910e-03
## [66] 1.973073e-03 1.670170e-03 1.413768e-03 1.196729e-03 1.013009e-03
## [71] 8.574939e-04 7.258529e-04 6.144212e-04 5.200963e-04 4.402521e-04
## [76] 3.726653e-04 3.154544e-04 2.670264e-04 2.260329e-04 1.913328e-04
## [81] 1.619597e-04 1.370959e-04 1.160492e-04 9.823351e-05 8.315287e-05
## [86] 7.038739e-05 5.958164e-05 5.043477e-05 4.269211e-05 3.613809e-05
## [91] 3.059023e-05 2.589407e-05 2.191886e-05 1.855391e-05 1.570555e-05
## [96] 1.329446e-05 1.125352e-05 9.525897e-06 8.063498e-06 6.825603e-06
## [101] 5.777749e-06 4.890759e-06 4.139938e-06 3.504382e-06 2.966395e-06
## [106] 2.510999e-06 2.125515e-06 1.799210e-06 1.522998e-06 1.289190e-06
## [111] 1.091276e-06 9.237450e-07 7.819332e-07 6.618922e-07 5.602796e-07
## [116] 4.742665e-07 4.014579e-07 3.398268e-07 2.876572e-07 2.434965e-07
## [121] 2.061154e-07
Mort_as_Factor_vec <- (Nt_vec1/Nt_vec1[1])
Abs_vector <- length(Lmat35[,1]) - (Nt_vec1/Nt_vec1[1] * length(Lmat35[,1]))
Abs_vector <- round(Abs_vector,0) # vector with absentees in each month
length(Abs_vector)
## [1] 121
Abs_vector[3] # number of absentees in month 3
## [1] 283
N_Mort <- 50 # number of absentees in a month
#adjust start column or Mort.
# replace "(j in 1:121)" by "(j in (1+12):121)", # "Abs_vector[j]" by "Abs_vector[j-12]"
for (j in (1+(4*12)):121) {
N_Mort <- Abs_vector[j-(4*12)]
N_Mort_lines_vec <-sample((seq(1:length(Lmat35[,1]))), N_Mort, replace = FALSE)
# replace absentees with NaN for column j
for (y in 1: N_Mort) {
Lmat35[(N_Mort_lines_vec[y]), j] <- NA
}
}
# View(Lmat35)
#plots!
plot(Lmat35[12,] ~ t_vec_10 )
# plot N(t) of the matrix # plot N_alive_vec
sum(is.na(Lmat35[,15])) # number ofNAs in column 15
## [1] 1000
N_alive <- length(Lmat35[,15]) - sum(is.na(Lmat35[,15])) # number ofNAs in column 15
N_alive_vec <- rep(0, ncol(Lmat35))
for (j in 1: ncol(Lmat35)) {
N_alive_vec[j] <- length(Lmat35[,j]) - sum(is.na(Lmat35[,j])) # number ofNAs in column j
}
plot(N_alive_vec ~ t_vec_10 , main = "N ind. in cohort")
text(8, 800, paste("Z_input=", Z_input))
N_alive_vec[(9*12)+1] # N alive at the onset of year 10
## [1] 0
# coh 1( n= 27 ind. when Z = 0.4)
# coh 2( n= 41 ind. when Z = 0.4)
# coh 3( n= 61 ind. when Z = 0.4)
# coh 4( n= 91 ind. when Z = 0.4)
# coh 5( n= 135 ind. when Z = 0.4)
# percent Coh 1 at the onset of year 10
Ncoh1_ini10 <- 27
NTOTAL_ini10 <-sum (27,41,61,91,135)
NTOTAL_ini10
## [1] 355
Ncoh1_ini10 / NTOTAL_ini10
## [1] 0.07605634
Ncoh1_ini10 / NTOTAL_ini10 * 100 # COH1 in percent
## [1] 7.605634
# 6.5.d. make 12 histograms, make a catch curve for the last year
# vector for "mids"
mids_OK_05 <- seq(0, 15, by = 0.5)
breaks_OK_05 <- seq(-0.2, 15.2, by = 0.5)
hist_mat35 <- matrix(0, nrow = (12), ncol = (length(breaks_OK_05)-1))
for (j in 110:121) {
hist(Lmat35[,j] , xlim = c(0, 15), ylim = c(0, 500), breaks = breaks_OK_05 )
obj <- hist(Lmat35[,j] , xlim = c(0, 15), ylim = c(0, 500), breaks = breaks_OK_05 )
hist_mat35[(j-109),(1:(length(obj$counts)))] <- obj$counts
}
# View(hist_mat35)
### COHORT 5 # 6.5 OK!
##
### Now COHORT 6 # 6.6
## check:
# 1.) increase t0 by 1
# 2.)increase by 1
#for (j in (1+(4*12)):121) {
#
# N_Mort <- Abs_vector[j-(4*12)]
### replace absentees with NaN for ALL columns
# define new matrix
Lmat36 <- matrix( NA , nrow = 1000, ncol = 121)
# 6.6.a. fill with growth trajectories
Lmat36 <- matrix( NA , nrow = 1000, ncol = 121)
for (i in 1:1000) {
for (w in 0:121) {
Lmat36[i,w] <- VBG.fun.full(t_vec_10[w], sample(Linfe_vec,1),
sample(Ke_vec,1), (5+ sample(t0e_vec,1)))
}
}
# 6.6.b. clean negative L values (first columns)
Lmat36 <- ifelse(Lmat36 < 0, NA,Lmat36)
# View(Lmat36)
plot(Lmat36[12,] ~ t_vec_10 )
## 6.6.c. apply Mortality
#INPUT: mortality vector Nt_vec1
plot(Nt_vec1 ~ t_vec_10)
Nt_vec1
## [1] 1.000000e+02 8.464817e+01 7.165313e+01 6.065307e+01 5.134171e+01
## [6] 4.345982e+01 3.678794e+01 3.114032e+01 2.635971e+01 2.231302e+01
## [11] 1.888756e+01 1.598797e+01 1.353353e+01 1.145588e+01 9.697197e+00
## [16] 8.208500e+00 6.948345e+00 5.881647e+00 4.978707e+00 4.214384e+00
## [21] 3.567399e+00 3.019738e+00 2.556153e+00 2.163737e+00 1.831564e+00
## [26] 1.550385e+00 1.312373e+00 1.110900e+00 9.403563e-01 7.959944e-01
## [31] 6.737947e-01 5.703549e-01 4.827950e-01 4.086771e-01 3.459377e-01
## [36] 2.928300e-01 2.478752e-01 2.098218e-01 1.776104e-01 1.503439e-01
## [41] 1.272634e-01 1.077261e-01 9.118820e-02 7.718914e-02 6.533920e-02
## [46] 5.530844e-02 4.681758e-02 3.963023e-02 3.354626e-02 2.839630e-02
## [51] 2.403695e-02 2.034684e-02 1.722323e-02 1.457915e-02 1.234098e-02
## [56] 1.044641e-02 8.842699e-03 7.485183e-03 6.336071e-03 5.363368e-03
## [61] 4.539993e-03 3.843021e-03 3.253047e-03 2.753645e-03 2.330910e-03
## [66] 1.973073e-03 1.670170e-03 1.413768e-03 1.196729e-03 1.013009e-03
## [71] 8.574939e-04 7.258529e-04 6.144212e-04 5.200963e-04 4.402521e-04
## [76] 3.726653e-04 3.154544e-04 2.670264e-04 2.260329e-04 1.913328e-04
## [81] 1.619597e-04 1.370959e-04 1.160492e-04 9.823351e-05 8.315287e-05
## [86] 7.038739e-05 5.958164e-05 5.043477e-05 4.269211e-05 3.613809e-05
## [91] 3.059023e-05 2.589407e-05 2.191886e-05 1.855391e-05 1.570555e-05
## [96] 1.329446e-05 1.125352e-05 9.525897e-06 8.063498e-06 6.825603e-06
## [101] 5.777749e-06 4.890759e-06 4.139938e-06 3.504382e-06 2.966395e-06
## [106] 2.510999e-06 2.125515e-06 1.799210e-06 1.522998e-06 1.289190e-06
## [111] 1.091276e-06 9.237450e-07 7.819332e-07 6.618922e-07 5.602796e-07
## [116] 4.742665e-07 4.014579e-07 3.398268e-07 2.876572e-07 2.434965e-07
## [121] 2.061154e-07
Mort_as_Factor_vec <- (Nt_vec1/Nt_vec1[1])
Abs_vector <- length(Lmat36[,1]) - (Nt_vec1/Nt_vec1[1] * length(Lmat36[,1]))
Abs_vector <- round(Abs_vector,0) # vector with absentees in each month
length(Abs_vector)
## [1] 121
Abs_vector[3] # number of absentees in month 3
## [1] 283
N_Mort <- 50 # number of absentees in a month
#adjust start column or Mort.
# replace "(j in 1:121)" by "(j in (1+12):121)", # "Abs_vector[j]" by "Abs_vector[j-12]"
for (j in (1+(5*12)):121) {
N_Mort <- Abs_vector[j-(5*12)]
N_Mort_lines_vec <-sample((seq(1:length(Lmat36[,1]))), N_Mort, replace = FALSE)
# replace absentees with NaN for column j
for (y in 1: N_Mort) {
Lmat36[(N_Mort_lines_vec[y]), j] <- NA
}
}
# View(Lmat36)
#plots!
plot(Lmat36[12,] ~ t_vec_10 )
# plot N(t) of the matrix # plot N_alive_vec
sum(is.na(Lmat36[,15])) # number ofNAs in column 15
## [1] 1000
N_alive <- length(Lmat36[,16]) - sum(is.na(Lmat36[,16])) # number ofNAs in column 16
N_alive_vec <- rep(0, ncol(Lmat36))
for (j in 1: ncol(Lmat36)) {
N_alive_vec[j] <- length(Lmat36[,j]) - sum(is.na(Lmat36[,j])) # number ofNAs in column j
}
plot(N_alive_vec ~ t_vec_10 , main = "N ind. in cohort")
text(8, 800, paste("Z_input=", Z_input))
N_alive_vec[(9*12)+1] # N alive at the onset of year 10
## [1] 0
# coh 1( n= 27 ind. when Z = 0.4)
# coh 2( n= 41 ind. when Z = 0.4)
# coh 3( n= 61 ind. when Z = 0.4)
# coh 4( n= 91 ind. when Z = 0.4)
# coh 5( n= 135 ind. when Z = 0.4)
# coh 6( n= 202 ind. when Z = 0.4)
# percent Coh 1 at the onset of year 10
Ncoh1_ini10 <- 27
NTOTAL_ini10 <-sum (27,41,61,91,135, 202)
NTOTAL_ini10
## [1] 557
Ncoh1_ini10 / NTOTAL_ini10
## [1] 0.04847397
Ncoh1_ini10 / NTOTAL_ini10 * 100 # COH1 in percent
## [1] 4.847397
# 6.6.d. make 12 histograms, make a catch curve for the last year
# vector for "mids"
mids_OK_05 <- seq(0, 15, by = 0.5)
breaks_OK_05 <- seq(-0.2, 15.2, by = 0.5)
hist_mat36 <- matrix(0, nrow = (12), ncol = (length(breaks_OK_05)-1))
for (j in 110:121) {
hist(Lmat36[,j] , xlim = c(0, 15), ylim = c(0, 500), breaks = breaks_OK_05 )
obj <- hist(Lmat36[,j] , xlim = c(0, 15), ylim = c(0, 500), breaks = breaks_OK_05 )
hist_mat36[(j-109),(1:(length(obj$counts)))] <- obj$counts
}
# View(hist_mat36)
### COHORT 6 # 6.6 OK!
###
# NOW COHORT 7 # 6.7
## check:
# 1.) increase t0 by 1
# 2.)increase by 1
#for (j in (1+(4*12)):121) {
#
# N_Mort <- Abs_vector[j-(4*12)]
### replace absentees with NaN for ALL columns
# define new matrix
Lmat37 <- matrix( NA , nrow = 1000, ncol = 121)
# 6.7.a. fill with growth trajectories
Lmat37 <- matrix( NA , nrow = 1000, ncol = 121)
for (i in 1:1000) {
for (w in 0:121) {
Lmat37[i,w] <- VBG.fun.full(t_vec_10[w], sample(Linfe_vec,1),
sample(Ke_vec,1), (6+ sample(t0e_vec,1)))
}
}
# 6.7.b. clean negative L values (first columns)
Lmat37 <- ifelse(Lmat37 < 0, NA,Lmat37)
# View(Lmat37)
plot(Lmat37[12,] ~ t_vec_10 )
## 6.7.c. apply Mortality
#INPUT: mortality vector Nt_vec1
plot(Nt_vec1 ~ t_vec_10)
Nt_vec1
## [1] 1.000000e+02 8.464817e+01 7.165313e+01 6.065307e+01 5.134171e+01
## [6] 4.345982e+01 3.678794e+01 3.114032e+01 2.635971e+01 2.231302e+01
## [11] 1.888756e+01 1.598797e+01 1.353353e+01 1.145588e+01 9.697197e+00
## [16] 8.208500e+00 6.948345e+00 5.881647e+00 4.978707e+00 4.214384e+00
## [21] 3.567399e+00 3.019738e+00 2.556153e+00 2.163737e+00 1.831564e+00
## [26] 1.550385e+00 1.312373e+00 1.110900e+00 9.403563e-01 7.959944e-01
## [31] 6.737947e-01 5.703549e-01 4.827950e-01 4.086771e-01 3.459377e-01
## [36] 2.928300e-01 2.478752e-01 2.098218e-01 1.776104e-01 1.503439e-01
## [41] 1.272634e-01 1.077261e-01 9.118820e-02 7.718914e-02 6.533920e-02
## [46] 5.530844e-02 4.681758e-02 3.963023e-02 3.354626e-02 2.839630e-02
## [51] 2.403695e-02 2.034684e-02 1.722323e-02 1.457915e-02 1.234098e-02
## [56] 1.044641e-02 8.842699e-03 7.485183e-03 6.336071e-03 5.363368e-03
## [61] 4.539993e-03 3.843021e-03 3.253047e-03 2.753645e-03 2.330910e-03
## [66] 1.973073e-03 1.670170e-03 1.413768e-03 1.196729e-03 1.013009e-03
## [71] 8.574939e-04 7.258529e-04 6.144212e-04 5.200963e-04 4.402521e-04
## [76] 3.726653e-04 3.154544e-04 2.670264e-04 2.260329e-04 1.913328e-04
## [81] 1.619597e-04 1.370959e-04 1.160492e-04 9.823351e-05 8.315287e-05
## [86] 7.038739e-05 5.958164e-05 5.043477e-05 4.269211e-05 3.613809e-05
## [91] 3.059023e-05 2.589407e-05 2.191886e-05 1.855391e-05 1.570555e-05
## [96] 1.329446e-05 1.125352e-05 9.525897e-06 8.063498e-06 6.825603e-06
## [101] 5.777749e-06 4.890759e-06 4.139938e-06 3.504382e-06 2.966395e-06
## [106] 2.510999e-06 2.125515e-06 1.799210e-06 1.522998e-06 1.289190e-06
## [111] 1.091276e-06 9.237450e-07 7.819332e-07 6.618922e-07 5.602796e-07
## [116] 4.742665e-07 4.014579e-07 3.398268e-07 2.876572e-07 2.434965e-07
## [121] 2.061154e-07
Mort_as_Factor_vec <- (Nt_vec1/Nt_vec1[1])
Abs_vector <- length(Lmat37[,1]) - (Nt_vec1/Nt_vec1[1] * length(Lmat37[,1]))
Abs_vector <- round(Abs_vector,0) # vector with absentees in each month
length(Abs_vector)
## [1] 121
Abs_vector[3] # number of absentees in month 3
## [1] 283
N_Mort <- 50 # number of absentees in a month
#adjust start column or Mort.
# replace "(j in 1:121)" by "(j in (1+12):121)", # "Abs_vector[j]" by "Abs_vector[j-12]"
for (j in (1+(6*12)):121) {
N_Mort <- Abs_vector[j-(6*12)]
N_Mort_lines_vec <-sample((seq(1:length(Lmat37[,1]))), N_Mort, replace = FALSE)
# replace absentees with NaN for column j
for (y in 1: N_Mort) {
Lmat37[(N_Mort_lines_vec[y]), j] <- NA
}
}
# View(Lmat37)
#plots!
plot(Lmat37[12,] ~ t_vec_10 )
# plot N(t) of the matrix # plot N_alive_vec
sum(is.na(Lmat37[,15])) # number ofNAs in column 15
## [1] 1000
N_alive <- length(Lmat37[,16]) - sum(is.na(Lmat37[,16])) # number ofNAs in column 16
N_alive_vec <- rep(0, ncol(Lmat37))
for (j in 1: ncol(Lmat37)) {
N_alive_vec[j] <- length(Lmat37[,j]) - sum(is.na(Lmat37[,j])) # number ofNAs in column j
}
plot(N_alive_vec ~ t_vec_10 , main = "N ind. in cohort")
text(8, 800, paste("Z_input=", Z_input))
N_alive_vec[(9*12)+1] # N alive at the onset of year 10
## [1] 2
# coh 1( n= 27 ind. when Z = 0.4)
# coh 2( n= 41 ind. when Z = 0.4)
# coh 3( n= 61 ind. when Z = 0.4)
# coh 4( n= 91 ind. when Z = 0.4)
# coh 5( n= 135 ind. when Z = 0.4)
# coh 6( n= 202 ind. when Z = 0.4)
# coh 7( n= 301 ind. when Z = 0.4)
# percent Coh 1 at the onset of year 10
Ncoh1_ini10 <- 27
NTOTAL_ini10 <-sum (27,41,61,91,135, 202, 301)
NTOTAL_ini10
## [1] 858
Ncoh1_ini10 / NTOTAL_ini10
## [1] 0.03146853
Ncoh1_ini10 / NTOTAL_ini10 * 100 # COH1 in percent
## [1] 3.146853
# 6.7.d. make 12 histograms, make a catch curve for the last year
# vector for "mids"
mids_OK_05 <- seq(0, 15, by = 0.5)
breaks_OK_05 <- seq(-0.2, 15.2, by = 0.5)
hist_mat37 <- matrix(0, nrow = (12), ncol = (length(breaks_OK_05)-1))
for (j in 110:121) {
hist(Lmat37[,j] , xlim = c(0, 15), ylim = c(0, 500), breaks = breaks_OK_05 )
obj <- hist(Lmat37[,j] , xlim = c(0, 15), ylim = c(0, 500), breaks = breaks_OK_05 )
hist_mat37[(j-109),(1:(length(obj$counts)))] <- obj$counts
}
# View(hist_mat37)
### COHORT 7 # 6.7 OK!
###
# Now COHORT 8
##### # 6.8
## check:
# 1.) increase t0 by 1
# 2.)increase by 1
#for (j in (1+(4*12)):121) {
#
# N_Mort <- Abs_vector[j-(4*12)]
### replace absentees with NaN for ALL columns
# define new matrix
Lmat38 <- matrix( NA , nrow = 1000, ncol = 121)
# 6.8.a. fill with growth trajectories
Lmat38 <- matrix( NA , nrow = 1000, ncol = 121)
for (i in 1:1000) {
for (w in 0:121) {
Lmat38[i,w] <- VBG.fun.full(t_vec_10[w], sample(Linfe_vec,1),
sample(Ke_vec,1), (7+ sample(t0e_vec,1)))
}
}
# 6.8.b. clean negative L values (first columns)
Lmat38 <- ifelse(Lmat38 < 0, NA,Lmat38)
# View(Lmat38)
plot(Lmat38[12,] ~ t_vec_10 )
## 6.8.c. apply Mortality
#INPUT: mortality vector Nt_vec1
plot(Nt_vec1 ~ t_vec_10)
Nt_vec1
## [1] 1.000000e+02 8.464817e+01 7.165313e+01 6.065307e+01 5.134171e+01
## [6] 4.345982e+01 3.678794e+01 3.114032e+01 2.635971e+01 2.231302e+01
## [11] 1.888756e+01 1.598797e+01 1.353353e+01 1.145588e+01 9.697197e+00
## [16] 8.208500e+00 6.948345e+00 5.881647e+00 4.978707e+00 4.214384e+00
## [21] 3.567399e+00 3.019738e+00 2.556153e+00 2.163737e+00 1.831564e+00
## [26] 1.550385e+00 1.312373e+00 1.110900e+00 9.403563e-01 7.959944e-01
## [31] 6.737947e-01 5.703549e-01 4.827950e-01 4.086771e-01 3.459377e-01
## [36] 2.928300e-01 2.478752e-01 2.098218e-01 1.776104e-01 1.503439e-01
## [41] 1.272634e-01 1.077261e-01 9.118820e-02 7.718914e-02 6.533920e-02
## [46] 5.530844e-02 4.681758e-02 3.963023e-02 3.354626e-02 2.839630e-02
## [51] 2.403695e-02 2.034684e-02 1.722323e-02 1.457915e-02 1.234098e-02
## [56] 1.044641e-02 8.842699e-03 7.485183e-03 6.336071e-03 5.363368e-03
## [61] 4.539993e-03 3.843021e-03 3.253047e-03 2.753645e-03 2.330910e-03
## [66] 1.973073e-03 1.670170e-03 1.413768e-03 1.196729e-03 1.013009e-03
## [71] 8.574939e-04 7.258529e-04 6.144212e-04 5.200963e-04 4.402521e-04
## [76] 3.726653e-04 3.154544e-04 2.670264e-04 2.260329e-04 1.913328e-04
## [81] 1.619597e-04 1.370959e-04 1.160492e-04 9.823351e-05 8.315287e-05
## [86] 7.038739e-05 5.958164e-05 5.043477e-05 4.269211e-05 3.613809e-05
## [91] 3.059023e-05 2.589407e-05 2.191886e-05 1.855391e-05 1.570555e-05
## [96] 1.329446e-05 1.125352e-05 9.525897e-06 8.063498e-06 6.825603e-06
## [101] 5.777749e-06 4.890759e-06 4.139938e-06 3.504382e-06 2.966395e-06
## [106] 2.510999e-06 2.125515e-06 1.799210e-06 1.522998e-06 1.289190e-06
## [111] 1.091276e-06 9.237450e-07 7.819332e-07 6.618922e-07 5.602796e-07
## [116] 4.742665e-07 4.014579e-07 3.398268e-07 2.876572e-07 2.434965e-07
## [121] 2.061154e-07
Mort_as_Factor_vec <- (Nt_vec1/Nt_vec1[1])
Abs_vector <- length(Lmat38[,1]) - (Nt_vec1/Nt_vec1[1] * length(Lmat38[,1]))
Abs_vector <- round(Abs_vector,0) # vector with absentees in each month
length(Abs_vector)
## [1] 121
Abs_vector[3] # number of absentees in month 3
## [1] 283
N_Mort <- 50 # number of absentees in a month
#adjust start column or Mort.
# replace "(j in 1:121)" by "(j in (1+12):121)", # "Abs_vector[j]" by "Abs_vector[j-12]"
for (j in (1+(7*12)):121) {
N_Mort <- Abs_vector[j-(7*12)]
N_Mort_lines_vec <-sample((seq(1:length(Lmat38[,1]))), N_Mort, replace = FALSE)
# replace absentees with NaN for column j
for (y in 1: N_Mort) {
Lmat38[(N_Mort_lines_vec[y]), j] <- NA
}
}
# View(Lmat38)
#plots!
plot(Lmat38[12,] ~ t_vec_10 )
# plot N(t) of the matrix # plot N_alive_vec
sum(is.na(Lmat38[,15])) # number ofNAs in column 15
## [1] 1000
N_alive <- length(Lmat38[,16]) - sum(is.na(Lmat38[,16])) # number ofNAs in column 16
N_alive_vec <- rep(0, ncol(Lmat38))
for (j in 1: ncol(Lmat38)) {
N_alive_vec[j] <- length(Lmat38[,j]) - sum(is.na(Lmat38[,j])) # number ofNAs in column j
}
plot(N_alive_vec ~ t_vec_10 , main = "N ind. in cohort")
text(8, 800, paste("Z_input=", Z_input))
N_alive_vec[(9*12)+1] # N alive at the onset of year 10
## [1] 18
# coh 1( n= 27 ind. when Z = 0.4)
# coh 2( n= 41 ind. when Z = 0.4)
# coh 3( n= 61 ind. when Z = 0.4)
# coh 4( n= 91 ind. when Z = 0.4)
# coh 5( n= 135 ind. when Z = 0.4)
# coh 6( n= 202 ind. when Z = 0.4)
# coh 7( n= 301 ind. when Z = 0.4)
# coh 8( n= 449 ind. when Z = 0.4)
# percent Coh 1 at the onset of year 10
Ncoh1_ini10 <- 27
NTOTAL_ini10 <-sum (27,41,61,91,135, 202, 301,449)
NTOTAL_ini10
## [1] 1307
Ncoh1_ini10 / NTOTAL_ini10
## [1] 0.020658
Ncoh1_ini10 / NTOTAL_ini10 * 100 # COH1 in percent
## [1] 2.0658
# 6.8.d. make 12 histograms, make a catch curve for the last year
# vector for "mids"
mids_OK_05 <- seq(0, 15, by = 0.5)
breaks_OK_05 <- seq(-0.2, 15.2, by = 0.5)
hist_mat38 <- matrix(0, nrow = (12), ncol = (length(breaks_OK_05)-1))
for (j in 110:121) {
hist(Lmat38[,j] , xlim = c(0, 15), ylim = c(0, 500), breaks = breaks_OK_05 )
obj <- hist(Lmat38[,j] , xlim = c(0, 15), ylim = c(0, 500), breaks = breaks_OK_05 )
hist_mat38[(j-109),(1:(length(obj$counts)))] <- obj$counts
}
# View(hist_mat38)
### COHORT 8 # 6.8 OK!
#####
### NOW COHORT 9 # 6.9
##
## check:
# 1.) increase t0 by 1
# 2.)increase by 1
#for (j in (1+(4*12)):121) {
#
# N_Mort <- Abs_vector[j-(4*12)]
### replace absentees with NaN for ALL columns
# define new matrix
Lmat39 <- matrix( NA , nrow = 1000, ncol = 121)
# 6.9.a. fill with growth trajectories
Lmat39 <- matrix( NA , nrow = 1000, ncol = 121)
for (i in 1:1000) {
for (w in 0:121) {
Lmat39[i,w] <- VBG.fun.full(t_vec_10[w], sample(Linfe_vec,1),
sample(Ke_vec,1), (8+ sample(t0e_vec,1)))
}
}
# 6.9.b. clean negative L values (first columns)
Lmat39 <- ifelse(Lmat39 < 0, NA,Lmat39)
# View(Lmat39)
plot(Lmat39[12,] ~ t_vec_10 )
## 6.9.c. apply Mortality
#INPUT: mortality vector Nt_vec1
plot(Nt_vec1 ~ t_vec_10)
Nt_vec1
## [1] 1.000000e+02 8.464817e+01 7.165313e+01 6.065307e+01 5.134171e+01
## [6] 4.345982e+01 3.678794e+01 3.114032e+01 2.635971e+01 2.231302e+01
## [11] 1.888756e+01 1.598797e+01 1.353353e+01 1.145588e+01 9.697197e+00
## [16] 8.208500e+00 6.948345e+00 5.881647e+00 4.978707e+00 4.214384e+00
## [21] 3.567399e+00 3.019738e+00 2.556153e+00 2.163737e+00 1.831564e+00
## [26] 1.550385e+00 1.312373e+00 1.110900e+00 9.403563e-01 7.959944e-01
## [31] 6.737947e-01 5.703549e-01 4.827950e-01 4.086771e-01 3.459377e-01
## [36] 2.928300e-01 2.478752e-01 2.098218e-01 1.776104e-01 1.503439e-01
## [41] 1.272634e-01 1.077261e-01 9.118820e-02 7.718914e-02 6.533920e-02
## [46] 5.530844e-02 4.681758e-02 3.963023e-02 3.354626e-02 2.839630e-02
## [51] 2.403695e-02 2.034684e-02 1.722323e-02 1.457915e-02 1.234098e-02
## [56] 1.044641e-02 8.842699e-03 7.485183e-03 6.336071e-03 5.363368e-03
## [61] 4.539993e-03 3.843021e-03 3.253047e-03 2.753645e-03 2.330910e-03
## [66] 1.973073e-03 1.670170e-03 1.413768e-03 1.196729e-03 1.013009e-03
## [71] 8.574939e-04 7.258529e-04 6.144212e-04 5.200963e-04 4.402521e-04
## [76] 3.726653e-04 3.154544e-04 2.670264e-04 2.260329e-04 1.913328e-04
## [81] 1.619597e-04 1.370959e-04 1.160492e-04 9.823351e-05 8.315287e-05
## [86] 7.038739e-05 5.958164e-05 5.043477e-05 4.269211e-05 3.613809e-05
## [91] 3.059023e-05 2.589407e-05 2.191886e-05 1.855391e-05 1.570555e-05
## [96] 1.329446e-05 1.125352e-05 9.525897e-06 8.063498e-06 6.825603e-06
## [101] 5.777749e-06 4.890759e-06 4.139938e-06 3.504382e-06 2.966395e-06
## [106] 2.510999e-06 2.125515e-06 1.799210e-06 1.522998e-06 1.289190e-06
## [111] 1.091276e-06 9.237450e-07 7.819332e-07 6.618922e-07 5.602796e-07
## [116] 4.742665e-07 4.014579e-07 3.398268e-07 2.876572e-07 2.434965e-07
## [121] 2.061154e-07
Mort_as_Factor_vec <- (Nt_vec1/Nt_vec1[1])
Abs_vector <- length(Lmat39[,1]) - (Nt_vec1/Nt_vec1[1] * length(Lmat39[,1]))
Abs_vector <- round(Abs_vector,0) # vector with absentees in each month
length(Abs_vector)
## [1] 121
Abs_vector[3] # number of absentees in month 3
## [1] 283
N_Mort <- 50 # number of absentees in a month
#adjust start column or Mort.
# replace "(j in 1:121)" by "(j in (1+12):121)", # "Abs_vector[j]" by "Abs_vector[j-12]"
for (j in (1+(8*12)):121) {
N_Mort <- Abs_vector[j-(8*12)]
N_Mort_lines_vec <-sample((seq(1:length(Lmat39[,1]))), N_Mort, replace = FALSE)
# replace absentees with NaN for column j
for (y in 1: N_Mort) {
Lmat39[(N_Mort_lines_vec[y]), j] <- NA
}
}
# View(Lmat39)
#plots!
plot(Lmat39[12,] ~ t_vec_10 )
# plot N(t) of the matrix # plot N_alive_vec
sum(is.na(Lmat39[,15])) # number ofNAs in column 15
## [1] 1000
N_alive <- length(Lmat39[,16]) - sum(is.na(Lmat39[,16])) # number ofNAs in column 16
N_alive_vec <- rep(0, ncol(Lmat39))
for (j in 1: ncol(Lmat39)) {
N_alive_vec[j] <- length(Lmat39[,j]) - sum(is.na(Lmat39[,j])) # number ofNAs in column j
}
plot(N_alive_vec ~ t_vec_10 , main = "N ind. in cohort")
text(2, 800, paste("Z_input=", Z_input))
N_alive_vec[(9*12)+1] # N alive at the onset of year 10
## [1] 135
# coh 1( n= 27 ind. when Z = 0.4)
# coh 2( n= 41 ind. when Z = 0.4)
# coh 3( n= 61 ind. when Z = 0.4)
# coh 4( n= 91 ind. when Z = 0.4)
# coh 5( n= 135 ind. when Z = 0.4)
# coh 6( n= 202 ind. when Z = 0.4)
# coh 7( n= 301 ind. when Z = 0.4)
# coh 8( n= 449 ind. when Z = 0.4)
# coh 8( n= 670 ind. when Z = 0.4)
# percent Coh 1 at the onset of year 10
Ncoh1_ini10 <- 27
NTOTAL_ini10 <-sum (27,41,61,91,135, 202, 301,449,670)
NTOTAL_ini10
## [1] 1977
Ncoh1_ini10 / NTOTAL_ini10
## [1] 0.01365706
Ncoh1_ini10 / NTOTAL_ini10 * 100 # COH1 in percent
## [1] 1.365706
# 6.9.d. make 12 histograms, make a catch curve for the last year
# vector for "mids"
mids_OK_05 <- seq(0, 15, by = 0.5)
breaks_OK_05 <- seq(-0.2, 15.2, by = 0.5)
hist_mat39 <- matrix(0, nrow = (12), ncol = (length(breaks_OK_05)-1))
for (j in 110:121) {
hist(Lmat39[,j] , xlim = c(0, 15), ylim = c(0, 500), breaks = breaks_OK_05 )
obj <- hist(Lmat39[,j] , xlim = c(0, 15), ylim = c(0, 500), breaks = breaks_OK_05 )
hist_mat39[(j-109),(1:(length(obj$counts)))] <- obj$counts
}
# View(hist_mat39)
### COHORT 9 # 6.9 OK!
### ###
### NOW COHORT 10 # 6.10
### ###
## check:
# 1.) increase t0 by 1
# 2.)increase by 1
#for (j in (1+(4*12)):121) {
#
# N_Mort <- Abs_vector[j-(4*12)]
### replace absentees with NaN for ALL columns
# define new matrix
Lmat310 <- matrix( NA , nrow = 1000, ncol = 121)
# 6.10.a. fill with growth trajectories
Lmat310 <- matrix( NA , nrow = 1000, ncol = 121)
for (i in 1:1000) {
for (w in 0:121) {
Lmat310[i,w] <- VBG.fun.full(t_vec_10[w], sample(Linfe_vec,1),
sample(Ke_vec,1), (9+ sample(t0e_vec,1)))
}
}
# 6.10.b. clean negative L values (first columns)
Lmat310 <- ifelse(Lmat310 < 0, NA,Lmat310)
# View(Lmat310)
plot(Lmat310[12,] ~ t_vec_10 )
## 6.10.c. apply Mortality
#INPUT: mortality vector Nt_vec1
plot(Nt_vec1 ~ t_vec_10)
Nt_vec1
## [1] 1.000000e+02 8.464817e+01 7.165313e+01 6.065307e+01 5.134171e+01
## [6] 4.345982e+01 3.678794e+01 3.114032e+01 2.635971e+01 2.231302e+01
## [11] 1.888756e+01 1.598797e+01 1.353353e+01 1.145588e+01 9.697197e+00
## [16] 8.208500e+00 6.948345e+00 5.881647e+00 4.978707e+00 4.214384e+00
## [21] 3.567399e+00 3.019738e+00 2.556153e+00 2.163737e+00 1.831564e+00
## [26] 1.550385e+00 1.312373e+00 1.110900e+00 9.403563e-01 7.959944e-01
## [31] 6.737947e-01 5.703549e-01 4.827950e-01 4.086771e-01 3.459377e-01
## [36] 2.928300e-01 2.478752e-01 2.098218e-01 1.776104e-01 1.503439e-01
## [41] 1.272634e-01 1.077261e-01 9.118820e-02 7.718914e-02 6.533920e-02
## [46] 5.530844e-02 4.681758e-02 3.963023e-02 3.354626e-02 2.839630e-02
## [51] 2.403695e-02 2.034684e-02 1.722323e-02 1.457915e-02 1.234098e-02
## [56] 1.044641e-02 8.842699e-03 7.485183e-03 6.336071e-03 5.363368e-03
## [61] 4.539993e-03 3.843021e-03 3.253047e-03 2.753645e-03 2.330910e-03
## [66] 1.973073e-03 1.670170e-03 1.413768e-03 1.196729e-03 1.013009e-03
## [71] 8.574939e-04 7.258529e-04 6.144212e-04 5.200963e-04 4.402521e-04
## [76] 3.726653e-04 3.154544e-04 2.670264e-04 2.260329e-04 1.913328e-04
## [81] 1.619597e-04 1.370959e-04 1.160492e-04 9.823351e-05 8.315287e-05
## [86] 7.038739e-05 5.958164e-05 5.043477e-05 4.269211e-05 3.613809e-05
## [91] 3.059023e-05 2.589407e-05 2.191886e-05 1.855391e-05 1.570555e-05
## [96] 1.329446e-05 1.125352e-05 9.525897e-06 8.063498e-06 6.825603e-06
## [101] 5.777749e-06 4.890759e-06 4.139938e-06 3.504382e-06 2.966395e-06
## [106] 2.510999e-06 2.125515e-06 1.799210e-06 1.522998e-06 1.289190e-06
## [111] 1.091276e-06 9.237450e-07 7.819332e-07 6.618922e-07 5.602796e-07
## [116] 4.742665e-07 4.014579e-07 3.398268e-07 2.876572e-07 2.434965e-07
## [121] 2.061154e-07
Mort_as_Factor_vec <- (Nt_vec1/Nt_vec1[1])
Abs_vector <- length(Lmat310[,1]) - (Nt_vec1/Nt_vec1[1] * length(Lmat310[,1]))
Abs_vector <- round(Abs_vector,0) # vector with absentees in each month
length(Abs_vector)
## [1] 121
Abs_vector[3] # number of absentees in month 3
## [1] 283
N_Mort <- 50 # number of absentees in a month
#adjust start column or Mort.
# replace "(j in 1:121)" by "(j in (1+12):121)", # "Abs_vector[j]" by "Abs_vector[j-12]"
for (j in (1+(9*12)):121) {
N_Mort <- Abs_vector[j-(9*12)]
N_Mort_lines_vec <-sample((seq(1:length(Lmat310[,1]))), N_Mort, replace = FALSE)
# replace absentees with NaN for column j
for (y in 1: N_Mort) {
Lmat310[(N_Mort_lines_vec[y]), j] <- NA
}
}
# View(Lmat310)
#plots!
plot(Lmat310[12,] ~ t_vec_10 )
# plot N(t) of the matrix # plot N_alive_vec
sum(is.na(Lmat310[,15])) # number ofNAs in column 15
## [1] 1000
N_alive <- length(Lmat310[,16]) - sum(is.na(Lmat310[,16])) # number ofNAs in column 16
N_alive_vec <- rep(0, ncol(Lmat310))
for (j in 1: ncol(Lmat310)) {
N_alive_vec[j] <- length(Lmat310[,j]) - sum(is.na(Lmat310[,j])) # number ofNAs in column j
}
plot(N_alive_vec ~ t_vec_10 , main = "N ind. in cohort")
text(2, 800, paste("Z_input=", Z_input))
N_alive_vec[(9*12)+1] # N alive at the onset of year 10
## [1] 712
# coh 1( n= 27 ind. when Z = 0.4)
# coh 2( n= 41 ind. when Z = 0.4)
# coh 3( n= 61 ind. when Z = 0.4)
# coh 4( n= 91 ind. when Z = 0.4)
# coh 5( n= 135 ind. when Z = 0.4)
# coh 6( n= 202 ind. when Z = 0.4)
# coh 7( n= 301 ind. when Z = 0.4)
# coh 8( n= 449 ind. when Z = 0.4)
# coh 9( n= 670 ind. when Z = 0.4)
# coh 9( n= 1000 ind. when Z = 0.4)
# percent Coh 1 at the onset of year 10
Ncoh1_ini10 <- 27
NTOTAL_ini10 <-sum (27,41,61,91,135, 202, 301,449,670, 1000)
NTOTAL_ini10
## [1] 2977
Ncoh1_ini10 / NTOTAL_ini10
## [1] 0.009069533
Ncoh1_ini10 / NTOTAL_ini10 * 100 # COH1 in percent
## [1] 0.9069533
# 6.10.d. make 12 histograms, make a catch curve for the last year
# vector for "mids"
mids_OK_05 <- seq(0, 15, by = 0.5)
breaks_OK_05 <- seq(-0.2, 15.2, by = 0.5)
hist_mat310 <- matrix(0, nrow = (12), ncol = (length(breaks_OK_05)-1))
for (j in 110:121) {
hist(Lmat310[,j] , xlim = c(0, 15), ylim = c(0, 500), breaks = breaks_OK_05 )
obj <- hist(Lmat310[,j] , xlim = c(0, 15), ylim = c(0, 500), breaks = breaks_OK_05 )
hist_mat310[(j-109),(1:(length(obj$counts)))] <- obj$counts
}
# View(hist_mat310)
### COHORT 10 # 6.10 OK!
# 7. summarize (integrate)
# 7.0 Plot ALL cohorts together
plot( Lmat3[1,1:121] ~ t_vec_10,
col=rgb(0, 0.1, 1, 0), pch=1, cex = 0.4,
ylim = c(0, 15), xlab ="t (years)", ylab = "Length(cm)", main = "N = 10,000")
for (i in 1:1000) {
points( Lmat3[i,1:121] ~ t_vec_10, col=rgb(0, 0.1, 1, 0.04), pch=1, cex = 0.5)
}
for (i in 1:1000) {
points( Lmat32[i,1:121] ~ t_vec_10, col=rgb(0, 0.1, 1, 0.04), pch=1, cex = 0.5)
}
for (i in 1:1000) {
points( Lmat33[i,1:121] ~ t_vec_10, col=rgb(0, 0.1, 1, 0.04), pch=1, cex = 0.5)
}
for (i in 1:1000) {
points( Lmat34[i,1:121] ~ t_vec_10, col=rgb(0, 0.1, 1, 0.04), pch=1, cex = 0.5)
}
for (i in 1:1000) {
points( Lmat35[i,1:121] ~ t_vec_10, col=rgb(0, 0.1, 1, 0.04), pch=1, cex = 0.5)
}
for (i in 1:1000) {
points( Lmat36[i,1:121] ~ t_vec_10, col=rgb(0, 0.1, 1, 0.04), pch=1, cex = 0.5)
}
for (i in 1:1000) {
points( Lmat37[i,1:121] ~ t_vec_10, col=rgb(0, 0.1, 1, 0.04), pch=1, cex = 0.5)
}
for (i in 1:1000) {
points( Lmat38[i,1:121] ~ t_vec_10, col=rgb(0, 0.1, 1, 0.04), pch=1, cex = 0.5)
}
for (i in 1:1000) {
points( Lmat39[i,1:121] ~ t_vec_10, col=rgb(0, 0.1, 1, 0.04), pch=1, cex = 0.5)
}
for (i in 1:1000) {
points( Lmat310[i,1:121] ~ t_vec_10, col=rgb(0, 0.1, 1, 0.04), pch=1, cex = 0.5)
}
text(1.6, 14, paste("Z:",Z_input,"K:", K_input), cex= 0.8)
text(7, 14, paste("C.V.Linf:",cvLinf_input,"C.V.K:", cvK_input, "CVt0:", cvt0_input ),cex= 0.8)
abline (h = 10, col=rgb(0, 0.7, 0.23, 0.3) , lwd=0.7,lty="dotted")
text(0.2, 11, paste("Linf"), cex= 0.8, col = "grey")
# 7.A sum of 10 matrices
hmat_list <- list(hist_mat3,hist_mat32,hist_mat33, hist_mat34, hist_mat35,
hist_mat36, hist_mat37,hist_mat38, hist_mat39, hist_mat310)
hist_mat_all <- Reduce('+', hmat_list)
hist(hist_mat_all[1,(1:30)] ) # [row, columns}
hist(hist_mat_all[4,(1:30)] ) # [row, columns}
hist(hist_mat_all[6,(1:30)] ) # [row, columns}
hist(hist_mat_all[8,(1:30)] ) # [row, columns}
hist(hist_mat_all[10,(1:30)] ) # [row, columns}
hist(hist_mat_all[12,(1:30)] ) # [row, columns}
# 7.B summarize all into one catch curve
colSums (hist_mat_all)
## [1] 77 164 205 262 350 393 425 465 401 394 429 326 338 259 214 181 145
## [18] 109 85 39 14 3 0 0 0 0 0 0 0 0
Catch_all <- colSums (hist_mat_all)
Catch_all <- Catch_all[1:31]
Len_all <- mids_OK_05
length(Catch_all)
## [1] 31
length(Len_all)
## [1] 31
plot(Catch_all~ Len_all)
# View(hist_mat_all)
# 7.C Apply Recruitment corve
# Recruitment + gear selection: logistic curve with midpont at L = 2 cm
recruit <- SSfpl(Len_all, 0, 1, 2,0.5) # recruitment + gear selection: logistic curve with midpont at L = 2 cm
length(recruit)
## [1] 31
plot(recruit ~ Len_all)
plot(recruit ~ Len_all, xlim = c(0,10), main= "Recruitment & Gear Selection", ylab= "Recruitment (0 to 1)", xlab = "L (cm)")
legend("bottomright", legend = "midpont=2 scale=0.5")
Catch_allR <- Catch_all[1:length(Len_all)]* recruit
length(Catch_allR)
## [1] 31
length(Len_all)
## [1] 31
plot(Catch_allR~ Len_all)
text(12, 300, paste("Z:",Z_input,"K:", K_input), cex= 0.8)
text(12, 200, paste("C.V.Linf:",cvLinf_input,"C.V.K:", cvK_input),cex= 0.8)
# 7. clean cut catch curve (no zeros)
Len_all <- Len_all[which(Catch_allR!=0,arr.ind = TRUE)]
Catch_allR <- Catch_allR[which(Catch_allR!=0,arr.ind = TRUE)]
length(Catch_allR)
## [1] 22
length(Len_all)
## [1] 22
plot(Catch_allR~ Len_all)
text(8.7, 500, paste("Z:",Z_input,"K:", K_input), cex= 0.8)
text(8.7, 400, paste("C.V.Linf:",cvLinf_input,"C.V.K:", cvK_input),cex= 0.8)
which(Catch_allR!=0,arr.ind = TRUE)
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22
# 8. Test with P-W methods
###
######## FUNCTION "mP_W.ffun"
#### II.a, FUNCTION mP-W plot as Function "mP_W.ffun"
## "modifed" version of the method ("mP-W plot", Pauly, 1986, Wetherall, 1986),
# with GAMMA selection mode+10% to 93%Lmax
mP_W.ffun <- function(Len, Cat) {
N_c <- length(Len) # N of length classes
int_Len <- (max(Len)-min(Len))/ (N_c-1) # interval (class width)
L_c <- Len - (int_Len/2) # cutoff lengths (L - 0.5)
# L_c <- Len - 0.5
# plot(Len, Cat)
### Calculate Mean lengths
#Catch*L
Catch_L <- Cat * Len
#Sum of "Catch*L" above each Lc
#write loop
Sum_of_Catch_times_L <- rep(0, N_c)
for(i in 1:N_c){
Sum_of_Catch_times_L[i] <- sum(Catch_L[i:N_c]);
}
#Sum of Catch above Lc
#write loop
Sum_of_Catch <- rep(0, N_c)
for(i in 1:N_c){
Sum_of_Catch[i] <- sum(Cat[i:N_c]);
}
# Calculate mean length for each Lc
mean_length <- Sum_of_Catch_times_L / Sum_of_Catch
#Pauly's modified P-W_plot method
# mP-W plot
# GAMMA selection: select data from Mode+20% to 90%Lmax
mean_length # mean length
L_c # cutoff length
df_1 <- data.frame(mean_length, L_c, Len, Cat)
d_mode <- df_1$Len[df_1$Cat == max(df_1$Cat)] # mode of original LFD data
df_gamma <- subset(df_1, Len > (1.1 * d_mode))
df_gamma <- subset(df_gamma, Len < (0.93 * max(Len)))
mean_length_G <- df_gamma$mean_length # Gamma selected mean lengths
L_c_G <- df_gamma$L_c # Gamma selected cutoff lengths
y_p <- (mean_length - L_c)
y_p_G <- (mean_length_G - L_c_G)
plot((mean_length - L_c) ~ L_c)
# plot((mean_length - L_c) ~ L_c, xlim = c(min(L_c), max(L_c)), ylim = c(min(y_p), max(y_p)))
points((mean_length_G - L_c_G) ~ L_c_G, pch = 16, col = "blue")
mP_W_mod1 <- lm((mean_length_G - L_c_G) ~ L_c_G)
abline(mP_W_mod1, col = "blue")
a.mP_W_mod1 <- summary(mP_W_mod1)$coefficients[1:1]
b.mP_W_mod1 <- summary(mP_W_mod1)$coefficients[2,1]
# Linf = a / -b and Z/K = (1 + b)/(- b)
L_inf_mP_W <- a.mP_W_mod1 / (-1 * b.mP_W_mod1)
Z_K_mP_W = (1 + b.mP_W_mod1)/(-1 * b.mP_W_mod1)
; c("LinF:", L_inf_mP_W, "Z/K:", Z_K_mP_W, "mP-W method") # output
} #### EOF # end of Function "mP_W.ffun"
###
# Z/K (INPUT is)
Z_input
## [1] 2
K_input
## [1] 1
(ZK_input <- Z_input / K_input)
## [1] 2
# OUTPUT mP-W plot an estimates of Z/K and Linf
mP_W.ffun(Len_all, Catch_allR)
## [1] "LinF:" "10.488269242757" "Z/K:"
## [4] "2.23452832911765" "mP-W method"
mPW.out1<- mP_W.ffun(Len_all, Catch_allR)
Linf_estim <- as.numeric(mPW.out1[2])
ZK_estim <- as.numeric(mPW.out1[4])
(Linfbias <- 100* (Linf_estim - Linf_input)/Linf_input )
## [1] 4.882692
(ZKbias <- 100* (ZK_estim - ZK_input)/ZK_input)
## [1] 11.72642