许多书籍通过以下图示来说明Fisher线性判别分析的概念（这张图来自于模式识别与机器学习，第188页）

我想知道如何在R语言（或其他任何语言）中重现这个图。我在R中进行了初步尝试，模拟了两组数据，并使用abline()函数绘制了线性判别。欢迎任何建议。

set.seed(2014)library(MASS)library(DiscriMiner) # For scatter matrices# 模拟二元正态分布，包含2个类别mu1 <- c(2, -4)mu2 <- c(2, 6)rho <- 0.8s1 <- 1s2 <- 3Sigma <- matrix(c(s1^2, rho * s1 * s2, rho * s1 * s2, s2^2), byrow = TRUE, nrow = 2)n <- 50X1 <- mvrnorm(n, mu = mu1, Sigma = Sigma)X2 <- mvrnorm(n, mu = mu2, Sigma = Sigma)y <- rep(c(0, 1), each = n)X <- rbind(x1 = X1, x2 = X2)X <- scale(X)# 散布矩阵B <- betweenCov(variables = X, group = y)W <- withinCov(variables = X, group = y)# 特征向量ev <- eigen(solve(W) %*% B)$vectorsslope <- - ev[1,1] / ev[2,1]intercept <- ev[2,1]par(pty = "s")plot(X, col = y + 1, pch = 16)abline(a = slope, b = intercept, lwd = 2, lty = 2)

我的（未完成的）工作

我将当前的解决方案粘贴在下面。主要问题是如何根据决策边界旋转（和移动）密度图。仍然欢迎任何建议。

require(ggplot2)library(grid)library(MASS)# 模拟参数mu1 <- c(5, -9)mu2 <- c(4, 9)rho <- 0.5s1 <- 1s2 <- 3Sigma <- matrix(c(s1^2, rho * s1 * s2, rho * s1 * s2, s2^2), byrow = TRUE, nrow = 2)n <- 50# 多元正态抽样X1 <- mvrnorm(n, mu = mu1, Sigma = Sigma)X2 <- mvrnorm(n, mu = mu2, Sigma = Sigma)# 组合成数据框y <- rep(c(0, 1), each = n)X <- rbind(x1 = X1, x2 = X2)X <- scale(X)X <- data.frame(X, class = y)# 应用lda()m1 <- lda(class ~ X1 + X2, data = X)m1.pred <- predict(m1)# 计算abline的截距和斜率gmean <- m1$prior %*% m1$meansconst <- as.numeric(gmean %*% m1$scaling)z <- as.matrix(X[, 1:2]) %*% m1$scaling - constslope <- - m1$scaling[1] / m1$scaling[2]intercept <- const / m1$scaling[2]# 投影值LD <- data.frame(predict(m1)$x, class = y)# 散点图p1 <- ggplot(X, aes(X1, X2, color=as.factor(class))) +   geom_point() +  theme_bw() +  theme(legend.position = "none") +  scale_x_continuous(limits=c(-5, 5)) +   scale_y_continuous(limits=c(-5, 5)) +  geom_abline(intecept = intercept, slope = slope)# 密度图 p2 <- ggplot(LD, aes(x = LD1)) +  geom_density(aes(fill = as.factor(class), y = ..scaled..)) +  theme_bw() +  theme(legend.position = "none")grid.newpage()print(p1)vp <- viewport(width = .7, height = 0.6, x = 0.5, y = 0.3, just = c("centre"))pushViewport(vp)print(p2, vp = vp)

回答：

基本上，你需要沿着分类器的方向投影数据，为每个类别绘制直方图，然后旋转直方图，使其x轴与分类器平行。为了获得良好的结果，需要对直方图进行一些试错调整。这里有一个在Matlab中如何做的示例，适用于朴素分类器（类别均值的差异）。对于Fisher分类器，过程当然是相似的，你只需使用不同的分类器w。我更改了你代码中的参数，使图表更接近你提供的图表。

rng('default')n = 1000;mu1 = [1,3]';mu2 = [4,1]';rho = 0.3;s1 = .8;s2  = .5;Sigma = [s1^2,rho*s1*s1;rho*s1*s1, s2^2];X1 = mvnrnd(mu1,Sigma,n);X2 = mvnrnd(mu2,Sigma,n);X = [X1; X2];Y = [zeros(n,1);ones(n,1)];scatter(X1(:,1), X1(:,2), [], 'b' );hold onscatter(X2(:,1), X2(:,2), [], 'r' );axis equalm1 = mean(X(1:n,:))';m2 = mean(X(n+1:end,:))';plot(m1(1),m1(2),'bx','markersize',18)plot(m2(1),m2(2),'rx','markersize',18)plot([m1(1),m2(1)], [m1(2),m2(2)],'g')%% 仅考虑均值的分类器w = m2 - m1; w = w / norm(w);% 将数据投影到w上X1_projected = X1 * w;X2_projected = X2 * w;% 绘制直方图并旋转它angle = 180/pi * atan(w(2)/w(1));[hy1, hx1] = hist(X1_projected);[hy2, hx2] = hist(X2_projected);hy1 = hy1 / sum(hy1); % 归一化hy2 = hy2 / sum(hy2); % 归一化scale = 4; % 手动设置h1 = bar(hx1, scale*hy1,'b');h2 = bar(hx2, scale*hy2,'r');set([h1, h2],'ShowBaseLine','off')% 围绕原点旋转rotate(get(h1,'children'),[0,0,1], angle, [0,0,0])rotate(get(h2,'children'),[0,0,1], angle, [0,0,0])