Description Usage Arguments Details Value Author(s) References See Also Examples
This function generates a posterior density sample for a Linear Dependent Bernstein-Dirichlet Process model for bounded conditional density estimation.
1 2 3 4 |
formula |
a two-sided linear formula object describing the
model fit, with the response on the
left of a |
xpred |
a matrix giving the covariate values where the predictive density is evaluated. |
prior |
a list giving the prior information. The list includes the following
parameter: |
mcmc |
a list giving the MCMC parameters. The list must include
the following elements: |
state |
a list giving the current value of the parameters. This list is used if the current analysis is the continuation of a previous analysis. |
status |
a logical variable indicating whether this run is new ( |
ngrid |
integer giving the number of grid points where the conditional density estimate is evaluated. The default is 100. |
grid |
vector of grid points where the conditional density estimate is evaluated. The default value is NULL and the grid is chosen according to the range of the data. |
compute.band |
logical variable indicating whether the credible band for the conditional density and mean function must be computed. |
type.band |
string indication the type of credible band to be computed; if equal to "HPD" or "PD" then the 95 percent pointwise HPD or PD band is computed, respectively. |
data |
data frame. |
na.action |
a function that indicates what should happen when the data
contain |
work.dir |
working directory. |
This generic function fits a Linear Dependent Dirichlet-Bernstein model (Barrientos, Jara and Quintana, 2010), given by:
yi | GXi ~ GXi
{GX: X in X} | k,alpha, G0 ~ LDBDP(k,alpha G0)
where, G0 = N(beta| mub, Sb). To complete the model specification, independent hyperpriors are assumed,
k | lambda ~ Poisson(lambda)I(k > 0)
mub | m0, S0 ~ N(m0,S0)
Sb | nu, psi ~ IW(nu,psi)
Note that the inverted-Wishart prior is parametrized such that if A ~ IWq(nu, psi) then E(A)= psiinv/(nu-q-1).
Note also that the LDBDP model is a extension of the Bernstein-Dirichlet model for density estimation (Petrone, 1999a, 1999b; Petrone and Waserman, 2002).
The computational implementation of the model is based on the finite approximation to the dependent Dirichlet process prior and on the use of conditional MCMC methods. The regression coefficients and stick-breaking parameters are updated jointly using multivariate Slice sampling (Neal, 2003). The degree of the Bernstein polynomial is updated using a Metropolis-Hasting algorithm.
An object of class LDBDPdensity
representing the LDBDP model fit.
Generic functions such as print
, plot
,
and summary
have methods to show the results of the fit. The results include
k
. mub
and Sb
.
The list state
in the output object contains the current value of the parameters
necessary to restart the analysis. If you want to specify different starting values
to run multiple chains set status=TRUE
and create the list state based on
this starting values. In this case the list state
must include the following objects:
k |
an integer giving the degree of the Bernstein polynomial. |
beta |
a matrix of dimension |
alpha |
giving the value of the precision parameter. |
mub |
giving the mean of the normal baseline distributions. |
Sb |
giving the covariance matrix the normal baseline distributions. |
v |
giving the |
Felipe Barrientos <afbarrie@.mat.puc.cl>
Alejandro Jara <atjara@uc.cl>
Barrientos, F., Jara, A., Quintana, F. (2010). Bounded density regression using dependent Bernstein polynomials. Technical Report, Department of Statistics, Pontificia Universidad Catolica de Chile.
Neal, R. (2003) Slice sampling. Anals of Statistics, 31: 705-767.
Petrone, S. (1999a) Random Bernstein polynomials. Scandinavian Journal of Statistics, 26: 373-393.
Petrone, S. (1999b) Bayesian density estimation using Bernstein polynomials. The Canadian Journal of Statistics, 27: 105-126.
Petrone, S. and Waserman, L. (2002) Consistency of Bernstein polynomial posterior. Journal of the Royal Statistical Society, Series B, 64: 79-100.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 | ## Not run:
########################
# Simulate data
########################
nrec <- 500
x <- runif(nrec)
y <- rep(0,nrec)
for(i in 1:nrec)
{
y[i] <- ifelse(runif(1) < (0.8-0.5*x[i]^2),
rbeta(1,22-(x[i]^2)*20,5+x[i]*20),
rbeta(1,8+x[i]*5,20))
}
# true model
true.dens <- function(grid,x)
{
(0.8-0.5*x^2)*dbeta(grid,22-(x^2)*20,5+x*20)+
(0.2+0.5*x^2)*dbeta(grid,8+x*5,20)
}
true.mean <- function(x)
{
(0.8-0.5*x^2)*(22-(x^2)*20)/(22-(x^2)*20+5+x*20)+
(0.2+0.5*x^2)*(8+x*5)/(8+x*5+20)
}
# predictions
grid <- seq(0,1,len=100)
npred <- 25
xpred <- matrix(1,ncol=2,nrow=npred)
xpred[,2] <- seq(0,1,len=npred)
# prior
prior <- list(maxn = 25,
alpha = 1,
lambda = 25,
nu = 4,
psiinv = diag(1000,2),
m0 = rep(0,2),
S0 = diag(1000,2))
# mcmc
mcmc <- list(nburn = 5000,
nskip = 3,
ndisplay = 100,
nsave = 5000)
# state
state <- NULL
# fitting the model
fit <- LDBDPdensity(formula=y~x,xpred=xpred,
prior=prior,
mcmc=mcmc,
state=NULL,status=TRUE,
grid=grid,
compute.band=TRUE,type.band="PD")
fit
summary(fit)
plot(fit)
# Plots for some predictions
# (conditional density and mean function)
par(mfrow=c(2,2))
plot(fit$grid,fit$densp.h[7,],type="l",lwd=2,
xlim=c(0,1),ylim=c(0,4),xlab="y",ylab="density",lty=2)
lines(fit$grid,fit$densp.m[7,],lwd=2,lty=1)
lines(fit$grid,fit$densp.l[7,],lwd=2,lty=2)
lines(fit$grid,true.dens(fit$grid,fit$xpred[7,2]),lwd=2,col="red")
plot(fit$grid,fit$densp.h[13,],type="l",lwd=2,
xlim=c(0,1),ylim=c(0,4),xlab="y",ylab="density",lty=2)
lines(fit$grid,fit$densp.m[13,],lwd=2,lty=1)
lines(fit$grid,fit$densp.l[13,],lwd=2,lty=2)
lines(fit$grid,true.dens(fit$grid,fit$xpred[13,2]),lwd=2,col="red")
plot(fit$grid,fit$densp.h[19,],type="l",lwd=2,
xlim=c(0,1),ylim=c(0,4),xlab="y",ylab="density",lty=2)
lines(fit$grid,fit$densp.m[19,],lwd=2,lty=1)
lines(fit$grid,fit$densp.l[19,],lwd=2,lty=2)
lines(fit$grid,true.dens(fit$grid,fit$xpred[19,2]),lwd=2,col="red")
plot(x,y)
lines(fit$xpred[,2],fit$meanfp.m,lwd=2,lty=1)
lines(fit$xpred[,2],fit$meanfp.l,lwd=2,lty=2)
lines(fit$xpred[,2],fit$meanfp.h,lwd=2,lty=2)
lines(fit$xpred[,2],true.mean(fit$xpred[,2]),lwd=2,lty=1,col="red")
## End(Not run)
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.