Jak określić znaczące główne komponenty za pomocą ładowania początkowego lub podejścia Monte Carlo?

Interesuje mnie określenie liczby znaczących wzorców pochodzących z analizy głównych składników (PCA) lub analizy empirycznej funkcji ortogonalnej (EOF). Jestem szczególnie zainteresowany zastosowaniem tej metody do danych klimatycznych. Pole danych jest macierzą MxN, gdzie M jest wymiarem czasowym (np. Dni), a N jest wymiarem przestrzennym (np. Lokalizacje lon / lat). Czytałem o możliwej metodzie bootstrap w celu ustalenia znaczących komputerów, ale nie byłem w stanie znaleźć bardziej szczegółowego opisu. Do tej pory stosowałem Regułę kciuka Northa (North i in ., 1982), aby ustalić tę granicę, ale zastanawiałem się, czy dostępna jest bardziej niezawodna metoda.

Jako przykład:

###Generate data
x <- -10:10
y <- -10:10
grd <- expand.grid(x=x, y=y)

#3 spatial patterns
sp1 <- grd$x^3+grd$y^2
tmp1 <- matrix(sp1, length(x), length(y))
image(x,y,tmp1)

sp2 <- grd$x^2+grd$y^2
tmp2 <- matrix(sp2, length(x), length(y))
image(x,y,tmp2)

sp3 <- 10*grd$y
tmp3 <- matrix(sp3, length(x), length(y))
image(x,y,tmp3)


#3 respective temporal patterns
T <- 1:1000

tp1 <- scale(sin(seq(0,5*pi,,length(T))))
plot(tp1, t="l")

tp2 <- scale(sin(seq(0,3*pi,,length(T))) + cos(seq(1,6*pi,,length(T))))
plot(tp2, t="l")

tp3 <- scale(sin(seq(0,pi,,length(T))) - 0.2*cos(seq(1,10*pi,,length(T))))
plot(tp3, t="l")


#make data field - time series for each spatial grid (spatial pattern multiplied by temporal pattern plus error)
set.seed(1)
F <- as.matrix(tp1) %*% t(as.matrix(sp1)) + 
as.matrix(tp2) %*% t(as.matrix(sp2)) + 
as.matrix(tp3) %*% t(as.matrix(sp3)) +
matrix(rnorm(length(T)*dim(grd)[1], mean=0, sd=200), nrow=length(T), ncol=dim(grd)[1]) # error term

dim(F)
image(F)


###Empirical Orthogonal Function (EOF) Analysis 
#scale field
Fsc <- scale(F, center=TRUE, scale=FALSE)

#make covariance matrix
C <- cov(Fsc)
image(C)

#Eigen decomposition
E <- eigen(C)

#EOFs (U) and associated Lambda (L) 
U <- E$vectors
L <- E$values

#projection of data onto EOFs (U) to derive principle components (A)
A <- Fsc %*% U

dim(U)
dim(A)

#plot of top 10 Lambda
plot(L[1:10], log="y")

#plot of explained variance (explvar, %) by each EOF
explvar <- L/sum(L) * 100
plot(explvar[1:20], log="y")


#plot original patterns versus those identified by EOF
layout(matrix(1:12, nrow=4, ncol=3, byrow=TRUE), widths=c(1,1,1), heights=c(1,0.5,1,0.5))
layout.show(12)

par(mar=c(4,4,3,1))
image(tmp1, main="pattern 1")
image(tmp2, main="pattern 2")
image(tmp3, main="pattern 3")

par(mar=c(4,4,0,1)) 
plot(T, tp1, t="l", xlab="", ylab="")
plot(T, tp2, t="l", xlab="", ylab="")
plot(T, tp3, t="l", xlab="", ylab="")

par(mar=c(4,4,3,1))
image(matrix(U[,1], length(x), length(y)), main="eof 1") 
image(matrix(U[,2], length(x), length(y)), main="eof 2")
image(matrix(U[,3], length(x), length(y)), main="eof 3")

par(mar=c(4,4,0,1)) 
plot(T, A[,1], t="l", xlab="", ylab="")
plot(T, A[,2], t="l", xlab="", ylab="")
plot(T, A[,3], t="l", xlab="", ylab="")

A oto metoda, której użyłem do określenia znaczenia komputera. Zasadniczo reguła jest taka, że ​​różnica między sąsiednimi Lambdas musi być większa niż związany z nimi błąd.

###Determine significant EOFs

