Statystyka, prognozowanie, ekonometria, data mining

bulva · Chorąży **Pomógł:** 5 razy Posty: 172 Skąd: Zgierz

Witam,
Chciałbym w R wykonać PCA w wersji uodpornionej na wartości odstające, w tym celu posługując się macierzą wariancji-kowariancji estymowana metodą MCD. Coś mi tu nie gra, ponieważ zwrócone w ten sposób wyniki nie różnią się niczym od klasycznej PCA. Oto mój kod:

Kod:

#Robust PCA
library(MASS)
mcd <- cov.rob(possum[,6:14],cor=T, nsamp="best", method="mcd")
rob.pca <- princomp(possum[,6:14], covlist=mcd, cor=T, na.action=na.exclude)
summary(rob.pca)

#Zwraca:
#Importance of components:
# Comp.1 Comp.2 Comp.3 Comp.4 Comp.5
#Standard deviation 1.9858729 1.4068546 0.94399072 0.85915143 0.75083196
#Proportion of Variance 0.4381879 0.2199155 0.09901316 0.08201569 0.06263874
#Cumulative Proportion 0.4381879 0.6581034 0.75711661 0.83913229 0.90177103
# Comp.6 Comp.7 Comp.8 Comp.9
#Standard deviation 0.55625841 0.51648331 0.40478944 0.37950998
#Proportion of Variance 0.03438038 0.02963944 0.01820605 0.01600309
#Cumulative Proportion 0.93615141 0.96579085 0.98399691 1.00000000

#Klasyczne PCA:

class.pca <- princomp(possum[,6:14],cor=T, na.action=na.exclude)
summary(class.pca)

#Zwraca
#Importance of components:
# Comp.1 Comp.2 Comp.3 Comp.4 Comp.5
#Standard deviation 1.9858729 1.4068546 0.94399072 0.85915143 0.75083196
#Proportion of Variance 0.4381879 0.2199155 0.09901316 0.08201569 0.06263874
#Cumulative Proportion 0.4381879 0.6581034 0.75711661 0.83913229 0.90177103
# Comp.6 Comp.7 Comp.8 Comp.9
#Standard deviation 0.55625841 0.51648331 0.40478944 0.37950998
#Proportion of Variance 0.03438038 0.02963944 0.01820605 0.01600309
#Cumulative Proportion 0.93615141 0.96579085 0.98399691 1.00000000

Google

Crunchy · Wysłany: 2010-02-27, 11:22

Rozumiem, że zapomniałeś dodać.

Kod:

library(DAAG)

Na windowsowej wersji zwraca błąd, bo nawet przy wybraniu kolumn zostaje row.names, zobacz jak jest u Ciebie, chyba że używasz innego zbioru. MCD jest jeszcze w kilku innych pakietach.

bulva · Chorąży **Pomógł:** 5 razy Posty: 172 Skąd: Zgierz

To sa troszke inne dane, fragment oryginalnego zbioru possum - np sa tylko wie lokacje. Doszedlem juz o co chodzi, taki kod dziala:

Kod:

library(MASS)
mcd <- cov.rob(possum[,6:14],cor=T, nsamp="best", method="mcd")
rob.pca <- princomp(covmat=mcd, cor=T, na.action=na.exclude)

Wydaje mi sie ze w R jest maly balagan z PCA, kazdy pakiet uzywa innego algorytmu a help ubogi. Analiza tym pakietem:

Kod:

library(rrcov)

pcaR <- PcaCov(possum[,6:14], corr=T) #CovMcd
summary(pcaR)

pcaC <- PcaClassic(possum[,6:14],corr=t)
summary(pcaC)

par(mfrow=c(1,2))
biplot(pcaC, main="Classic", xlabs=possum$site)
biplot(pcaR, main="Robust", xlabs=possum$site)
par(mfrow=c(1,1))

Daje kompletnie inne wyniki, dopiero przez eksperymenty doszedlem iz uzywa on dekompozycji SVD (jak funkcja prcomp) a nie standardowych wartosci wlasnych macierzy kowariancji (jak princomp) - w help ani slowa. To samo jest z MCD, wyniki roznia sie miedzy pakietami (nawet pomimo ustawienia seed przed analizami). Pakiet MASS jest "znany i uznany", zdaje sie ze nawet jest w podstawowej instalacji R, wiec trzymam sie jego.

Druga sprawa, to wykresy diagnostyczne, ani prcomp ani princomp nie maja domyslnej funkcji, w innych pakietach sa ale znowu obszary krytyczne sie roznia, narazie mam to:

Kod:

diagn.plot.pca <- function (pca.data, pca.object, title, a=2)
{
if (is.null(a))
{
a = ncol(pca.object$loadings)
}
Score.Dist = sqrt(apply(t(t(pca.object$scores[, 1:a]^2)/pca.object$sdev[1:a]^2), 1, sum))
Orth.Dist = sqrt(apply((pca.data - pca.object$scores[, 1:a] %*% t(pca.object$loadings[, 1:a]))^2, 1, sum))

#z google, sa inne:
cutoff.SD = sqrt(qchisq(0.975, a))
cutoff.OD = (median(Orth.Dist^(2/3)) + mad(Orth.Dist^(2/3)) * qnorm(0.975))^(3/2)

plot(Score.Dist, Orth.Dist, xlab=paste("Score distance",a,"PC's" ), ylab=paste("Orthogonal Distance",a,"PC's"),
col = "blue", main=title)
abline(h=cutoff.OD, col="red")
abline(v=cutoff.SD, col="red")
}

