if (FALSE) {
# EXAMPLE 1 (INTERFACE=FORMULA): For this example, we load Giovanni
# Baiocchi's Italian GDP panel (see Italy for details), and compute the
# likelihood cross-validated bandwidths (default) using a second-order
# Gaussian kernel (default). Note - this may take a minute or two
# depending on the speed of your computer.
data("Italy")
attach(Italy)
# First, compute the bandwidths... note that this may take a minute or
# two depending on the speed of your computer.
bw <- npcdensbw(formula=gdp~ordered(year))
# Next, compute the condensity object...
fhat <- npcdens(bws=bw)
# The object fhat now contains results such as the estimated conditional
# density function (fhat$condens) and so on...
summary(fhat)
# Call the plot() function to visualize the results (-C will
# interrupt on *NIX systems, will interrupt on MS Windows
# systems).
plot(bw)
detach(Italy)
# EXAMPLE 1 (INTERFACE=DATA FRAME): For this example, we load Giovanni
# Baiocchi's Italian GDP panel (see Italy for details), and compute the
# likelihood cross-validated bandwidths (default) using a second-order
# Gaussian kernel (default). Note - this may take a minute or two
# depending on the speed of your computer.
data("Italy")
attach(Italy)
# First, compute the bandwidths... note that this may take a minute or
# two depending on the speed of your computer.
# Note - we cast `X' and `y' as data frames so that plot() can
# automatically grab names (this looks like overkill, but in
# multivariate settings you would do this anyway, so may as well get in
# the habit).
X <- data.frame(year=ordered(year))
y <- data.frame(gdp)
bw <- npcdensbw(xdat=X, ydat=y)
# Next, compute the condensity object...
fhat <- npcdens(bws=bw)
# The object fhat now contains results such as the estimated conditional
# density function (fhat$condens) and so on...
summary(fhat)
# Call the plot() function to visualize the results (-C will
# interrupt on *NIX systems, will interrupt on MS Windows systems).
plot(bw)
detach(Italy)
# EXAMPLE 2 (INTERFACE=FORMULA): For this example, we load the old
# faithful geyser data from the R `datasets' library and compute the
# conditional density function.
library("datasets")
data("faithful")
attach(faithful)
# Note - this may take a few minutes depending on the speed of your
# computer...
bw <- npcdensbw(formula=eruptions~waiting)
summary(bw)
# Plot the density function (-C will interrupt on *NIX systems,
# will interrupt on MS Windows systems).
plot(bw)
detach(faithful)
# EXAMPLE 2 (INTERFACE=DATA FRAME): For this example, we load the old
# faithful geyser data from the R `datasets' library and compute the
# conditional density function.
library("datasets")
data("faithful")
attach(faithful)
# Note - this may take a few minutes depending on the speed of your
# computer...
# Note - we cast `X' and `y' as data frames so that plot() can
# automatically grab names (this looks like overkill, but in
# multivariate settings you would do this anyway, so may as well get in
# the habit).
X <- data.frame(waiting)
y <- data.frame(eruptions)
bw <- npcdensbw(xdat=X, ydat=y)
summary(bw)
# Plot the density function (-C will interrupt on *NIX systems,
# will interrupt on MS Windows systems)
plot(bw)
detach(faithful)
# EXAMPLE 3 (INTERFACE=FORMULA): Replicate the DGP of Klein & Spady
# (1993) (see their description on page 405, pay careful attention to
# footnote 6 on page 405).
set.seed(123)
n <- 1000
# x1 is chi-squared having 3 df truncated at 6 standardized by
# subtracting 2.348 and dividing by 1.511
x <- rchisq(n, df=3)
x1 <- (ifelse(x < 6, x, 6) - 2.348)/1.511
# x2 is normal (0, 1) truncated at +- 2 divided by 0.8796
x <- rnorm(n)
x2 <- ifelse(abs(x) < 2 , x, 2) / 0.8796
# y is 1 if y* > 0, 0 otherwise.
y <- ifelse(x1 + x2 + rnorm(n) > 0, 1, 0)
# Generate data-driven bandwidths (likelihood cross-validation). Note -
# this may take a few minutes depending on the speed of your computer...
bw <- npcdensbw(formula=factor(y)~x1+x2)
# Next, create the evaluation data in order to generate a perspective
# plot
x1.seq <- seq(min(x1), max(x1), length=50)
x2.seq <- seq(min(x2), max(x2), length=50)
X.eval <- expand.grid(x1=x1.seq,x2=x2.seq)
data.eval <- data.frame(y=factor(rep(1, nrow(X.eval))),x1=X.eval[,1],x2=X.eval[,2])
# Now evaluate the conditional probability for y=1 and for the
# evaluation Xs
fit <- fitted(npcdens(bws=bw,newdata=data.eval))
# Finally, coerce the data into a matrix for plotting with persp()
fit.mat <- matrix(fit, 50, 50)
# Generate a perspective plot similar to Figure 2 b of Klein and Spady
# (1993)
persp(x1.seq,
x2.seq,
fit.mat,
col="white",
ticktype="detailed",
expand=0.5,
axes=FALSE,
box=FALSE,
main="Estimated Nonparametric Probability Perspective",
theta=310,
phi=25)
# EXAMPLE 3 (INTERFACE=DATA FRAME): Replicate the DGP of Klein & Spady
# (1993) (see their description on page 405, pay careful attention to
# footnote 6 on page 405).
set.seed(123)
n <- 1000
# x1 is chi-squared having 3 df truncated at 6 standardized by
# subtracting 2.348 and dividing by 1.511
x <- rchisq(n, df=3)
x1 <- (ifelse(x < 6, x, 6) - 2.348)/1.511
# x2 is normal (0, 1) truncated at +- 2 divided by 0.8796
x <- rnorm(n)
x2 <- ifelse(abs(x) < 2 , x, 2) / 0.8796
# y is 1 if y* > 0, 0 otherwise.
y <- ifelse(x1 + x2 + rnorm(n) > 0, 1, 0)
# Create the X matrix
X <- cbind(x1, x2)
# Generate data-driven bandwidths (likelihood cross-validation). Note -
# this may take a few minutes depending on the speed of your computer...
bw <- npcdensbw(xdat=X, ydat=factor(y))
# Next, create the evaluation data in order to generate a perspective
# plot
x1.seq <- seq(min(x1), max(x1), length=50)
x2.seq <- seq(min(x2), max(x2), length=50)
X.eval <- expand.grid(x1=x1.seq,x2=x2.seq)
# Now evaluate the conditional probability for y=1 and for the
# evaluation Xs
fit <- fitted(npcdens(exdat=X.eval,
eydat=factor(rep(1, nrow(X.eval))),
bws=bw))
# Finally, coerce the data into a matrix for plotting with persp()
fit.mat <- matrix(fit, 50, 50)
# Generate a perspective plot similar to Figure 2 b of Klein and Spady
# (1993)
persp(x1.seq,
x2.seq,
fit.mat,
col="white",
ticktype="detailed",
expand=0.5,
axes=FALSE,
box=FALSE,
main="Estimated Nonparametric Probability Perspective",
theta=310,
phi=25)
}
Run the code above in your browser using DataLab