#North's Rule of Thumb
Lambda_err <- sqrt(2/dim(F)[2])*L
upper.lim <- L+Lambda_err
lower.lim <- L-Lambda_err
NORTHok=0*L
for(i in seq(L)){
    Lambdas <- L
    Lambdas[i] <- NaN
    nearest <- which.min(abs(L[i]-Lambdas))
    if(nearest > i){
        if(lower.lim[i] > upper.lim[nearest]) NORTHok[i] <- 1
    }
    if(nearest < i){
        if(upper.lim[i] < lower.lim[nearest]) NORTHok[i] <- 1
    }
}
n_sig <- min(which(NORTHok==0))-1

plot(L[1:10],log="y", ylab="Lambda (dots) and error (vertical lines)", xlab="EOF")
segments(x0=seq(L), y0=L-Lambda_err, x1=seq(L), y1=L+Lambda_err)
abline(v=n_sig+0.5, col=2, lty=2)
text(x=n_sig, y=mean(L[1:10]), labels="North's Rule of Thumb", srt=90, col=2)

Uważam, że sekcja rozdziału autorstwa Björnssona i Venegasa ( 1997 ) na temat testów istotności jest pomocna - odnoszą się one do trzech kategorii testów, z których prawdopodobnie mam nadzieję, że użyję dominującej wariancji . Odnoszą się do rodzaju podejścia Monte Carlo polegającego na tasowaniu wymiaru czasu i ponownym obliczaniu Lambdas na podstawie wielu permutacji. von Storch i Zweiers (1999) również odnoszą się do testu, który porównuje widmo Lambda z referencyjnym widmem „szumowym”. W obu przypadkach nie jestem pewien, jak można to zrobić, a także w jaki sposób wykonuje się test istotności, biorąc pod uwagę przedziały ufności określone przez permutacje.

Dzięki za pomoc.

Odniesienia: Björnsson, H. and Venegas, SA (1997). „Podręcznik analiz EOF i SVD danych klimatycznych”, McGill University, Raport CCGCR nr 97-1, Montréal, Québec, 52 str. http://andvari.vedur.is/%7Efolk/halldor/PICKUP/eof.pdf

GR North, TL Bell, RF Cahalan i FJ Moeng. (1982). Błędy próbkowania w estymacji empirycznych funkcji ortogonalnych. Pon. Wea. Rev. 110: 699–706.

von Storch, H, Zwiers, FW (1999). Analiza statystyczna w badaniach klimatu. Cambridge University Press.

Spróbuję tu nieco przyspieszyć dialog, chociaż to moje pytanie. Minęło 6 miesięcy, odkąd o to zapytałem i niestety nie otrzymano pełnych odpowiedzi. Spróbuję podsumować to, co do tej pory zebrałem, i zobaczę, czy ktokolwiek może rozwinąć pozostałe kwestie. Proszę wybaczyć długą odpowiedź, ale nie widzę innego wyjścia ...

Najpierw przedstawię kilka podejść przy użyciu możliwie lepszego syntetycznego zestawu danych. Pochodzi z artykułu Beckersa i Rixona ( 2003 ) ilustrującego użycie algorytmu do przeprowadzania EOF na nieszczelnych danych. Algorytm odtworzyłem w języku R, jeśli ktoś jest zainteresowany ( link ).

Zestaw danych syntetycznych:

#color palette
pal <- colorRampPalette(c("blue", "cyan", "yellow", "red"))

#Generate data
m=50
n=100
frac.gaps <- 0.5 # the fraction of data with NaNs
N.S.ratio <- 0.25 # the Noise to Signal ratio for adding noise to data

x <- (seq(m)*2*pi)/m
t <- (seq(n)*2*pi)/n


#True field
Xt <- 
 outer(sin(x), sin(t)) + 
 outer(sin(2.1*x), sin(2.1*t)) + 
 outer(sin(3.1*x), sin(3.1*t)) +
 outer(tanh(x), cos(t)) + 
 outer(tanh(2*x), cos(2.1*t)) + 
 outer(tanh(4*x), cos(0.1*t)) + 
 outer(tanh(2.4*x), cos(1.1*t)) + 
 tanh(outer(x, t, FUN="+")) + 
 tanh(outer(x, 2*t, FUN="+"))

Xt <- t(Xt)
image(Xt, col=pal(100))

#Noise field
set.seed(1)
RAND <- matrix(runif(length(Xt), min=-1, max=1), nrow=nrow(Xt), ncol=ncol(Xt))
R <- RAND * N.S.ratio * Xt

#True field + Noise field
Xp <- Xt + R
image(Xp, col=pal(100))

