NOp<-FALSE
N1<-10000
N2<-20
a<-1000
A<-c(rep(1,a),rep(2,N1-a))
A<-sample(A)
xlim<-ylim<-c(0,1)
bs<-seq(from=1,to=N2-1,by=1)
t<-seq(from=0,to=1,length=100)*2*pi
ra<-0.2
rb<-ra
max<-0.25
may<-0.5
mbx<-0.75
mby<-0.5
p.out<-rep(0,length(bs))
for(i in 1:length(bs)){
b<-bs[i]
B<-c(rep(1,b),rep(2,N2-b))
B<-sample(B)
plot(ra*cos(t)+max,ra*sin(t)+may,type="l",xlim=xlim,ylim=ylim,main="",xlab="",ylab="")
par(new=TRUE)
plot(rb*cos(t)+mbx,rb*sin(t)+mby,type="l",xlim=xlim,ylim=ylim,main="",xlab="",ylab="")
par(new=TRUE)
rs<-runif(N1)*ra*0.9
ts<-runif(N1)*2*pi
plot(rs*cos(ts)+max,rs*sin(ts)+may,pch=19,cex=0.1,col=A,xlim=xlim,ylim=ylim,main="",xlab="",ylab="")
par(new=TRUE)
rs<-runif(N2)*rb*0.9
ts<-runif(N2)*2*pi
v1<-c(rep(1,length(A)),rep(2,length(B)))
v2<-c(A,B)
print(table(v1,v2))
p.f<-fisher.test(table(v1,v2))$p.value
p.out[i]<-p.f
if(NOp)p.f=""
plot(rs*cos(ts)+mbx,rs*sin(ts)+mby,pch=19,cex=1,col=B,xlim=xlim,ylim=ylim,main=p.f,xlab="",ylab="")
par(ask=TRUE)
}
plot(bs,p.out)