r - R の NA を最も近い値に置き換える

Question

na.locf()パッケージに似たものを探していますが、以前の非値をzoo常に使用する代わりに、最も近い非値を使用したいと思います。いくつかのサンプルデータ:NANA

dat <- c(1, 3, NA, NA, 5, 7)

(3 は繰り越される)NAに置き換えます。na.locf

library(zoo)
na.locf(dat)
# 1 3 3 3 5 7

に設定するとna.locf( 5 は後方に繰り越されます):fromLastTRUE

na.locf(dat, fromLast = TRUE)
# 1 3 5 5 5 7

しかし、最も近い非NA値が使用されることを望みます。私の例では、これは 3 を最初のNAに繰り越し、5 を 2 番目に繰り戻す必要があることを意味しNAます。

1 3 3 5 5 7

解決策をコード化しましたが、車輪の再発明ではないことを確認したかったのです。すでに何かが浮かんでいますか？

参考までに、私の現在のコードは次のとおりです。おそらく、他に何もないとしても、誰かがそれをより効率的にする方法を提案できます。これを改善するための明らかな方法が欠けているように感じます：

  na.pos <- which(is.na(dat))
  if (length(na.pos) == length(dat)) {
    return(dat)
  }
  non.na.pos <- setdiff(seq_along(dat), na.pos)
  nearest.non.na.pos <- sapply(na.pos, function(x) {
    return(which.min(abs(non.na.pos - x)))
  })
  dat[na.pos] <- dat[non.na.pos[nearest.non.na.pos]]

以下のsmciの質問に答えるには：

1. いいえ、どのエントリも NA にできます
2. すべてが NA の場合は、そのままにしておきます
3. いいえ。私の現在の解決策はデフォルトで左側の最も近い値ですが、問題ではありません
4. 通常、これらの行は数十万の要素であるため、理論的には上限は数十万になります。実際には、あちこちに数個しかなく、通常は 1 つです。

更新したがって、私たちはまったく別の方向に進んでいることがわかりましたが、これは依然として興味深い議論でした。皆さんありがとう！

Answer 1

これは非常に速いものです。元のデータfindIntervalのそれぞれについて考慮すべき2つの位置を見つけるために使用します。NA

f1 <- function(dat) {
  N <- length(dat)
  na.pos <- which(is.na(dat))
  if (length(na.pos) %in% c(0, N)) {
    return(dat)
  }
  non.na.pos <- which(!is.na(dat))
  intervals  <- findInterval(na.pos, non.na.pos,
                             all.inside = TRUE)
  left.pos   <- non.na.pos[pmax(1, intervals)]
  right.pos  <- non.na.pos[pmin(N, intervals+1)]
  left.dist  <- na.pos - left.pos
  right.dist <- right.pos - na.pos

  dat[na.pos] <- ifelse(left.dist <= right.dist,
                        dat[left.pos], dat[right.pos])
  return(dat)
}

そしてここで私はそれをテストします：

# sample data, suggested by @JeffAllen
dat <- as.integer(runif(50000, min=0, max=10))
dat[dat==0] <- NA

# computation times
system.time(r0 <- f0(dat))    # your function
# user  system elapsed 
# 5.52    0.00    5.52
system.time(r1 <- f1(dat))    # this function
# user  system elapsed 
# 0.01    0.00    0.03
identical(r0, r1)
# [1] TRUE

Answer 2

以下のコード。最初の質問は完全に明確に定義されていませんでした。私はこれらの説明を求めていました。

1. 少なくとも最初および/または最後のエントリが非 NA であることが保証されていますか? [いいえ]
2. 行のすべてのエントリが NA の場合はどうすればよいですか? 【そのまま放置】
3. タイがどのように分割されるか、つまり の真ん中の NA をどのように扱うか気にします1 3 NA NA NA 5 7か? [どうでもいい・左]
4. 連続する NA の最長連続スパンに上限 (S) はありますか? （Sが小さい場合は再帰的なソリューションを考えています。または、ifelseSが大きく、行と列の数が多い場合のデータフレームソリューション。）[最悪の場合、Sは病的に大きくなる可能性があるため、再帰を使用しないでください]

geoffjentry、あなたのソリューションでは、ボトルネックはシリアル計算nearest.non.na.posとシリアル割り当てになりますdat[na.pos] <- dat[non.na.pos[nearest.non.na.pos]] 長さ G の大きなギャップの場合、実際に計算する必要があるのは、最初の (G/2、切り上げ) アイテムが左から塗りつぶされることだけです。右から休みます。(使用して回答を投稿できifelseますが、似たようなものになります。) 基準は、ランタイム、big-O 効率、一時メモリ使用量、またはコードの読みやすさですか?

カップラ可能な微調整:

• N <- length(dat)一度だけ計算する必要があります
• 一般的な速度向上: if (length(na.pos) == 0)NA がないため、行をスキップします
• if (length(na.pos) == length(dat)-1)非 NA エントリが 1 つしかない (まれな) ケースで、行全体を入力します。