Trzecia sprawa to iż ciężko coś wyciągnąc z przestrzeni nazw tych funkcji. Lubie sobie robic sliczne tabelki (polecam swoją drogą, naprawde przyspiesza i ulatwia prace), ktore potem wstawiam do latexa, natomiast ten kod zamiast wariancji, uparcie zwraca mi tabelke z "loadings"

Kod:

library(Hmisc)
latex(summary(rob.pca)[[2]] ,title="Robust PCA",
file="my.tex", append=T, first.hline.double=F,table.env=F)

No to sie pozalilem :)

Crunchy · Wysłany: 2010-02-27, 15:37

Cytat:

R is free software and comes with ABSOLUTELY NO WARRANTY.

Na pocieszenie powiem, że całkiem niedawno odkryłem, że w pewnym komercyjnym pakiecie program zwraca zamiast dystrybuanty, wartość punktową...

Kod:

latex((summary(rob.pca$sdev)^2) ,title="Robust PCA",
file="my.tex", append=T, first.hline.double=F,table.env=F)

Czy to działa? (bo nie sprawdzałem...)
Przy R helpa niestety trzeba szukać gdzie indziej...

bulva · Chorąży **Pomógł:** 5 razy Posty: 172 Skąd: Zgierz

Hehe, no tak :)
U mnie się ostatnio pokazuje jeszcze :

Cytat:

Type 'revo()' to visit www.revolution-computing.com for the latest
REvolution R news, 'forum()' for the community forum, or 'readme()'
for release notes.

Ale to chyba jakaś reklama jest bo Revolution jest komercyjną wersją R o ile sie nie myle?

Jak bym chciał żeby ta tabelka zwróciła mi "Proportion of Variance" i "Cumulative Proportion" i do spółki ze screeplotem pokazywała ile całkowitej wariancji łapią poszczególne PC... No nic niech stracę, przepisze sobie raz taką tabelkę, w wolnej chwili może opracuje jakiś kod który to będzie zwracał jak trzeba.

Crunchy · Wysłany: 2010-02-28, 18:36

Revolution ma kilka wersji, ale ta ogólnie dostępna jest zbudowana na R-2.7.X chyba, mi zbyt wiele nie dała, ale właśnie próbuję GotoBLAS 2 + optymalizacje pod swój procesor przez Cygwina, zanim się przesiądę na Archa...

bulva · Chorąży **Pomógł:** 5 razy Posty: 172 Skąd: Zgierz

Polecam, ale mnie nie pasowało domyślne KDE jako środowisko... Zresztą sam sobie go dopieścisz :)

[ Dodano: 2010-03-02, 20:54 ]
Dla potomnosci, jak to zrobic zeby sie nie narobic:

Kod:

Eigen <- function (pca.obj)
{
N=length(pca.obj$sdev)
Prop.var=rep(NA,N)
n=rep(NA,N)
for (i in 1:N)
{
Prop.var[i]=pca.obj$sdev[i]^2/sum(pca.obj$sdev^2)
}
Cumulative.var=cumsum(Prop.var)
eigen <- data.frame(rbind(Prop.var, Cumulative.var))
for (i in 1:N)
{
n[i] <- paste("Comp.",i,"")
}
names(eigen)[1:N]<-n
return(eigen)
}

latex(round(Eigen(rob.pca),4) ,title="Robust PCA",
file="my.tex", append=T, first.hline.double=F,table.env=F)

Crunchy · Wysłany: 2010-03-03, 22:35

A tak swoją drogą, czy nie wiesz jak ustawić flagi dla pakietów globalnie (budując R ze źródeł) zamiast dla każdego oddzielnie..?

bulva · Chorąży **Pomógł:** 5 razy Posty: 172 Skąd: Zgierz

Hmmm w Archu chyba się ustawiało domyślne flagi w /etc/makepkg.conf...

Crunchy · Wysłany: 2010-03-04, 18:41

Trauma Gentoo wróciła :-D

Ale raczej nie chodzi mi o to

Cytat:

The R system and package-specific compilation flags can be overridden or added to by setting the appropriate Make variables in the personal file HOME/.R/Makevars-R_PLATFORM (but HOME/.R/Makevars.win on Windows), or if that does not exist, HOME/.R/Makevars, where ‘R_PLATFORM’ is the platform for which R was built, as available in the platform component of the R variable R.version.

Flagi definiowane są w kilku różnych miejscach, a każdy pakiet ma swoje w Make'ach...

bulva · Chorąży **Pomógł:** 5 razy Posty: 172 Skąd: Zgierz

Tu jestem bezsilny. Ale jak się dowiesz chętnie poczytam :)

Crunchy · Wysłany: 2010-03-08, 18:12

/R-*/etc/Makeconf :mrgreen:

Optymalizacja pakietów mi dała średnio 7% krótszy czas obliczeń, samego R nawet nie 1%, wywalenie BLASa dodatkowo 13% 8-)

Czyli przy odrobinie chęci można przyśpieszyć niektóre rzeczy nawet o więcej niż 20%!

bulva · Chorąży **Pomógł:** 5 razy Posty: 172 Skąd: Zgierz

Jak ze wszystkim opłaca się być blisko bitów :) - następny system jaki postawie będzie z Ubuntu minimal + własny build kernela. Przyjemne z pożytecznym, repozytoria Ubuntu i kontrola nad KAŻDYM pakietem w systemie. Tylko co nie tak z BLASem?

Crunchy · Wysłany: 2010-03-09, 21:42

BLAS to nie jest najwydajniejsza biblioteka, po prostu... najprościej, podmień na ATLAS i się przekonasz, dla AMD można jeszcze zobaczyć ACML, ale z tego co wiem Goto BLAS daje najlepsze wyniki (z Texas Advanced Computing Center).