2012-07-18 75 views
7

我試圖實現Chebyshev濾波器來平滑時間序列,但不幸的是,在數據序列中有NAs。R濾波器()處理NA

例如,

t <- seq(0, 1, len = 100)      
x <- c(sin(2*pi*t*2.3) + 0.25*rnorm(length(t)),NA, cos(2*pi*t*2.3) + 0.25*rnorm(length(t))) 

我用的切比雪夫濾波器:cf1 = cheby1(5, 3, 1/44, type = "low")

我試圖過濾器系列排除NAS上的時間,但不會弄亂訂單/位置。所以,我已經試過na.rm=T,但似乎沒有這樣的說法。 Then

z <- filter(cf1, x) # apply filter 

謝謝你們。

回答

1

您可以使用compelete.cases函數預先刪除NAs。你也可以考慮輸入缺失的數據。查看mtsdi或Amelia II軟件包。

編輯:

下面是與RCPP的解決方案。這可能是有幫助的是速度是非常重要的:

require(inline) 
require(Rcpp) 
t <- seq(0, 1, len = 100) 
set.seed(7337) 
x <- c(sin(2*pi*t*2.3) + 0.25*rnorm(length(t)),NA, cos(2*pi*t*2.3) + 0.25*rnorm(length(t))) 
NAs <- x 
x2 <- x[!is.na(x)] 
#do something to x2 
src <- ' 
Rcpp::NumericVector vecX(vx); 
Rcpp::NumericVector vecNA(vNA); 
int j = 0; //counter for vx 
for (int i=0;i<vecNA.size();i++) { 
    if (!(R_IsNA(vecNA[i]))) { 
    //replace and update j 
    vecNA[i] = vecX[j]; 
    j++; 
    } 
} 
return Rcpp::wrap(vecNA); 
' 
fun <- cxxfunction(signature(vx="numeric", 
          vNA="numeric"), 
        src,plugin="Rcpp") 
if (identical(x,fun(x2,NAs))) 
    print("worked") 
# [1] "worked" 
+0

我只是想知道complete.case和na.omit是否一樣。另外,由於我使用的是觀測到的SST時間序列,我不確定輸入缺失值是否是一個好主意。 – 2012-07-18 13:15:40

+0

希望這個更新解決了這個問題。 – chandler 2012-07-19 11:37:15

2

嘗試使用x <- x[!is.na(x)]刪除NA,然後運行過濾器。

+0

對不起,但如果我使用na.omit,它將大量增加訂單,我只想在過濾器剩餘NA之後的NA值,但所有其他非NA值都可以通過。 – 2012-07-18 12:52:07

+0

對不起,我不確定你在問什麼。你想保持你的價值觀在同一個位置,並且在你的數據中有空格嗎? – 2012-07-18 12:54:50

+0

將tseries包中的na.remove()或na.remove.ts()做你想要的嗎? – 2012-07-18 13:01:15

1

我不知道,如果ts對象可以有缺失值,但如果你只是想重新插入NA值,你可以使用?insertR.utils。可能有更好的方法來做到這一點。

install.packages(c('R.utils', 'signal')) 
require(R.utils) 
require(signal) 
t <- seq(0, 1, len = 100)      
set.seed(7337) 
x <- c(sin(2*pi*t*2.3) + 0.25*rnorm(length(t)), NA, NA, cos(2*pi*t*2.3) + 0.25*rnorm(length(t)), NA) 
cf1 = cheby1(5, 3, 1/44, type = "low") 
xex <- na.omit(x) 
z <- filter(cf1, xex) # apply 
z <- as.numeric(z) 
for (m in attributes(xex)$na.action) { 
    z <- insert(z, ats = m, values = NA) 
} 
all.equal(is.na(z), is.na(x)) 
?insert 
0

這裏是一個函數,你可以用它來過濾一個信號,其中含有NAs。 NA被忽略而不是被零代替。

然後,您可以指定NA在濾波後信號的任意點可能採取的最大權重百分比。如果在特定點有太多的NAs(並且實際數據太少),則濾波後的信號本身將被設置爲NA。

# This function applies a filter to a time series with potentially missing data 

filter_with_NA <- function(x, 
          window_length=12,       # will be applied centrally 
          myfilter=rep(1/window_length,window_length), # a boxcar filter by default 
          max_percentage_NA=25)      # which percentage of weight created by NA should not be exceeded 
{ 
    # make the signal longer at both sides 
    signal <- c(rep(NA,window_length),x,rep(NA,window_length)) 
    # see where data are present and not NA 
    present <- is.finite(signal) 

    # replace the NA values by zero 
    signal[!is.finite(signal)] <- 0 
    # apply the filter 
    filtered_signal <- as.numeric(filter(signal,myfilter, sides=2)) 

    # find out which percentage of the filtered signal was created by non-NA values 
    # this is easy because the filter is linear 
    original_weight <- as.numeric(filter(present,myfilter, sides=2)) 
    # where this is lower than one, the signal is now artificially smaller 
    # because we added zeros - compensate that 
    filtered_signal <- filtered_signal/original_weight 
    # but where there are too few values present, discard the signal 
    filtered_signal[100*(1-original_weight) > max_percentage_NA] <- NA 

    # cut away the padding to left and right which we previously inserted 
    filtered_signal <- filtered_signal[((window_length+1):(window_length+length(x)))] 
    return(filtered_signal) 
}