MixNRMI2: Normalized Random Measures Mixture of Type II

Description Usage Arguments Details Value Warning Author(s) References See Also Examples

Description

Bayesian nonparametric estimation based on normalized measures driven mixtures for locations and scales.

Usage

1
2
3
4
MixNRMI2(x, probs = c(0.025, 0.5, 0.975), Alpha = 1, Beta = 0, 
    Gama = 0.4, distr.k = 1, distr.py0 = 1, distr.pz0 = 2, mu.pz0 = 3, 
    sigma.pz0 = sqrt(10), delta = 4, kappa = 2, Delta = 2, Meps = 0.01, 
    Nx = 150, Nit = 1500, Pbi = 0.1, epsilon = NULL, printtime = TRUE, extras = FALSE) 

Arguments

x

Numeric vector. Data set to which the density is fitted.

probs

Numeric vector. Desired quantiles of the density estimates.

Alpha

Numeric constant. Total mass of the centering measure. See details.

Beta

Numeric positive constant. See details.

Gama

Numeric constant. 0 <= Gama <=1. See details.

distr.k

Integer number identifying the mixture kernel: 1 = Normal; 2 = Gamma; 3 = Beta; 4 = Double Exponential; 5 = Lognormal.

distr.py0

Integer number identifying the centering measure for locations: 1 = Normal; 2 = Gamma; 3 = Beta.

distr.pz0

Integer number identifying the centering measure for scales: 2 = Gamma. For more options use MixNRMI2cens.

mu.pz0

Numeric constant. Prior mean of the centering measure for scales.

sigma.pz0

Numeric constant. Prior standard deviation of the centering measure for scales.

delta

Numeric positive constant. Metropolis-Hastings proposal variation coefficient for sampling the scales.

kappa

Numeric positive constant. Metropolis-Hastings proposal variation coefficient for sampling the location parameters.

Delta

Numeric positive constant. Metropolis-Hastings proposal variation coefficient for sampling the latent U.

Meps

Numeric constant. Relative error of the jump sizes in the continuous component of the process. Smaller values imply larger number of jumps.

Nx

Integer constant. Number of grid points for the evaluation of the density estimate.

Nit

Integer constant. Number of MCMC iterations.

Pbi

Numeric constant. Burn-in period proportion of Nit.

epsilon

Numeric constant. Extension to the evaluation grid range. See details.

printtime

Logical. If TRUE, prints out the execution time.

extras

Logical. If TRUE, gives additional objects: means, sigmas, weights and Js.

Details

This generic function fits a normalized random measure (NRMI) mixture model for density estimation (James et al. 2009). Specifically, the model assumes a normalized generalized gamma (NGG) prior for both, locations (means) and standard deviations, of the mixture kernel, leading to a fully nonparametric mixture model.

The details of the model are:

X_i|Y_i,Z_i ~ k(.|Y_i,Z_i)

(Y_i,Z_i)|P ~ P, i=1,...,n

P ~ NGG(Alpha, Beta, Gama; P_0)

where, X_i's are the observed data, (Y_i,Z_i)'s are bivariate latent (location and scale) vectors, k is a parametric kernel parameterized in terms of mean and standard deviation, (Alpha, Beta, Gama; P_0) are the parameters of the NGG prior with a bivariate P_0 being the centering measure with independent components, that is, P_0(Y,Z) = P_0(Y)*P_0(Z). The parameters of P_0(Y) are assigned vague hyper prior distributions and (mu.pz0,sigma.pz0) are the hyper-parameters of P_0(Z). In particular, NGG(Alpha, 1, 0; P_0) defines a Dirichlet process; NGG(1, Beta, 1/2;P_0) defines a Normalized inverse Gaussian process; and NGG(1, 0, Gama; P_0) defines a normalized stable process. The evaluation grid ranges from min(x) - epsilon to max(x) + epsilon. By default epsilon=sd(x)/4.

Value

The function returns a list with the following components:

xx

Numeric vector. Evaluation grid.

qx

Numeric array. Matrix of dimension Nx x (length(probs)+1) with the posterior mean and the desired quantiles input in probs.

cpo

Numeric vector of length(x) with conditional predictive ordinates.

R

Numeric vector of length(Nit*(1-Pbi)) with the number of mixtures components (clusters).

U

Numeric vector of length(Nit*(1-Pbi)) with the values of the latent variable U.

Allocs

List of length(Nit*(1-Pbi)) with the clustering allocations.

means

List of length(Nit*(1-Pbi)) with the cluster means (locations). Only if extras = TRUE.

sigmas

Numeric vector of length(Nit*(1-Pbi)) with the cluster standard deviations. Only if extras = TRUE.

weights

List of length(Nit*(1-Pbi)) with the mixture weights. Only if extras = TRUE.

Js

List of length(Nit*(1-Pbi)) with the unnormalized weights (jump sizes). Only if extras = TRUE.

Nm

Integer constant. Number of jumps of the continuous component of the unnormalized process.

Nx

Integer constant. Number of grid points for the evaluation of the density estimate.

Nit

Integer constant. Number of MCMC iterations.

Pbi

