-
Notifications
You must be signed in to change notification settings - Fork 5
/
matrixOperators.r
53 lines (39 loc) · 1.53 KB
/
matrixOperators.r
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
# Matrix operators for population models
## maximum lambda function
max.lambda <- function(x) Re((eigen(x)$values)[1]) ## where 'x' is a Leslie matrix
## Maximum r function
max.r <- function(x) log(Re((eigen(x)$values)[1])) ## where 'x' is a Leslie matrix
## Stable stage distribution
stable.stage.dist <- function(x) ((x %*% (Re((eigen(x)$vectors)[,1])))/(sum((x %*% (Re((eigen(x)$vectors)[,1]))))))[,1]
## Generation length function
R.val <- function(X,age.max) ## reproductive value (R0) where X = Leslie matrix; age.max = maximum age of females
{
## define the transition matrix
T <- X[1:age.max,1:age.max]
T[1,1:(age.max)] <- 0
## define the fertility matrix
F <- X[1:age.max,1:age.max]
diag(F[2:age.max,1:(age.max-1)]) <- 0
## define the identity matrix
I <- matrix(data<-0,nrow<-age.max,ncol<-age.max)
diag(I) <- 1
## define the fundamental matrix
library(MASS)
N.fund <- ginv(I-T)
## define the reproductive matrix
R <- F %*% N.fund
## define R0 (number of female offspring produced per female during lifetime)
R0 <- Re((eigen(R)$values)[1])
## output
print("number of female offspring produced per female during its lifetime")
print("_________________________________________________________________")
print(R0)
}
## Mean generation time function
G.val <- function (X,age.max) ## where X is a Leslie Matrix
{
G <- (log(R.val(X,age.max)))/(log(Re((eigen(X)$values)[1])))
print("mean generation time")
print("____________________")
print(G)
}