【问题标题】:R filter() dealing with NAsR filter() 处理 NA
【发布时间】:2012-07-17 11:41:23
【问题描述】:

我正在尝试实现 Chebyshev 滤波器来平滑时间序列,但不幸的是,数据序列中有 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")

我正在尝试过滤排除 NA 的时间序列,但不会弄乱订单/仓位。所以,我已经尝试过na.rm=T,但似乎没有这样的论点。 那么

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

谢谢你们。

【问题讨论】:

    标签: r filter signals time-series na


    【解决方案1】:

    尝试使用x &lt;- x[!is.na(x)] 删除 NA,然后运行过滤器。

    【讨论】:

    • 对不起,但是如果我使用 na.omit,它会堆积订单,我只想要过滤器剩余 NA 后的 NA 值,但所有其他非 NA 值都可以通过。
    • 抱歉,我不确定你在问什么。您想将您的值保持在同一位置,并在您的数据中有空格吗?
    • tseries 包中的 na.remove() 或 na.remove.ts() 会做你想做的事吗?
    • 非常感谢,但我想知道如果我使用 x
    • 你不能使用 z
    【解决方案2】:

    您可以使用 compelete.cases 函数预先删除 NA。您也可以考虑估算缺失的数据。查看 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"
    

    【讨论】:

    • 我只是想知道 complete.case 是否与 na.omit 相同。另外,由于我使用的是观察到的 SST 时间序列,我不确定输入缺失值是否是个好主意。
    【解决方案3】:

    我不知道ts 对象是否可以有缺失值,但如果您只想重新插入NA 值,您可以使用R.utils 中的?insert。可能有更好的方法来做到这一点。

    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
    

    【讨论】:

      【解决方案4】:

      这是一个函数,您可以使用它来过滤带有 NA 的信号。 NA 被忽略而不是被零替换。

      然后,您可以指定 NA 在滤波信号的任何点上可以采用的最大权重百分比。如果某个特定点的 NA 太多(而实际数据太少),则滤波后的信号本身将设置为 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)
      }
      

      【讨论】:

        猜你喜欢
        • 1970-01-01
        • 2023-01-19
        • 2017-06-23
        • 2017-04-17
        • 2020-01-25
        • 2012-01-17
        • 1970-01-01
        • 1970-01-01
        • 2011-10-29
        相关资源
        最近更新 更多