################################################### ### chunk number 1: ################################################### setwd("~/Documents/Projects/Book/draft/Chap02") rm(list=ls()) badfonts = FALSE .detach.stuff = function() { s1 = grep("package|Autoloads",search()) nattach = length(search()) xattach = search()[c(-1,-s1)] for (i in xattach) eval(substitute(detach(i),list(i=i))) } .detach.stuff() ## ps.options(family="NimbusSan") ps.options(colormodel="cmyk") options(width=75) ## palette(gray((0:8)/8)) library(lattice) trellis.par.set(canonical.theme(color=FALSE)) library(primer) ################################################### ### chunk number 2: ################################################### M <- matrix( 1:4 , nr=2, byrow=T) M N <- matrix( c(10, 20, 30, 40), nr=2) N ################################################### ### chunk number 3: ################################################### 1 * 10 + 2 * 20 ################################################### ### chunk number 4: ################################################### M %*% N ################################################### ### chunk number 5: ################################################### A <- matrix(c(0, 0.5, 20, 0.3, 0, 0, 0, 0.5, 0.9), nr=3, byrow=TRUE) N0 <- matrix(c(100, 250, 50), ncol=1) ################################################### ### chunk number 6: ################################################### N1 <- A%*%N0 N1 ################################################### ### chunk number 7: ################################################### years <- 6 N.projections <- matrix( 0, nrow=nrow(A), ncol=years+1) N.projections[,1] <- N0 ################################################### ### chunk number 8: ################################################### for(i in 1:years) N.projections[,i+1] <- A%*%N.projections[,i] ################################################### ### chunk number 9: Mallgrow ################################################### matplot(0:years, t(N.projections), type="l", lty=1:3, col=1, ylab="Stage Abundance", xlab="Year") legend('topleft', legend=c("Seeds", "Small Adult","Large Adult"), lty=1:3, col=1, bty='n') ################################################### ### chunk number 10: ################################################### N.totals <- apply(N.projections, 2, sum) ################################################### ### chunk number 11: ################################################### Rs <- N.totals[-1] / N.totals[-(years+1)] ################################################### ### chunk number 12: MallR ################################################### plot(0:(years-1), Rs, type="b", xlab="Year", ylab="R") ################################################### ### chunk number 13: ################################################### eigs.A <- eigen(A) eigs.A ################################################### ### chunk number 14: ################################################### dom.pos <- which.max( eigs.A[["values"]] ) ################################################### ### chunk number 15: ################################################### L1 <- Re(eigs.A[["values"]][dom.pos]) L1 ################################################### ### chunk number 16: ################################################### t <- 20 Nt <- N0/sum(N0) ################################################### ### chunk number 17: ################################################### R.t <- numeric(t) for(i in 1:t) R.t[i] <- {Nt1 <- A%*%Nt R <- sum(Nt1)/sum(Nt) Nt <- Nt1/sum(Nt1) R } ################################################### ### chunk number 18: Converge ################################################### par(mar=c(5,4,3,2)) plot(1:t, R.t, type='b', main=quote("Convergence Toward "*lambda)) points(t, L1, pch=19, cex=1.5) ################################################### ### chunk number 19: ################################################### w <- Re(eigs.A[["vectors"]][,dom.pos]) ssd <- w/sum(w) round(ssd, 3) ################################################### ### chunk number 20: ################################################### M <- eigen(t(A)) v <- Re(M$vectors[,which.max(Re(M$values))]) RV <- v/v[1] RV ################################################### ### chunk number 21: ################################################### vw.s <- v %*% t(w) ################################################### ### chunk number 22: ################################################### (S <- vw.s/as.numeric(v%*%w)) ################################################### ### chunk number 23: ################################################### elas <- (A/L1) * S round(elas, 3) ################################################### ### chunk number 24: ################################################### data(stagedat) data(fruitdat) data(seeddat) ################################################### ### chunk number 25: ################################################### str(stagedat) ################################################### ### chunk number 26: ################################################### table(stagedat[["Y2004"]]) ################################################### ### chunk number 27: ################################################### str(fruitdat); ################################################### ### chunk number 28: ################################################### table(fruitdat[["Stage"]], fruitdat[["Y2004"]]) ################################################### ### chunk number 29: ################################################### table(seeddat) ################################################### ### chunk number 30: ################################################### mat1 <- matrix(0, nrow=5, ncol=5) ################################################### ### chunk number 31: ################################################### ferts <- tapply(fruitdat$Y2004,fruitdat$Stage, mean ) / 2 ferts ################################################### ### chunk number 32: ################################################### mat1[1,4] <- ferts[1] mat1[1,5] <- ferts[2] ################################################### ### chunk number 33: ################################################### seed.freqs <- table(seeddat[,1]) seedfates <- seed.freqs / length(seeddat[,1]); seedfates ################################################### ### chunk number 34: ################################################### mat1[1,1] <- seedfates[2] mat1[2,1] <- seedfates[3] ################################################### ### chunk number 35: ################################################### for(i in 2:5) { for(j in 2:5) mat1[i, j] <- { x <-subset(stagedat, stagedat$Y2003 == j) jT <- nrow(x) iT <- sum(x$Y2004 == i) iT/jT } } ################################################### ### chunk number 36: ################################################### round(mat1,2) ################################################### ### chunk number 37: ################################################### #source("DemoFunc.R") ################################################### ### chunk number 38: eval=FALSE ################################################### ## ProjMat(stagedat, fruitdat, seeddat) ################################################### ### chunk number 39: ################################################### str( DemoInfo(mat1) ) ################################################### ### chunk number 40: ################################################### nL <- nrow(stagedat) nF <- nrow(fruitdat) nS <- nrow(seeddat) ################################################### ### chunk number 41: ################################################### n <- 5 ################################################### ### chunk number 42: ################################################### n <- 5 out <- lapply(1:n, function(i) { stageR <- stagedat[sample(1:nL, nL, replace=TRUE),] fruitR <- fruitdat[sample(1:nF, nF, replace=TRUE),] seedR <- as.data.frame(seeddat[sample(1:nS,nS, replace=TRUE),]) matR <- ProjMat(stagedat=stageR, fruitdat=fruitR, seeddat=seedR) DemoInfo( matR ) } ) ################################################### ### chunk number 43: ################################################### sapply(out, function(x) x$lambda) ################################################### ### chunk number 44: ################################################### args(DemoBoot) ################################################### ### chunk number 45: ################################################### estims <- DemoInfo( ProjMat(stagedat,fruitdat,seeddat)) estims$lambda ################################################### ### chunk number 46: ################################################### round( estims$Elasticities, 4) ################################################### ### chunk number 47: ################################################### system.time( out.boot <- DemoBoot(stagedat,fruitdat,seeddat,n=1000) ) lambdas <- sapply(out.boot, function(out) out$lambda) ################################################### ### chunk number 48: ################################################### myQ(mr=c(5,4,3,2)) ################################################### ### chunk number 49: ################################################### hist(lambdas, prob=T) lines(density(lambdas)) ################################################### ### chunk number 50: ################################################### dev.print(pdf, "lambdaCI.pdf") ################################################### ### chunk number 51: ################################################### alpha <- 0.05 quantile(lambdas, c(alpha/2, .5, 1-alpha/2)) ################################################### ### chunk number 52: ################################################### bias <- mean(lambdas)-estims$lambda bias ################################################### ### chunk number 53: ################################################### quantile(lambdas-bias, c(alpha/2, .5, 1-alpha/2))