#generalized Polya Urn 

aX<-10             # aX=\alpha(X) 

aN<-function() {  # alpha()/alpha(X), the normalized measure
    rnorm(1,0,1)  # put in different types of distributions here
                }

# Part 1  : convergence of empirical distribution 
X.samp<-NULL
n.samp<-2000
for(n in 1:n.samp) {

new.samp<-rbinom(1,1, aX/(aX+n-1)  )

if(new.samp)  { X.samp<-c(X.samp,aN()) }
if(!new.samp) { X.samp<-c(X.samp,sample(X.samp,1)) }

                   }

par(mfrow=c(3,2))
for(n.part  in c(25,50,100,500,1000,2000)) {
plot(sort(X.samp[1:n.part]),(1:n.part)/n.part,
     type="l",xlab="x",ylab="") 
mtext( paste("n=",n.part),side=3)
                                     }

# Part 2  :  range of Dirichlet Processes
for( aX in c(10,100) ) {
par(mfrow=c(1,1))
plot(c(-2.5,2.5),c(0,1),type="n",xlab="x",ylab="")
for(j in 1:50) {
X.samp<-NULL
n.samp<-2000
for(n in 1:n.samp) {
new.samp<-rbinom(1,1, aX/(aX+n-1)  )
if(new.samp)  { X.samp<-c(X.samp,aN()) }
if(!new.samp) { X.samp<-c(X.samp,sample(X.samp,1)) }
                   }
lines(sort(X.samp[1:n.samp]),(1:n.samp)/n.samp, col="orange",lty=2)
                 }
x<-seq(-2.5,2.5,length=100)
lines( x, pnorm(x), col="red",lwd=2)
mtext( paste( "alpha(X) = ", aX) ,side=3) 
                          }