ND<-500              #number of draws to take from each Polya urn 
NU<-100 
par(mfrow=c(2,1))
plot(c(1,ND),c(0,1),type="n",xlab="draw",ylab="p(Red)")
abline(h= (aRG[1,1]-1)/( aRG[1,1] + aRG[1,2] -2 ),col="blue")
PS<-NULL
for(nsim in 1:NU){   #run NU independent urn schemes, each with ND draws
aRG<-matrix(c(6,3),nrow=1,ncol=2)
for(n in 1:ND) {
x<-rbinom(1,1, aRG[n,1]/(sum(aRG[n,])))            # sample the ball 
aRG<-rbind(aRG, c( aRG[n,1]+x , aRG[n,2]+(1-x) ) ) # add a ball of approp 
                                                   # color
                 }
lines( aRG[,1]/( aRG[,1]+aRG[,2]))
PS<-c(PS,aRG[n,1]/( aRG[n,1]+aRG[n,2] ) )          # keep sample proportion 
                                                   # of red balls at the 
                                                   # end of each urn scheme
                   }

hist(PS,prob=T)                                    # plot out the dist 
x<-seq(0,1,length=100)                             # of the "limiting" 
lines( x, dbeta(x,aRG[1,1],aRG[1,2]) ,col="blue" ) # proportion for each urn