Tak więc prawdziwe pole danych Xtskłada się z 9 sygnałów i dodałem do niego trochę szumu, aby utworzyć obserwowane pole Xp, które zostanie wykorzystane w poniższych przykładach. EOF określa się jako takie:

EOF

#make covariance matrix
C <- t(Xp) %*% Xp #cov(Xp)
image(C)

#Eigen decomposition
E <- svd(C)

#EOFs (U) and associated Lambda (L) 
U <- E$u
L <- E$d

#projection of data onto EOFs (U) to derive principle components (A)
A <- Xp %*% U

Idąc za przykładem, którego użyłem w moim oryginalnym przykładzie, określę „znaczące” EOF na podstawie praktycznej reguły Northa.

Reguła kciuka Northa

Lambda_err <- sqrt(2/dim(Xp)[2])*L
upper.lim <- L+Lambda_err
lower.lim <- L-Lambda_err
NORTHok=0*L
for(i in seq(L)){
    Lambdas <- L
    Lambdas[i] <- NaN
    nearest <- which.min(abs(L[i]-Lambdas))
    if(nearest > i){
        if(lower.lim[i] > upper.lim[nearest]) NORTHok[i] <- 1
    }
    if(nearest < i){
        if(upper.lim[i] < lower.lim[nearest]) NORTHok[i] <- 1
    }
}
n_sig <- min(which(NORTHok==0))-1
n_sig

plot(L[1:20],log="y", ylab="Lambda (dots) and error (vertical lines)", xlab="EOF")
segments(x0=seq(L), y0=L-Lambda_err, x1=seq(L), y1=L+Lambda_err)
abline(v=n_sig+0.5, col=2, lty=2)
text(x=n_sig, y=mean(L[1:10]), labels="North's Rule of Thumb", srt=90, col=2)

Ponieważ wartości Lambda 2: 4 są bardzo blisko siebie pod względem amplitudy, są one uważane za nieistotne zgodnie z zasadą kciuka - tj. Ich odpowiednie wzorce EOF mogą nakładać się i mieszać, biorąc pod uwagę ich podobne amplitudy. Jest to niefortunne, biorąc pod uwagę, że wiemy, że 9 sygnałów faktycznie istnieje w terenie.

Bardziej subiektywnym podejściem jest przeglądanie przekształconych logarytmicznie wartości Lambda („wykres Scree”), a następnie dopasowanie regresji do wartości końcowych. Następnie można określić wizualnie, na jakim poziomie wartości lambda leżą powyżej tej linii:

Działka piaskowa

ntrail <- 35
tail(L, ntrail)
fit <- lm(log(tail(L, ntrail)) ~ seq(length(L)-ntrail+1, length(L)))
plot(log(L))
abline(fit, col=2)

Zatem 5 wiodących EOF leży powyżej tej linii. Próbowałem tego przykładu, gdy Xpnie dodano żadnego dodatkowego szumu, a wyniki ujawniają wszystkie 9 oryginalnych sygnałów. Tak więc nieistotność EOF 6: 9 wynika z faktu, że ich amplituda jest mniejsza niż hałas w polu.

Bardziej obiektywną metodą są kryteria „Reguły N” autorstwa Overland i Preisendorfer (1982). W wqpakiecie znajduje się implementacja , którą pokazuję poniżej.

Reguła N

library(wq)
eofNum(Xp, distr = "normal", reps = 99)

RN <- ruleN(nrow(Xp), ncol(Xp), type = "normal", reps = 99)
RN
eigs <- svd(cov(Xp))$d
plot(eigs, log="y")
lines(RN, col=2, lty=2)

Zasada N zidentyfikowała 4 znaczące EOF. Osobiście muszę lepiej zrozumieć tę metodę; Dlaczego można zmierzyć poziom błędu na podstawie losowego pola, które nie używa tego samego rozkładu, co w Xp? Jedną z odmian tej metody byłoby ponowne próbkowanie danych Xp, aby każda kolumna była losowo przetasowana. W ten sposób zapewniamy, że całkowita wariancja pola losowego jest taka sama jak Xp. Wielokrotnie próbkując, jesteśmy w stanie obliczyć podstawowy błąd rozkładu.

Monte Carlo z polem losowym (tj. Porównanie modelu zerowego)

iter <- 499
LAMBDA <- matrix(NaN, ncol=iter, nrow=dim(Xp)[2])

