####################################
# Bivariate example:
# Censored data is artificially
# created
####################################
data(airquality)
attach(airquality)
ozone <- Ozone**(1/3)
radiation <- Solar.R
y <- na.omit(cbind(radiation,ozone))
# create censored-data
xxlim <- seq(0,300,50)
yylim <- seq(1.5,5.5,1)
left <- matrix(0,nrow=nrow(y),ncol=2)
right <- matrix(0,nrow=nrow(y),ncol=2)
for(i in 1:nrow(y))
{
left[i,1] <- NA
right[i,1] <- NA
if(y[i,1] < xxlim[1]) right[i,1] <- xxlim[1]
for(j in 1:length(xxlim))
{
if(y[i,1] >= xxlim[j]) left[i,1] <- xxlim[j]
if(y[i,1] >= xxlim[j]) right[i,1] <- xxlim[j+1]
}
left[i,2] <- NA
right[i,2] <- NA
if(y[i,2] < yylim[1]) right[i,2] <- yylim[1]
for(j in 1:length(yylim))
{
if(y[i,2] >= yylim[j]) left[i,2] <- yylim[j]
if(y[i,2] >= yylim[j]) right[i,2] <- yylim[j+1]
}
}
# Prior information
s2 <- matrix(c(10000,0,0,1),ncol=2)
m2 <- c(180,3)
psiinv2 <- diag(c(1/10000,1),2)
prior <- list(alpha=1,nu1=4,nu2=4,s2=s2,
m2=m2,psiinv2=psiinv2,tau1=0.01,tau2=0.01)
# Initial state
state <- NULL
# MCMC parameters
nburn <- 5000
nsave <- 5000
nskip <- 3
ndisplay <- 1000
mcmc <- list(nburn=nburn,
nsave=nsave,
nskip=nskip,
ndisplay=ndisplay)
# Fitting the model
fit1 <- DPMdencens(left=left,right=right,ngrid=100,
prior=prior,mcmc=mcmc,
state=state,status=TRUE)
fit1
summary(fit1)
# Plot the estimated density
plot(fit1)
# Extracting the univariate density estimates
cbind(fit1$grid[,1],fit1$funi[[1]])
cbind(fit1$grid[,2],fit1$funi[[2]])
# Extracting the bivariate density estimates
fit1$grid[,1]
fit1$grid[,2]
fit1$fbiv[[1]]
# Plot of the estimated density along with the
# true data points and censoring limits
contour(fit1$grid[,1],fit1$grid[,2],fit1$fbiv[[1]])
points(y)
abline(v=xxlim)
abline(h=yylim)Run the code above in your browser using DataLab