Numeric constant. Burn-in period proportion of Nit.

procTime

Numeric vector with execution time provided by proc.time function.

Warning

The function is computing intensive. Be patient.

Author(s)

Barrios, E., Lijoi, A., Nieto-Barajas, L.E. and Prüenster, I.

References

1.- Barrios, E., Lijoi, A., Nieto-Barajas, L. E. and Prüenster, I. (2013). Modeling with Normalized Random Measure Mixture Models. Statistical Science. Vol. 28, No. 3, 313-334.

2.- James, L.F., Lijoi, A. and Prüenster, I. (2009). Posterior analysis for normalized random measure with independent increments. Scand. J. Statist 36, 76-97.

See Also

MixNRMI1, MixNRMI1cens, MixNRMI2cens

Examples

 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
96
97
98
99
## Not run: 
### Example 1
# Data
data(acidity)
x <- acidity
# Fitting the model under default specifications
out <- MixNRMI2(x)
# Plotting density estimate + 95
attach(out)
m <- ncol(qx)
ymax <- max(qx[,m])
par(mfrow=c(1,1))
hist(x,probability=TRUE,breaks=20,col=grey(.9),ylim=c(0,ymax))
lines(xx,qx[,1],lwd=2)
lines(xx,qx[,2],lty=3,col=4)
lines(xx,qx[,m],lty=3,col=4)
detach()

## End(Not run)

### Example 2
## Do not run
# set.seed(150520)
# data(enzyme)
# x <- enzyme
#  Enzyme2.out <- MixNRMI2(x, Alpha = 1, Beta = 0.007, Gama = 0.5, 
#                          distr.k = 2, distr.py0 = 2,
#                          distr.pz0 = 2, mu.pz0 = 1, sigma.pz0 = 1, Meps=0.005,
#                          Nit = 5000, Pbi = 0.2)
# The output of this run is already loaded in the package
# To show results run the following
# Data
data(enzyme)
x <- enzyme
data(Enzyme2.out)
attach(Enzyme2.out)
# Plotting density estimate + 95% credible interval
m <- ncol(qx)
ymax <- max(qx[,m])
par(mfrow=c(1,1))
hist(x,probability=TRUE,breaks=20,col=grey(.9),ylim=c(0,ymax))
lines(xx,qx[,1],lwd=2)
lines(xx,qx[,2],lty=3,col=4)
lines(xx,qx[,m],lty=3,col=4)
# Plotting number of clusters
par(mfrow=c(2,1))
plot(R,type="l",main="Trace of R")
hist(R,breaks=min(R-0.5):max(R+0.5),probability=TRUE)
# Plotting u
par(mfrow=c(2,1))
plot(U,type="l",main="Trace of U")
hist(U,nclass=20,probability=TRUE,main="Histogram of U")
# Plotting cpo
par(mfrow=c(2,1))
plot(cpo,main="Scatter plot of CPO's")
boxplot(cpo,horizontal=TRUE,main="Boxplot of CPO's")
print(paste('Average log(CPO)=',round(mean(log(cpo)),4)))
print(paste('Median log(CPO)=',round(median(log(cpo)),4)))
detach()

### Example 3
## Do not run
# set.seed(150520)
# data(galaxy)
# x <- galaxy
#  Galaxy2.out <- MixNRMI2(x, Alpha = 1, Beta = 0.015, Gama = 0.5, 
#                          distr.k = 1, distr.py0 = 2,
#                          distr.pz0 = 2, mu.pz0 = 1, sigma.pz0 = 1,  Meps=0.005,
#                          Nit = 5000, Pbi = 0.2)
# The output of this run is already loaded in the package
# To show results run the following
# Data
data(galaxy)
x <- galaxy
data(Galaxy2.out)
attach(Galaxy2.out)
# Plotting density estimate + 95% credible interval
m <- ncol(qx)
ymax <- max(qx[,m])
par(mfrow=c(1,1))
hist(x,probability=TRUE,breaks=20,col=grey(.9),ylim=c(0,ymax))
lines(xx,qx[,1],lwd=2)
lines(xx,qx[,2],lty=3,col=4)
lines(xx,qx[,m],lty=3,col=4)
# Plotting number of clusters
par(mfrow=c(2,1))
plot(R,type="l",main="Trace of R")
hist(R,breaks=min(R-0.5):max(R+0.5),probability=TRUE)
# Plotting u
par(mfrow=c(2,1))
plot(U,type="l",main="Trace of U")
hist(U,nclass=20,probability=TRUE,main="Histogram of U")
# Plotting cpo
par(mfrow=c(2,1))
plot(cpo,main="Scatter plot of CPO's")
boxplot(cpo,horizontal=TRUE,main="Boxplot of CPO's")
print(paste('Average log(CPO)=',round(mean(log(cpo)),4)))
print(paste('Median log(CPO)=',round(median(log(cpo)),4)))
detach()


Search within the BNPdensity package
Search all R packages, documentation and source code

Questions? Problems? Suggestions? or email at ian@mutexlabs.com.

Please suggest features or report bugs with the GitHub issue tracker.

All documentation is copyright its authors; we didn't write any of that.