plot.plda <- function(result, Xnew, y, n.seq = 100, threshold = 5, ...) {
N <- nrow(Xnew)
P <- ncol(Xnew)
is.pca <- P != 2
stopifnot(result$parameters$P == P)
# Use first 2 PC's to plot
if (is.pca) {
pca <- prcomp(Xnew, center = TRUE, scale. = TRUE)
X.pca <- cbind(pca$x[, 1], pca$x[, 2])
fit.pca <- plda(X.pca, y, type = result$parameters$type, prior = result$parameters$prior)
X <- X.pca
model <- fit.pca
} else {
X <- Xnew
model <- result
}
# Generate boundries
Xmin <- min(X[, 1]) - threshold
Xmax <- max(X[, 1]) + threshold
Ymin <- min(X[, 2]) - threshold
Ymax <- max(X[, 2]) + threshold
# Generate grid
X.grid <- point.grid(Xmin, Xmax, Ymin, Ymax, n.seq)
y.grid <- predict(model, X.grid)
# Dataframes
df.data <- data.frame(x1 = X[, 1], x2 = X[, 2], y = y)
df.grid <- data.frame(x1 = X.grid[, 1], x2 = X.grid[, 2], y = y.grid)
# Plot deciscion boundry
ggplot2::ggplot() +
ggplot2::geom_jitter(data = df.data, ggplot2::aes(x = x1, y = x2, color = y)) +
ggplot2::geom_point(data = df.grid, ggplot2::aes(x = x1, y = x2, color = y), alpha = 0.2) +
ggplot2::theme_bw() +
ggplot2::xlab("feature 1") +
ggplot2::ylab("feature 2") +
ggplot2::xlim(Xmin, Xmax) +
ggplot2::ylim(Ymin, Ymax) +
ggplot2::scale_x_continuous(expand = c(0, 0)) +
ggplot2::scale_y_continuous(expand = c(0, 0)) +
ggplot2::theme(legend.position = "none")
#ggplot2::theme(legend.title = ggplot2::element_blank())
}
predict.plda <- function(result, Xnew, ...) {
stopifnot(is.matrix(Xnew), ncol(Xnew) == result$parameters$P)
stopifnot(class(result) == "plda")
if (result$parameters$type %in% c("linear", "quadratic")) {
densities <- estimate.normal.densities(
Xnew
, mu.hat = result$estimates$mu.hat
, sigma.hat = result$estimates$sigma.hat
, type = result$parameters$type
)
}
if (result$parameters$type == "poisson") {
densities <- estimate.poisson.densities(
Xnew
, lambda.hat = result$estimate$lambda.hat
)
}
if (result$parameters$type == "poisson.seq") {
s.estimate <- estimate.size.factor(
Xnew
, X.train.parameter = result$estimates$X.train.parameter
, type = result$parameters$size.factor
)
s.hat <- s.estimate$s.hat
N.hat <- matrix(s.hat) %*% matrix(result$estimates$g.hat, nrow = 1)
densities <- estimate.poisson.seq.densities(
Xnew
, N.hat = N.hat
, d.hat = result$estimates$d.hat
)
}
posteriors <- bayes.rule(Xnew, densities, result$estimates$pi.hat)
fitted.posteriors <- fit.posteriors(posteriors, result$parameters$levels)
return(fitted.posteriors)
}
print.plda <- function(result, ...) {
dist <- ifelse(result$parameters$type == "poisson", "poisson", "normal")
cat(sprintf("X(%i x %i) with K = %i\n", result$parameters$N, result$parameters$P, result$parameters$K))
cat(sprintf("Type: **%s** discriminant analysis", result$parameters$type), "\n")
cat(sprintf("Dist: predictors have a **%s** distribution", dist), "\n")
cat("\npi.hat:\n")
printCoefmat(t(as.matrix(result$estimates$pi.hat)))
if (dist == "normal") {
cat("\nmu.hat:\n")
printCoefmat(result$estimates$mu.hat)
}
if (result$parameters$type == "linear") {
cat("\nsigma.hat:\n")
printCoefmat(result$estimates$sigma.hat)
}
if (dist == "poisson") {
cat("\nN.hat (N by p, summarized by column averages):\n")
means <- t(as.matrix(colMeans(as.matrix(result$estimates$N.hat))))
printCoefmat(means, has.Pvalue = FALSE)
cat("\nd.hat:\n")
printCoefmat(result$estimates$d.hat, has.Pvalue = FALSE)
}
cat("\n")
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.