R Exercise for 10/11 These are three functions that you can copy/paste into your R environment and then run to see what happens. They work by simulating a Brownian motion (by interpolating it by a Gaussian random walk). 1. The first, "brownianup" simulates a Brownian Motion on an axis, for 1000 steps of 0.001 in time, with variance 0.001 per step (so total variance is 1.0) for all 1000 steps. It plots it with the variable horizontally, with time vertically. It is similar to what you ran last time. plotbrownianup <- function(){ x <- rnorm(1000,0,sd=sqrt(0.001)); y <- cumsum(x); yl <- 2.5; t <- seq(0.001,1,0.001); plot(y,t,type="l",xlim=c(-yl,yl),xlab="x",ylab="time"); } 2. The second simulates a phylogeny, with a lineage undergoing Brownian motion for 1000 steps as above, but with the lineage splitting after a fraction (a) of that time, where (a) is a parameter that you choose. The two lineages start from that place and then evolve independently by the Brownian motion. The lineages are plotted, with time vertically, for a single tree. plottwotree <- function(a){ # argument (a) is the fraction of time that the two lineages share ancestry sd1 <- sqrt(0.001); xl <- 3.0; shtime <- trunc(a*1000); indtime <- 1000-shtime; rx1 <- rnorm(shtime,0,sd1); rx2 <- rnorm(indtime,0,sd1); rx3 <- rnorm(indtime,0,sd1); x <- cumsum(rx1); y <- x[shtime]+cumsum(rx2); z <- x[shtime]+cumsum(rx3); t <- seq(0.001,1,0.001); plot(x,t[1:shtime],type="l",xlim=c(-xl,xl),ylim=c(0,1),xlab="x+y or x+z",ylab="time"); lines(y,t[(shtime+1):1000]); lines(z,t[(shtime+1):1000]); } 3. The third function does the same, with the same parameter (a), but now it instead simulates 1000 trees, and plots a scatter plot of the results after the 1000 steps, with 400 replicate trees. It shows the phenotypes in the two lineages, one being x and one being y. How does the correlation between them depend on the value of (a) ? plottwocorr <- function(a){ # argument (a) is fraction of time that the two lineages share ancestry reps <- 400; sd1 <- sqrt(0.001); x <- rep(0,reps); y <- rep(0,reps); xl <- 5.0; shtime <- trunc(a*1000); indtime <- 1000-shtime; for(i in 1:reps) { rx1 <- rnorm(shtime,0,sd1); rx2 <- rnorm(indtime,0,sd1); rx3 <- rnorm(indtime,0,sd1); x[i] <- cumsum(c(rx1,rx2))[1000]; y[i] <- cumsum(c(rx1,rx3))[1000]; } plot(x,y,type="p",xlim=c(-xl,xl),ylim=c(-xl,xl), xlab="final x+y",ylab="final x+z"); }