################################################### ### chunk number 1: ################################################### n=100 x = runif(n,min=0,max=10) a=1 b=0.5 s=3 y_det=a*x*exp(-b*x) y = rgamma(n,shape=s,scale=y_det/s) ################################################### ### chunk number 2: ################################################### plot(x,y) curve(a*x*exp(-b*x),add=TRUE) ################################################### ### chunk number 3: ################################################### set.seed(1001) nparents = 50 noffspr = 10 L = 30 parent_x = runif(nparents,min=0,max=L) parent_y = runif(nparents,min=0,max=L) angle = runif(nparents*noffspr,min=0,max=2*pi) dist = rexp(nparents*noffspr,0.5) offspr_x = rep(parent_x,each=noffspr)+cos(angle)*dist offspr_y = rep(parent_y,each=noffspr)+sin(angle)*dist dist = sqrt((outer(offspr_x,offspr_x,"-"))^2+(outer(offspr_y,offspr_y,"-"))^2) nbrcrowd = apply(dist<2,1,sum)-1 ################################################### ### chunk number 4: ################################################### m = mean(nbrcrowd) s2 = var(nbrcrowd) ################################################### ### chunk number 5: ################################################### k.est = m/(s2/m-1) ################################################### ### chunk number 6: ################################################### b1 = barplot(table(factor(nbrcrowd,levels=0:max(nbrcrowd)))/length(nbrcrowd), xlab="Number of neighbors",ylab="Proportion") points(b1,dnbinom(0:max(nbrcrowd),mu=m,size=k.est),pch=16) ################################################### ### chunk number 7: ################################################### ci = nbrcrowd*3 M=2.3 alpha=0.49 mass_det=M/(1+ci) mass = rgamma(length(mass_det),scale=mass_det,shape=alpha) b = 271.6; k= 0.569 seed_det = b*mass seed = rnbinom(length(seed_det),mu=seed_det,size=k) ################################################### ### chunk number 8: ################################################### logxvec = seq(-7,1,length=100) xvec = 10^logxvec med = qnbinom(0.5,mu=b*xvec,size=k) ################################################### ### chunk number 9: ################################################### plot(mass,1+seed,log="xy",xlab="Mass",ylab="1+Seed set") curve(b*x+1,add=TRUE) lower = qnbinom(0.025,mu=b*xvec,size=k) upper = qnbinom(0.975,mu=b*xvec,size=k) lines(xvec,lower+1,lty=2,type="s") lines(xvec,upper+1,lty=2,type="s") lines(xvec,med+1,lwd=2,type="s") ################################################### ### chunk number 10: ################################################### rzinbinom = function(n,mu,size,zprob) { ifelse(runif(n)immig) & (powsimresults[,"slope.lo"]null.value) slope.pow = tapply(reject.null,faclist,sum)/nsim ################################################### ### chunk number 22: ################################################### par(mfrow=c(2,2)) matplot(nvec,slope.mean,type="b") matplot(nvec,slope.sd,type="b") matplot(nvec,slope.cov,type="b") matplot(nvec,slope.pow,type="b") ################################################### ### chunk number 23: ################################################### nvec=c(3,5,7,10,15,20) nsim=500 powsimresults = matrix(nrow=length(nvec)*nsim,ncol=5) colnames(powsimresults) = c("n","sim","quad","quad.lo","quad.hi") ctr=1 for (i in 1:length(nvec)) { nt = nvec[i] tvec=1:nt for (sim in 1:nsim) { current.sim = immigsim(nt=nt,N0=N0,surv=p,immig=immig) lm1=lm(current.sim~tvec+I(tvec^2)) quad=coef(lm1)[3] ci.quad=confint(lm1)[3,] powsimresults[ctr,] = c(nt,sim,quad,ci.quad) ctr = ctr+1 } } ################################################### ### chunk number 24: ################################################### quad.mean = tapply(powsimresults[,"quad"],nfac,mean) quad.sd = tapply(powsimresults[,"quad"],nfac,sd) nsim=500 null.value=0 reject.null = (powsimresults[,"quad.hi"]null.value) quad.pow = tapply(reject.null,nfac,sum)/nsim ################################################### ### chunk number 25: ################################################### op = par(mfrow=c(2,2)) plot(nvec,quad.mean,type="b") plot(nvec,quad.sd,type="b") plot(nvec,quad.pow,type="b") par(op)