アウトライン ソリューション:

残念ながら、na.locf はデータフレーム全体では機能しません。行ごとに sapply を使用する必要があります。

na.fill_from_nn <- function(x) {
  row.na <- is.na(x)
  fillFromLeft <- na.locf(x, na.rm=FALSE) 
  fillFromRight <- na.locf(x, fromLast=TRUE, na.rm=FALSE)

  disagree <- rle(fillFromLeft!=fillFromRight)
  for (loc in (disagree)) { ...  resolve conflicts, row-wise }
}

sapply(dat, na.fill_from_nn)

または、あなたが言うように、連続した NA はまれであるためifelse、左から孤立した NA を埋めるために高速でダムを実行します。これはデータフレームごとに動作します=>一般的なケースを高速にします。次に、行単位の for ループで他のすべてのケースを処理します。(これは、NA の長いスパンの中間要素のタイブレークに影響しますが、気にしないと言います。)

Answer 3

明らかな単純な解決策は思いつきませんが、提案 (特にsmciの使用に関する提案rle) を見て、より効率的であると思われる複雑な関数を思いつきました。

これがコードです。以下で説明します。

# Your function
your.func = function(dat) {
  na.pos <- which(is.na(dat))
  if (length(na.pos) == length(dat)) {
    return(dat)
  }
  non.na.pos <- setdiff(seq_along(dat), na.pos)
  nearest.non.na.pos <- sapply(na.pos, function(x) which.min(abs(non.na.pos - x)))
  dat[na.pos] <- dat[non.na.pos[nearest.non.na.pos]]
  dat
}

# My function
my.func = function(dat) {
    nas=is.na(dat)
    if (!any(!nas)) return (dat)
    t=rle(nas)
    f=sapply(t$lengths[t$values],seq)
    a=unlist(f)
    b=unlist(lapply(f,rev))
    x=which(nas)
    l=length(dat)
    dat[nas]=ifelse(a>b,dat[ ifelse((x+b)>l,x-a,x+b) ],dat[ifelse((x-a)<1,x+b,x-a)])
    dat
}


# Test
n = 100000
test.vec = 1:n
set.seed(1)
test.vec[sample(test.vec,n/4)]=NA

system.time(t1<-my.func(test.vec))
system.time(t2<-your.func(test.vec)) # 10 times speed improvement on my machine

# Verify
any(t1!=t2)

私の機能はに依存していrleます。上記のコメントを読んでいますがrle、NA. 小さな例で説明するのが最も簡単です。

ベクトルから始める場合:

dat=c(1,2,3,4,NA,NA,NA,8,NA,10,11,12,NA,NA,NA,NA,NA,18)

次に、すべての NA の位置を取得します。

x=c(5,6,7,8,13,14,15,16,17)

次に、NA の「実行」ごとに、1 から実行の長さまでのシーケンスを作成します。

a=c(1,2,3,1,1,2,3,4,5)

次に、もう一度実行しますが、シーケンスを逆にします。

b=c(3,2,1,1,5,4,3,2,1)

これで、ベクトル a と b を比較できます。a<=b の場合は、振り返って xa の値を取得します。a>b の場合、先を見越して x+b の値を取得します。残りは、ベクトルの最後または最初にすべての NA または NA 実行がある場合のコーナー ケースを処理するだけです。

もっと簡単な解決策があるかもしれませんが、これで始められることを願っています。

Answer 4

これが私の突き刺しです。R で for ループを見るのは好きではありませんが、まばらな NA ベクトルの場合、実際にはより効率的であるように見えます (以下のパフォーマンス メトリック)。コードの要点は以下です。

  #get the index of all NA values
  nas <- which(is.na(dat))

  #get the Boolean map of which are NAs, used later to determine which values can be used as a replacement, and which are just filled-in NA values
  namask <- is.na(dat)

  #calculate the maximum size of a run of NAs
  length <- getLengthNAs(dat);

  #the furthest away an NA value could be is half of the length of the maximum NA run
  windowSize <- ceiling(length/2)

  #loop through all NAs
  for (thisIndex in nas){
    #extract the neighborhood of this NA
    neighborhood <- dat[(thisIndex-windowSize):(thisIndex+windowSize)]
    #any already-filled-in values which were NA can be replaced with NAs
    neighborhood[namask[(thisIndex-windowSize):(thisIndex+windowSize)]] <- NA

    #the center of this neighborhood
    center <- windowSize + 1

    #compute the difference within this neighborhood to find the nearest non-NA value
    delta <- center - which(!is.na(neighborhood))

    #find the closest replacement
    replacement <- delta[abs(delta) == min(abs(delta))]
    #in case length > 1, just pick the first
    replacement <- replacement[1]

    #replace with the nearest non-NA value.
    dat[thisIndex] <- dat[(thisIndex - (replacement))]
  }

