z3mao <zhangyij <at> gmail.com> writes:
>
>
> Hi, this is my first time using R. I want to simulate the following process:
> "in a population of size N, there are i individuals bearing genotype A, the
> number of those bearing A is j in the next generation, which following a
> binominal distribution (choose j from 2*N, the p is i/2*N), to plot the
> probability of the next generations, my script is as follows. It cannot run
> successfully, declaring that the "ylim should be limited. " I wonder where
> the bug is. Thanks very much!
>
> generation<-function(i,N)
> {
> m<-1;gen<-numeric();
> for(m in 1:50)
> {
> testp<-runif(1,0,1);
> j<-0; sump<-0;
> while(sump < testp)
> { sump<-sump+dbinom(j,2*N,i/(2*N));
> j<-j+1;
> }
> i<-j; gen[m]<-j/(2*N); m<-m+1;
> }
> plot(m, gen[m]);
> }
## basic sim
N <- 50
ngen <- 500
nA <- numeric(ngen)
nA[1] <- 25
for (i in 2:ngen) {
fixed <- nA[i-1]==0 || nA[i-1]==N
if (fixed) break
nA[i] <- rbinom(1,prob=nA[i-1]/N,size=N)
}
## wrap it in a function
binsim <- function(N,ngen,A0) {
nA <- numeric(ngen)
nA[1] <- A0
for (i in 2:ngen) {
fixed <- nA[i-1]==0 || nA[i-1]==N
if (fixed) break
nA[i] <- rbinom(1,prob=nA[i-1]/N,size=N)
}
nA
}
## replicate it
m <- replicate(200,binsim(50,100,25))
matplot(m,col="grey",type="l",lty=1,pch=1)
means <- rowMeans(m)
quants <- t(apply(m,1,quantile,c(0.025,0.975)))
lines(means,lwd=2)
lines(rowMeans(m)-,lwd=2)
matlines(quants,lty=2,lwd=2,col=1)
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
