# Hal Caswell: Matrix Population Models, Chapter 2, Examples ############################################################################### # 2.1 Konstantní projekční matice KPM <- function(T, # délka projekce A=array(c(0,0.3,0,1,0,0.5,5,0,0),c(3,3)), # projekční matice n0=array(c(1,0,0),c(3,1)) # počáteční stav ){ n <- array(NA,c(3,T+1)) n[,1] <- n0 for (i in 1:T){n[,i+1] <- A %*% n[,i]} print(A);print(n) split.screen(c(2,1)) screen(1) plot(c(0,T),c(0,max(n)),type="n",xlab="Time",ylab="Numbers",bty="l") points(seq(0,T),n[1,],type="l",col="red") points(seq(0,T),n[2,],type="l",col="green") points(seq(0,T),n[3,],type="l",col="blue") text(T,n[1,T+1]-0.05*max(n),"1",col="red") text(T,n[2,T+1]-0.05*max(n),"2",col="green") text(T,n[3,T+1]-0.05*max(n),"3",col="blue") N <- array(NA,T+1) for (i in 1:(T+1)){N[i] <- sum(n[,i])} screen(2) plot(c(0,T),c(0,1),type="n",xlab="Time",ylab="Proportions",bty="l") points(seq(0,T),n[1,]/N,type="l",col="red") points(seq(0,T),n[2,]/N,type="l",col="green") points(seq(0,T),n[3,]/N,type="l",col="blue") text(T,n[1,T+1]/N[T+1]-0.05,"1",col="red") text(T,n[2,T+1]/N[T+1]-0.05,"2",col="green") text(T,n[3,T+1]/N[T+1]-0.05,"3",col="blue") close.screen(all=TRUE) } ############################################################################### # 2.2 Vliv počátečních podmínek VPP <- function(T=40, # délka projekce A=array(c(0,0.3,0,1,0,0.5,5,0,0),c(3,3)), # projekční matice NoSim=10, # počet simulací Census=FALSE # zobrazit součet? ){ print(A) n <- array(NA,c(3,T+1)) N <- array(NA,T+1) # Odhad maximální velikosti populace n[,1] <- array(c(0,0,1.75),c(3,1)) N[1] <- 1.75 for (i in 1:T){n[,i+1] <- A %*% n[,i] N[i+1] <- sum(n[,i+1])} mx <- max(n); MX <- max(N) if (Census){split.screen(c(2,2))} else{split.screen(c(2,1))} # Simulace for (S in 1:NoSim){ n[,1] <- array(runif(3,0,1),c(3,1)) N[1] <- sum(n[,1]) for (i in 1:T){n[,i+1] <- A %*% n[,i] N[i+1] <- sum(n[,i+1])} screen(1) plot(c(0,T),c(0,mx),type="n",xlab="Time",ylab="Numbers",bty="l") points(seq(0,T),n[1,],type="l",col="red") points(seq(0,T),n[2,],type="l",col="green") points(seq(0,T),n[3,],type="l",col="blue") screen(2) plot(c(0,T),c(0,1),type="n",xlab="Time",ylab="Proportions",bty="l") points(seq(0,T),n[1,]/N,type="l",col="red") points(seq(0,T),n[2,]/N,type="l",col="green") points(seq(0,T),n[3,]/N,type="l",col="blue") if (Census){ screen(3) plot(c(0,T),c(0,MX),type="n",xlab="Time",ylab="Census",bty="l") points(seq(0,T),N,type="l") screen(4) plot(c(0,T),c(0.1,MX),type="n",xlab="Time",ylab="Census",log="y",bty="l") points(seq(0,T),N,type="l")} } close.screen(all=TRUE) } ############################################################################### # 2.3 Vliv změny koeficientů VZK <- function(T=40, # délka projekce A=array(c(0,0.3,0,1,0,0.5,5,0,0),c(3,3)), # projekční matice n0=array(c(1,0,0),c(3,1)), # počáteční stav P=-0.1 # perturbace ){ print(A);print(n0) Ap <- array(NA,c(3,3,6)) n <- array(NA,c(3,T+1,6)) N <- array(NA,c(T+1,6)) for (k in 1:5){Ap[,,k] <- A} Ap[1,2,2] <- A[1,2]*(1+P) Ap[1,3,3] <- A[1,3]*(1+P) Ap[2,1,4] <- A[2,1]*(1+P) Ap[3,2,5] <- A[3,2]*(1+P) Ap[,,6]<- A*(1+P) for (k in 1:6){ n[,1,k] <- n0 N[1,k] <- sum(n0) for (i in 1:T){n[,i+1,k] <- Ap[,,k] %*% n[,i,k] N[i+1,k] <- sum(n[,i+1,k])} } plot(seq(0,T),N[,1],type="l",xlab="Time",ylab="Census",bty="l", ylim=c(-0.05*max(N),max(N))) points(seq(0,T),N[,6],type="l",lty=2) text(0,-0.05*max(N),"F_2",col="orange") points(seq(0,T),N[,2],type="l",col="orange") text(0.1*T,-0.05*max(N),"F_3",col="brown") points(seq(0,T),N[,3],type="l",col="brown") text(0.2*T,-0.05*max(N),"P_1",col="violet") points(seq(0,T),N[,4],type="l",col="violet") text(0.3*T,-0.05*max(N),"P_2",col="cyan") points(seq(0,T),N[,5],type="l",col="cyan") } ############################################################################### # 2.4 Vliv nákodné změny koeficientů VNZ <- function(T=40, # délka projekce A=array(c(0,0.3,0,1,0,0.5,5,0,0),c(3,3)), # projekční matice n0=array(c(1,0,0),c(3,1)), # počáteční stav NoSim=1, # počet simulací Ampl=2 # relativní odchylka ){ print(A);print(n0) N <- array(NA,c(T+1,NoSim)) n <- array(NA,c(3,T+1)) n[,1] <- n0 for (k in 1:NoSim){ h <- runif(T,-1,1) N[1,k] <- sum(n0) for (i in 1:T){Ap <- A; Ap[1,] <- A[1,]*(Ampl^sign(h[i])) n[,i+1] <- Ap %*% n[,i] N[i+1,k] <- sum(n[,i+1])} } Nav <- array(sum(n0),T+1) for (i in 1:T){n[,i+1] <- A %*% n[,i] Nav[i+1] <- sum(n[,i+1])} if (NoSim > 1){Mean <- array(NA,T+1) Var <- array(NA,T+1) for (i in 1:(T+1)){Mean[i] <- mean(N[i,]) Var[i] <- var(N[i,])} split.screen(c(2,1)) screen(1)} plot(seq(0,T),N[,1],type="l",xlab="Time",ylab="Census",log="y",bty="l", ylim=c(min(N),max(N))) if (NoSim > 1){for (k in 1:NoSim){points(seq(0,T),N[,k],type="l")} points(seq(0,T),Mean,type="l",col="orangered",lwd=2)} points(seq(0,T),Nav,type="l",col="cyan",lwd=2) if (NoSim > 1){screen(2) yl <- c(min(Mean,Var[which(Var>0)],Nav),max(Mean,Var,Nav)) plot(seq(0,T),Mean,type="l",col="orangered",xlab="Time",ylab="Mean, Variance", bty="l",ylim=yl,log="y") points(seq(1,T),Var[1:T],type="l",col="magenta") points(seq(0,T),Nav,type="l",col="cyan") text(0,yl[2],adj=0,"Variance",col="magenta") text(0,yl[2]^0.85,adj=0,"Mean",col="orangered") text(0,yl[2]^0.65,adj=0,"Unperturbed matrix",col="cyan") close.screen(all=TRUE)} } ############################################################################### # 2.5 Vliv velikosti populace (density dependence) VVP <- function(T=40, # délka projekce A=array(c(0,0.3,0,1,0,0.5,5,0,0),c(3,3)), # projekční matice n0=array(c(1,0,0),c(3,1)), # počáteční stav b=0.005,R=1 # parametry funkce g ){ print(A); print(n0) n <- array(NA,c(3,T+1)) N <- array(NA,T+1) n[,1] <- n0 N[1] <- sum(n[,1]) for (i in 1:T){g <- R*exp(-b*N[i]) Ap <- A; Ap[1,2] <- A[1,2]*g; Ap[1,3] <- A[1,3]*g n[,i+1] <- Ap %*% n[,i] N[i+1] <- sum(n[,i+1])} split.screen(c(2,1)) screen(1) plot(seq(0,T),N,type="l",xlab="Time",ylab="Numbers",bty="l",ylim=c(0,max(N))) points(seq(0,T),n[1,],type="l",col="red") points(seq(0,T),n[2,],type="l",col="green") points(seq(0,T),n[3,],type="l",col="blue") screen(2) plot(c(0,T),c(0,1),type="n",xlab="Time",ylab="Proportions",bty="l") points(seq(0,T),n[1,]/N,type="l",col="red") points(seq(0,T),n[2,]/N,type="l",col="green") points(seq(0,T),n[3,]/N,type="l",col="blue") close.screen(all=TRUE) }