あなたが提案したコードは気に入りましたが、マトリックス内のすべての NA 値と他のすべての非 NA インデックスの間のデルタを計算していることに気付きました。これが最大のパフォーマンス豚だったと思います。代わりに、各 NA の周囲の最小サイズの近傍またはウィンドウを抽出し、そのウィンドウ内で最も近い非 NA 値を見つけます。

そのため、パフォーマンスは、NA の数とウィンドウ サイズに比例して変化します。ここで、ウィンドウ サイズは、NA の最大実行の長さの半分 (の上限) です。NA の最大ランの長さを計算するには、次の関数を使用できます。

getLengthNAs <- function(dat){
  nas <- which(is.na(dat))
  spacing <- diff(nas)
  length <- 1;
  while (any(spacing == 1)){        
    length <- length + 1;
    spacing <- diff(which(spacing == 1))
  }
    length
}

性能比較

#create a test vector with 10% NAs and length 50,000.
dat <- as.integer(runif(50000, min=0, max=10))
dat[dat==0] <- NA

#the a() function is the code posted in the question
a <- function(dat){
  na.pos <- which(is.na(dat))
    if (length(na.pos) == length(dat)) {
        return(dat)
    }
    non.na.pos <- setdiff(seq_along(dat), na.pos)
    nearest.non.na.pos <- sapply(na.pos, function(x) {
        return(which.min(abs(non.na.pos - x)))
    })
    dat[na.pos] <- dat[non.na.pos[nearest.non.na.pos]]
    dat
}

#my code
b <- function(dat){
    #the same code posted above, but with some additional helper code to sanitize the input
    if(is.null(dat)){
      return(NULL);
    }

    if (all(is.na(dat))){
      stop("Can't impute NAs if there are no non-NA values.")
    }

    if (!any(is.na(dat))){
      return(dat);
    }

    #starts with an NA (or multiple), handle these
    if (is.na(dat[1])){
      firstNonNA <- which(!is.na(dat))[1]
      dat[1:(firstNonNA-1)] <- dat[firstNonNA]
    }

    #ends with an NA (or multiple), handle these
    if (is.na(dat[length(dat)])){
      lastNonNA <- which(!is.na(dat))
      lastNonNA <- lastNonNA[length(lastNonNA)]
      dat[(lastNonNA+1):length(dat)] <- dat[lastNonNA]
    }

    #get the index of all NA values
    nas <- which(is.na(dat))

    #get the Boolean map of which are NAs, used later to determine which values can be used as a replacement, and which are just filled-in NA values
    namask <- is.na(dat)

    #calculate the maximum size of a run of NAs
    length <- getLengthNAs(dat);

    #the furthest away an NA value could be is half of the length of the maximum NA run
    #if there's a run at the beginning or end, then the nearest non-NA value could possibly be `length` away, so we need to keep the window large for that case.
    windowSize <- ceiling(length/2)

    #loop through all NAs
    for (thisIndex in nas){
      #extract the neighborhood of this NA
      neighborhood <- dat[(thisIndex-windowSize):(thisIndex+windowSize)]
      #any already-filled-in values which were NA can be replaced with NAs
      neighborhood[namask[(thisIndex-windowSize):(thisIndex+windowSize)]] <- NA

      #the center of this neighborhood
      center <- windowSize + 1

      #compute the difference within this neighborhood to find the nearest non-NA value
      delta <- center - which(!is.na(neighborhood))

      #find the closest replacement
      replacement <- delta[abs(delta) == min(abs(delta))]
      #in case length > 1, just pick the first
      replacement <- replacement[1]

      #replace with the nearest non-NA value.
      dat[thisIndex] <- dat[(thisIndex - (replacement))]
    }
    dat
}

#nograpes' answer on this question
c <- function(dat){
  nas=is.na(dat)
  if (!any(!nas)) return (dat)
  t=rle(nas)
  f=sapply(t$lengths[t$values],seq)
  a=unlist(f)
  b=unlist(lapply(f,rev))
  x=which(nas)
  l=length(dat)
  dat[nas]=ifelse(a>b,dat[ ifelse((x+b)>l,x-a,x+b) ],dat[ifelse((x-a)<1,x+b,x-a)])
  dat
}

#run 10 times each to get average performance.
sum <- 0; for (i in 1:10){ sum <- sum + system.time(a(dat))["elapsed"];}; cat ("A: ", sum/10)
A:  5.059
sum <- 0; for (i in 1:10){ sum <- sum + system.time(b(dat))["elapsed"];}; cat ("B: ", sum/10)
B:  0.126
sum <- 0; for (i in 1:10){ sum <- sum + system.time(c(dat))["elapsed"];}; cat ("C: ", sum/10)
C:  0.287

したがって、このコードのように（少なくともこれらの条件下では）、質問に投稿された元のコードから約 40 倍のスピードアップを提供し、以下の @nograpes の回答よりも 2.2 倍のスピードアップを提供します（rleソリューションは確かに高速になると思いますいくつかの状況 -- NA が豊富なベクトルを含む)。