set.seed(1)
for(i in seq(iter)){
    #i=1

    #random reorganize dimensions of scaled field
    Xp.tmp <- NaN*Xp
    for(j in seq(dim(Xp.tmp)[2])){
        #j=1
        Xp.tmp[,j] <- Xp[,j][sample(nrow(Xp))]
    }

    #make covariance matrix
    C.tmp <- t(Xp.tmp) %*% Xp.tmp #cov(Xp.tmp)

    #SVD decomposition
    E.tmp <- svd(C.tmp)

    #record Lambda (L) 
    LAMBDA[,i] <- E.tmp$d

    print(paste(round(i/iter*100), "%", " completed", sep=""))
}

boxplot(t(LAMBDA), log="y", col=8, border=2, outpch="")
points(L)

Ponownie 4 EOF są powyżej rozkładów dla pól losowych. Martwię się tym podejściem i zasadą N, że nie odnoszą się one do przedziałów ufności wartości Lambda; np. wysoka pierwsza wartość Lambda automatycznie spowoduje mniejszą wariancję, którą należy wyjaśnić przez końcowe. Zatem Lambda obliczona na podstawie losowych pól zawsze będzie miała mniej strome nachylenie i może spowodować wybranie zbyt małej liczby EOF. [UWAGA: eofNum()Funkcja zakłada, że ​​EOF są obliczane z macierzy korelacji. Ta liczba może być inna, jeśli np. Używasz macierzy kowariancji (dane wyśrodkowane, ale nie skalowane).]

Wreszcie @ tomasz74 wspomniał o pracy Babamoradi i in. (2013), na który krótko spojrzałem. Jest to bardzo interesujące, ale wydaje się, że bardziej koncentruje się na obliczaniu CI ładunków i współczynników EOF, niż Lambda. Niemniej jednak uważam, że można przyjąć ocenę błędu Lambda przy użyciu tej samej metodologii. Ponowne próbkowanie bootstrap odbywa się w polu danych poprzez ponowne próbkowanie wierszy, aż do utworzenia nowego pola. Ten sam wiersz może być ponownie próbkowany więcej niż raz, co jest podejściem nieparametrycznym i nie trzeba przyjmować założeń dotyczących dystrybucji danych.

Bootstrap wartości Lambda

B <- 40 * nrow(Xp)
LAMBDA <- matrix(NaN, nrow=length(L), ncol=B)
for(b in seq(B)){
    samp.b <- NaN*seq(nrow(Xp))
    for(i in seq(nrow(Xp))){
        samp.b[i] <- sample(nrow(Xp), 1)
    }
    Xp.b  <- Xp[samp.b,]
    C.b  <- t(Xp.b) %*% Xp.b 
    E.b  <- svd(C.b)
    LAMBDA[,b] <- E.b$d
    print(paste(round(b/B*100), "%", " completed", sep=""))
}
boxplot(t(LAMBDA), log="y", col=8, outpch="", ylab="Lambda [log-scale]")
points(L, col=4)
legend("topright", legend=c("Original"), pch=1, col=4)

Chociaż może to być bardziej solidna niż ogólna zasada Northa do obliczania błędu wartości Lambda, uważam teraz, że kwestia znaczenia EOF sprowadza się do różnych opinii na temat tego, co to znaczy. W przypadku metody kciuka i metody ładowania początkowego północy wydaje się, że znaczenie jest bardziej oparte na tym, czy teere nakładają się na siebie wartości Lambda. Jeśli tak, to te EOF mogą być mieszane w swoich sygnałach i nie reprezentować „prawdziwych” wzorców. Z drugiej strony, te dwa EOF mogą opisywać znaczną wariancję (w porównaniu z rozkładem pola losowego - np. Reguła N). Więc jeśli ktoś jest zainteresowany odfiltrowaniem szumu (tj. Poprzez obcięcie EOF), wystarczy reguła N. Jeśli ktoś jest zainteresowany określeniem rzeczywistych wzorców w zbiorze danych, bardziej rygorystyczne kryteria nakładania się lambda mogą być bardziej niezawodne.

Ponownie nie jestem ekspertem w tych kwestiach, więc nadal mam nadzieję, że ktoś z większym doświadczeniem może dodać do tego wyjaśnienia.

Bibliografia:

Beckers, Jean-Marie i M. Rixen. „Obliczenia EOF i wypełnianie danych z niekompletnych zestawów danych oceanograficznych”. Journal of Atmospheric and Oceanic Technology 20.12 (2003): 1839-1856.

Overland, J. i R. Preisendorfer, Test istotności głównych składników zastosowany w klimatologii cyklonów, pon. Wea. Rev., 110, 1-4, 1982.