1

線形判別分析でいくつかの演習を行っているときに発生する、コーディングに関する質問があります。アイリスデータを使用しています：

## 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)

ご覧のとおり、虹彩の 1 列目と 2 列目を削除します。私がやりたいことは、線形判別分析を使用したこのデータのブートストラップです。ここに私のコードがあります:

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

これにより、使用したいインデックスが生成されます。1000 などの大きな数に注意してくださいB。apply を使用したいのですが、次のコードが機能しないのはなぜですか?

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

ここで、Species、Petal.Length などはアイリスのデータです。for ループを使用すると、すべて正常に動作しますが、もちろん、このよりエレガントな方法で実装したいと考えています。

2 番目の質問は について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)

最後のプロットには、3 つのクラスの予想平均を含む 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 に答える