1

我有一些编码问题，在线性判别分析中做一些练习。我们正在使用 Iris 数据：

## Read in dataset, set seed, load package
Iris <- iris[,-(1:2)]
grIris <- as.integer(iris[,"Species"])
set.seed(16)
library(MASS)

## Read n
n <- nrow(Iris)

如您所见，我们删除了 iris 的第一列和第二列。我想要做的是使用线性判别分析来引导这些数据，这是我的代码：

ind <- replicate(B,sample(seq(1:n),n,replace=TRUE))

这会生成我想要使用的索引。注意B是一些很大的数字，例如1000。现在我想使用apply，但是为什么下面的代码不起作用？

bst.sample <- apply(ind,2,lda(Species~Petal.Length+Petal.Width,data=Iris[ind,]))

其中 Species、Petal.Length 等是来自 iris 的数据。如果我使用 for 循环，一切正常，但我当然想以这种更优雅的方式实现。

我的第二个问题是关于points. 我还想计算估计的平均值，我通过以下代码完成了

est.lda <- vector("list",B)
est.qda <- vector("list",B)
mu_hat_1 <- mu_hat_2 <- mu_hat_3 <- matrix(0,ncol=B,nrow=2)
for (i in 1:B){
  est.lda[[i]] <- lda(Species~Petal.Length+Petal.Width,data=Iris[ind[,i],])
  mu_hat_1[,i] <- est.lda[[i]]$means[1,]
  mu_hat_2[,i] <- est.lda[[i]]$means[2,]
  mu_hat_3[,i] <- est.lda[[i]]$means[3,]
  est.qda[[i]] <- qda(Species~Petal.Length+Petal.Width,data=Iris[ind[,i],])

}

plot(mu_hat_1[1,],mu_hat_1[2,],pch=4)
points(mu_hat_2[1,],mu_hat_2[2,],pch=4,col=2)
points(mu_hat_3[1,],mu_hat_3[2,],pch=4,col=3)

最后的图应该显示三个区域，三个区域的预期平均值。然而，只显示了第一个图。

感谢您的帮助。

4

1 回答 1

2 
B <- 10
ind <- replicate(B,sample(seq(1:n),n,replace=TRUE))

#you need to pass a function to apply
bst.sample <- apply(ind,2, 
                function(i) lda(Species~Petal.Length+Petal.Width,data=Iris[i,]))
#extract means
bst.means <- lapply(bst.sample,function(x) x$means)

#bind means into array
library(abind)
bst.means <- do.call(function(...) abind(..., along=3), bst.means)

#you need to make sure that alle points are inside the axis limits
plot(bst.means[1,1,],bst.means[1,2,], 
     xlim=range(bst.means[,1,]), ylim=range(bst.means[,2,]), 
     xlab=dimnames(bst.means)[[2]][1],ylab=dimnames(bst.means)[[2]][2],
     col=1)
points(bst.means[2,1,],bst.means[2,2,], col=2)
points(bst.means[3,1,],bst.means[3,2,], col=3)
legend("topleft", legend=dimnames(bst.means)[[1]], col=1:3, pch=1)

在此处输入图像描述

于 2013-08-04T12:02:50.660 回答