<\body> > Help: you can write a series of R commands into a file, and then execute them all in a row by doing <\input| >> source("filename") > \; 1a. By simulating 10000 times flipping a coin 10 times, and storing the result in a 10> matrix, make a table of the probability to have heads show up 0 times, 1 time, 2 times, and so on. <\output> \; <\input| >> flip=matrix(sample(1:2,10000*10,rep=T),10000,10) <\input| >> count=apply(flip==1,1,sum) <\input| >> hist(count,breaks=seq(-0.5,10.5,by=1));v() <\output> |ps>||||||> <\input| >> sapply(0:10,function(i) sum(count==i)) <\output> \ [1] \ \ 13 \ 121 \ 463 1116 1998 2494 2031 1210 \ 456 \ \ 88 \ \ 10 <\input| >> table(count)/length(count) <\output> count \ \ \ \ \ 0 \ \ \ \ \ 1 \ \ \ \ \ 2 \ \ \ \ \ 3 \ \ \ \ \ 4 \ \ \ \ \ 5 \ \ \ \ \ 6 \ \ \ \ \ 7 \ \ \ \ \ 8 \ \ \ \ \ 9 \ \ \ \ 10\ 0.0013 0.0121 0.0463 0.1116 0.1998 0.2494 0.2031 0.1210 0.0456 0.0088 0.0010\ <\input| >> \; > \; 1b. If someone flipped a coin, and got 1 head, and 9 tails, we would like to test the hypothesis that the coin is unbiased. We will do this by calculating the probability that an unbiased coin flipped 10 times gives 0 or 1 heads, or 0 or 1 tails. Calculate this probability in the simulation above. <\input| >> sum(count \= 1 \| count \= 10-1)/length(count) <\output> [1] 0.0232 <\input| >> \; > \; 1c. How many heads/tails out of 10 would cause one to reject the hypothesis that the coin is unbiased at the 5% level? This means the following: what are the 5% most extreme coin flips: how much frequency is in 0-10, 10-0, how much in 0-10, 10-0, 9-1, 1-9, and how much if you also add 2-8, 8-2, and so on, and when does this get bigger than 5%. If you see a flip in that extreme group that you found, you can say that such an extreme bias has only a 5% or lower chance to happen.\ <\input| >> Tab=table(count)/length(count) <\input| >> Tab <\output> count \ \ \ \ \ 0 \ \ \ \ \ 1 \ \ \ \ \ 2 \ \ \ \ \ 3 \ \ \ \ \ 4 \ \ \ \ \ 5 \ \ \ \ \ 6 \ \ \ \ \ 7 \ \ \ \ \ 8 \ \ \ \ \ 9 \ \ \ \ 10\ 0.0013 0.0121 0.0463 0.1116 0.1998 0.2494 0.2031 0.1210 0.0456 0.0088 0.0010\ <\input| >> sapply(0:5,function(i) sum( count \= i \| count \= 10-i ) )/length(count) <\output> [1] 0.0023 0.0232 0.1151 0.3477 0.7506 1.0000 <\input| >> \; > \; 1d. With the command \ you can simulate a coin flip of a biased coin. Out of 1000 coin flips of length 10 that are biased as above, how many come up significant at the 5% level? I.e. how often do you get as biased as what you calculated in 1c? Do the same for an unbiased coin - how many come up significant for an unbiased coin ? How about a more biased coin, with probabilities ? \; <\input| >> bflip=matrix( sample(0:1, 10*1000, rep=T, prob=c(0.3,0.7) ), 1000, 10) <\input| >> bcount=apply(bflip,1,sum) <\input| >> \; > Since we reject at the 5% level 0:10, 1:9, 9:1, and 10:0, we count how often these occur: <\input| >> sum(bcount \= 1 \| bcount \= 10-1)/length(bcount) <\output> [1] 0.147 <\input| >> \; > So in 14% of the cases we recognize the biased coin as biased. Notice that an unbiased coin we would mark as biased in 5% of the cases. <\input| >> bflip=matrix( sample(0:1, 10*1000, rep=T, prob=c(0.1,0.9) ), 1000, 10) <\input| >> bcount=apply(bflip,1,sum) <\input| >> sum(bcount \= 1 \| bcount \= 10-1)/length(bcount) <\output> [1] 0.75 <\input| >> \; > In this case, we recognize the coin as biased in 74% of the cases. We can try to do a graph of this: <\input| >> recog=function(p) {\ \ bflip=matrix( sample(0:1, 10*1000, rep=T, prob=c(p,1-p) ), 1000, 10) \ bcount=apply(bflip,1,sum) \ sum(bcount \= 1 \| bcount \= 10-1)/length(bcount) } <\input| >> recog(0.1) <\output> [1] 0.739 <\input| >> x=seq(0,1,length=50) <\input| >> y=sapply(x,recog) <\input| >> plot(x,y,main="Chance of recognition of a biased coin");v() <\output> |ps>||||||> <\input| >> \; > \; \; 2. If the sex ratio at birth is 50%, but people adopt the following strategy: have children until you have one son, and then stop, unless you have 12 daughters (in which case you stop, too). What would be the average sex ratio in the population? (Calculate by simulation!) <\with|prog-language|r|prog-session|default> <\session> <\input| >> kids=matrix(sample(0:1,10000*12,rep=T),10000,12) <\input| >> children=apply(kids,1,function(x) {\ \ \ \ if( sum(x)==0 ) {\ \ \ \ \ \ c(daughters=12,sons=0)\ \ \ \ } else {\ \ \ \ \ \ c(daughters=min( which( x==1 ) )-1, sons=1) \ \ \ } } ) <\input|>> \; <\input| \ >> dim(children) <\output> [1] \ \ \ \ 2 10000 <\input| >> children[,1:10] <\output> \ \ \ \ \ \ \ \ \ \ [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] daughters \ \ \ 2 \ \ \ 1 \ \ \ 3 \ \ \ 1 \ \ \ 0 \ \ \ 0 \ \ \ 1 \ \ \ 0 \ \ \ 1 \ \ \ \ 1 sons \ \ \ \ \ \ \ \ 1 \ \ \ 1 \ \ \ 1 \ \ \ 1 \ \ \ 1 \ \ \ 1 \ \ \ 1 \ \ \ 1 \ \ \ 1 \ \ \ \ 1 <\input|>> \; <\input| >> sons=children["sons",];daughters=children["daughters",] <\input| >> table(sons,daughters) <\output> \ \ \ \ daughters sons \ \ \ 0 \ \ \ 1 \ \ \ 2 \ \ \ 3 \ \ \ 4 \ \ \ 5 \ \ \ 6 \ \ \ 7 \ \ \ 8 \ \ \ 9 \ \ 10 \ \ 11 \ \ 12 \ \ \ 0 \ \ \ 0 \ \ \ 0 \ \ \ 0 \ \ \ 0 \ \ \ 0 \ \ \ 0 \ \ \ 0 \ \ \ 0 \ \ \ 0 \ \ \ 0 \ \ \ 0 \ \ \ 0 \ \ \ 5 \ \ \ 1 5044 2491 1225 \ 648 \ 278 \ 156 \ \ 83 \ \ 40 \ \ 18 \ \ \ 7 \ \ \ 4 \ \ \ 1 \ \ \ 0 <\input| >> sum(sons)/sum(daughters) <\output> [1] 1.012357 <\input| >> \; \; \; We see that there seem to be equal numbers of sons and daughters. Let us test this: <\input| >> prop.test(sum(sons),n=sum(sons)+sum(daughters),p=0.5) <\output> \; \ \ \ \ \ \ \ \ 1-sample proportions test with continuity correction \; data: \ sum(sons) out of sum(sons) + sum(daughters), null probability 0.5\ X-squared = 0.7369, df = 1, p-value = 0.3907 alternative hypothesis: true p is not equal to 0.5\ 95 percent confidence interval: \ 0.4960928 0.5100465\ sample estimates: \ \ \ \ \ \ \ \ p\ 0.5030703\ <\input| >> \; > We don't see a significant difference from equal probabilities. \; \; 3a. Similar to what we did in the first exercise, start with a population of 100 individuals. We will assume that they all have different alleles, so that their alleles are 1 for the first individual, 2 for the second, 3 for the third, and so on. Now simulate neutral evolution. Create a new population by sampling from this one 100 individual, then another one by sampling from the second, and so on. Do this for 300 generations. How many of the original 100 alleles are left? <\input| >> pop=1:100 <\input| >> new.gen=function(pop){sample(pop,rep=T)} <\input| >> for(i in 1:300){ pop=new.gen(pop) } <\input| >> pop <\output> \ \ [1] 72 72 72 72 72 72 72 72 72 72 72 72 72 72 72 72 72 72 72 72 72 72 72 72 72 \ [26] 72 72 72 72 72 72 72 72 72 72 72 72 72 72 72 72 72 72 72 72 72 72 72 72 72 \ [51] 72 72 72 72 72 72 72 72 72 72 72 72 72 72 72 72 72 72 72 72 72 72 72 72 72 \ [76] 72 72 72 72 72 72 72 72 72 72 72 72 72 72 72 72 72 72 72 72 72 72 72 72 72 <\input| >> table(pop) <\output> pop \ 72\ 100\ <\input| >> length(unique(pop)) <\output> [1] 1 <\input| >> \; > Only one genotype survived. unique(x) shows how many different numbers (or strings) are in x. Here is a cool way to display this: <\input| >> pop=1:100 <\input| >> gens=matrix(0,300,100);gens[1,]=pop <\input| >> for(i in 1:299){gens[i+1,]=new.gen(gens[i,])} <\input| >> \; > Now we have the whole timeline. Now lets sort it: <\input| >> gens=t(apply(gens,1,sort)) <\input| >> image(gens,col=rainbow(100));v() <\output> |ps>||||||> <\input| >> \; > 3b. Now represent the population as a matrix, so that the rows are individuals, and the columns are loci \ (i.e. individuals are haploids). create a new population by chosing a mother and a father for each individual, and then chosing the alleles randomly from the father and the mother for each locus. Start with a population of 20 individuals with a genome length of 100. Alleles from how many individuals are present after 300 generations? <\input| >> pop=matrix(1:20,20,100) <\input| >> pop[1:3,1:12] <\output> \ \ \ \ \ [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [1,] \ \ \ 1 \ \ \ 1 \ \ \ 1 \ \ \ 1 \ \ \ 1 \ \ \ 1 \ \ \ 1 \ \ \ 1 \ \ \ 1 \ \ \ \ 1 \ \ \ \ 1 \ \ \ \ 1 [2,] \ \ \ 2 \ \ \ 2 \ \ \ 2 \ \ \ 2 \ \ \ 2 \ \ \ 2 \ \ \ 2 \ \ \ 2 \ \ \ 2 \ \ \ \ 2 \ \ \ \ 2 \ \ \ \ 2 [3,] \ \ \ 3 \ \ \ 3 \ \ \ 3 \ \ \ 3 \ \ \ 3 \ \ \ 3 \ \ \ 3 \ \ \ 3 \ \ \ 3 \ \ \ \ 3 \ \ \ \ 3 \ \ \ \ 3 <\input| >> \; > We created the population. Each row is an individual, each column is a locus. At each locus each individual is distinguishable. We want to see from how many individuals do we still see alleles at the end. <\input| >> mom.i=sample(1:20,rep=T) <\input| >> dad.i=sample(1:20,rep=T) <\input| >> mom=t( sapply( mom.i, function(i) pop[i,] ) ) <\input| >> dad=t( sapply( dad.i, function(i) pop[i,] ) ) <\input| >> \; > Now we have a matrix with the genotypes of the mothers, and a matrix with the genotypes of the fathers. To recombine we will randomly chose either from one, or from the other: <\input| >> rec=matrix(sample(c(T,F), 20*100,rep=T),20,100) <\input| >> next.gen=mom <\input| >> next.gen[rec]=dad[rec] <\input| >> next.gen[1:3,1:12] <\output> \ \ \ \ \ [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [1,] \ \ 12 \ \ \ 3 \ \ 12 \ \ 12 \ \ \ 3 \ \ 12 \ \ 12 \ \ \ 3 \ \ 12 \ \ \ \ 3 \ \ \ 12 \ \ \ \ 3 [2,] \ \ \ 4 \ \ \ 4 \ \ \ 4 \ \ 14 \ \ 14 \ \ \ 4 \ \ 14 \ \ \ 4 \ \ \ 4 \ \ \ 14 \ \ \ 14 \ \ \ \ 4 [3,] \ \ 14 \ \ \ 4 \ \ \ 4 \ \ 14 \ \ 14 \ \ \ 4 \ \ 14 \ \ 14 \ \ 14 \ \ \ \ 4 \ \ \ 14 \ \ \ \ 4 <\input| >> image(t(next.gen),col=rainbow(100));v() <\output> |ps>||||||> <\input| >> \; > That's it! You can see how every line has input from two colors. Now we need to do this 300 times. Let us write it as a function. <\input| >> new.gen2=function(pop){ \ mom.i=sample(1:20,rep=T); dad.i=sample(1:20,rep=T) \ mom=t( sapply( mom.i, function(i) pop[i,] ) ) \ dad=t( sapply( dad.i, function(i) pop[i,] ) ) \ rec=matrix(sample(c(T,F), 20*100,rep=T),20,100) \ next.gen=mom \ next.gen[rec]=dad[rec] \ next.gen } <\input|>> \; <\input| >> image( t(new.gen2(pop)), col=rainbow(100) );v() <\output>