#
# Nonparametric inference for a Dirichlet Process example.
# 
# Consider data generated from a beta distribution.
#
n <- 20
x <- rbeta(n,2,1)
xgrid <- c(0:10000)/10000
#
#
#
F0 <- xgrid
a <- 5
Fhat <- rep(0,length(xgrid))
Fpost <- rep(0,length(xgrid))
Ftrue <- pbeta(xgrid,2,1)
for (i in 1:10001){
  Fhat[i] <- sum(x<=xgrid[i])/n
  Fpost[i] <- (a*F0[i]+n*Fhat[i])/(a+n)
}
plot(xgrid,F0,type='l',lwd=5,col='green',xlab='x',ylab='F')
lines(xgrid,Fhat,type='l',lwd=5,col='blue')
lines(xgrid,Ftrue,type='l',lwd=5,col='black')
#
# Now a discrete example with a Dirichlet process prior.  
# In this case we can estimate the probability function exactly.
# 
# Consider data generated from a Poisson distribution.
#
n <- 20
r <- 10
p <- 0.7
x <- rnbinom(n,r,p)
plot(table(x)/length(x),xlim=c(0,(max(x)+2)),ylab='p',lwd=5)
xgrid <- c(0:(max(x)+2))
#
# Use a Poisson(lambda) prior mean.
#
a <- 5
lambda <- 10
p0 <- dpois(xgrid,lambda)
F0 <- ppois(xgrid,lambda)
lines(xgrid-0.2,p0,type='h',col='blue',lwd=5)
Fhat <- rep(0,length(xgrid))
Fpost <- rep(0,length(xgrid))
for (i in 1:(max(x)+3)){
  Fhat[i] <- sum(x<=(i-1))/n
  Fpost[i] <- (a*F0[i]+n*Fhat[i])/(a+n)
}
ppost <- rep(0,length(xgrid))
ppost[1] <- Fpost[1]
for (i in 2:(max(x)+3)){
  ppost[i]<- Fpost[i]-Fpost[i-1]
}
lines(xgrid+0.2,ppost,type='h',col='red',lwd=5)