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.
źródło
Odpowiedzi:
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:
Tak więc prawdziwe pole danych
Xt
składa się z 9 sygnałów i dodałem do niego trochę szumu, aby utworzyć obserwowane poleXp
, które zostanie wykorzystane w poniższych przykładach. EOF określa się jako takie:EOF
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
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
Zatem 5 wiodących EOF leży powyżej tej linii. Próbowałem tego przykładu, gdy
Xp
nie 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
wq
pakiecie znajduje się implementacja , którą pokazuję poniżej.Reguła N
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 danychXp
, aby każda kolumna była losowo przetasowana. W ten sposób zapewniamy, że całkowita wariancja pola losowego jest taka sama jakXp
. Wielokrotnie próbkując, jesteśmy w stanie obliczyć podstawowy błąd rozkładu.Monte Carlo z polem losowym (tj. Porównanie modelu zerowego)
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
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.
źródło