R

在估算數據中使用鄰居信息或查找非數據(在 R 中)

  • May 28, 2014

我有數據集,假設最近的鄰居是最好的預測器。只是可視化的雙向漸變的完美示例-

在此處輸入圖像描述

假設我們有一些缺失值的情況,我們可以根據鄰居和趨勢輕鬆預測。

在此處輸入圖像描述

R中的相應數據矩陣(鍛煉的虛擬示例):

miss.mat <- matrix (c(5:11, 6:10, NA,12, 7:13, 8:14, 9:12, NA, 14:15, 10:16),ncol=7, byrow = TRUE)
miss.mat 
   [,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,]    5    6    7    8    9   10   11
[2,]    6    7    8    9   10   NA   12
[3,]    7    8    9   10   11   12   13
[4,]    8    9   10   11   12   13   14
[5,]    9   10   11   12   NA   14   15
[6,]   10   11   12   13   14   15   16

注: (1)缺失值的性質假定為隨機的,它可以在任何地方發生。

(2) 所有數據點均來自單個變量,但假設它們的值受neighbors與其相鄰的行和列的影響。所以矩陣中的位置很重要,可以被視為其他變量。

我希望在某些情況下我可以預測一些偏離值(可能是錯誤)並糾正偏差(例如,讓我們在虛擬數據中生成這樣的錯誤):

> mat2 <- matrix (c(4:10, 5, 16, 7, 11, 9:11, 6:12, 7:13, 8:14, 9:13, 4,15, 10:11, 2, 13:16),ncol=7, byrow = TRUE)
> mat2

   [,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,]    4    5    6    7    8    9   10
[2,]    5   16    7   11    9   10   11
[3,]    6    7    8    9   10   11   12
[4,]    7    8    9   10   11   12   13
[5,]    8    9   10   11   12   13   14
[6,]    9   10   11   12   13    4   15
[7,]   10   11    2   13   14   15   16

上面的例子只是說明性的(可以直觀地回答),但真實的例子可能更令人困惑。我正在尋找是否有可靠的方法來進行此類分析。我認為這應該是可能的。執行此類分析的合適方法是什麼?任何 R 程序/包建議來做這種類型的分析?

在此處輸入圖像描述

該問題要求以穩健的方式使用最近鄰**來識別和糾正局部異常值。為什麼不這樣做呢?

該過程是計算魯棒的局部平滑,評估殘差,並將任何太大的歸零。這直接滿足了所有要求,並且足夠靈活以適應不同的應用程序,因為可以改變本地鄰域的大小和識別異常值的閾值。

(為什麼靈活性如此重要?因為任何此類程序都有很好的機會將某些局部行為識別為“異常”。因此,所有此類程序都可以被認為是更平滑的。它們將消除一些細節以及明顯的異常值。分析師需要對保留細節和未能檢測到局部異常值之間的權衡進行一些控制。)

此過程的另一個優點是它不需要值的矩形矩陣。事實上,它甚至可以通過使用適合此類數據的局部平滑器來應用於不規則數據。

R以及大多數功能齊全的統計數據包,內置了幾個強大的本地平滑器,例如loess. 以下示例是使用它處理的。矩陣有行和列——幾乎條目。它表示一個複雜的函數,具有多個局部極值以及一整行不可微分的點(“摺痕”)。略高於在這些點中——被認為是“異常”的比例非常高——被添加了高斯誤差,其標準偏差僅為原始數據的標準差。因此,這個合成數據集呈現了現實數據的許多具有挑戰性的特徵。

數據

請注意(根據R慣例)矩陣行被繪製為垂直條。除殘差外,所有圖像都進行了山體陰影處理,以幫助顯示其值的微小變化。沒有這個,幾乎所有的局部異常值都將不可見!

通過比較“估算”(修復)和“真實”(原始未污染)圖像,很明顯,去除異常值已經平滑了一些但不是全部的摺痕(從向下; 它在“殘差”圖中顯示為淺青色斜條紋)。

“殘差”圖中的斑點顯示了明顯孤立的局部異常值。此圖還顯示可歸因於基礎數據的其他結構(例如斜條紋)。人們可以通過使用數據的空間模型(通過地統計方法)來改進這一過程,但在這裡描述和說明它會讓我們走得太遠。

順便說一句,此代碼報告僅發現的引入的異常值。這不是程序的失敗。 因為異常值是正態分佈的,所以大約有一半接近於零——與範圍超過——他們在表面上沒有發現明顯的變化。

#
# Create data.
#
set.seed(17)
rows <- 2:80; cols <- 2:50
y <- outer(rows, cols, 
          function(x,y) 100 * exp((abs(x-y)/50)^(0.9)) * sin(x/10) * cos(y/20))
y.real <- y
#
# Contaminate with iid noise.
#
n.out <- 200
cat(round(100 * n.out / (length(rows)*length(cols)), 2), "% errors\n", sep="")
i.out <- sample.int(length(rows)*length(cols), n.out)
y[i.out] <- y[i.out] + rnorm(n.out, sd=0.05 * sd(y))
#
# Process the data into a data frame for loess.
#
d <- expand.grid(i=1:length(rows), j=1:length(cols))
d$y <- as.vector(y)
#
# Compute the robust local smooth.
# (Adjusting `span` changes the neighborhood size.)
#
fit <- with(d, loess(y ~ i + j, span=min(1/2, 125/(length(rows)*length(cols)))))
#
# Display what happened.
#
require(raster)
show <- function(y, nrows, ncols, hillshade=TRUE, ...) {
 x <- raster(y, xmn=0, xmx=ncols, ymn=0, ymx=nrows)
 crs(x) <- "+proj=lcc +ellps=WGS84"
 if (hillshade) {
   slope <- terrain(x, opt='slope')
   aspect <- terrain(x, opt='aspect')
   hill <- hillShade(slope, aspect, 10, 60)
   plot(hill, col=grey(0:100/100), legend=FALSE, ...)
   alpha <- 0.5; add <- TRUE
 } else {
   alpha <- 1; add <- FALSE
 }
 plot(x, col=rainbow(127, alpha=alpha), add=add, ...)
}

par(mfrow=c(1,4))
show(y, length(rows), length(cols), main="Data")

y.res <- matrix(residuals(fit), nrow=length(rows))
show(y.res, length(rows), length(cols), hillshade=FALSE, main="Residuals")
#hist(y.res, main="Histogram of Residuals", ylab="", xlab="Value")

# Increase the `8` to find fewer local outliers; decrease it to find more.
sigma <- 8 * diff(quantile(y.res, c(1/4, 3/4)))
mu <- median(y.res)
outlier <- abs(y.res - mu) > sigma
cat(sum(outlier), "outliers found.\n")

# Fix up the data (impute the values at the outlying locations).
y.imp <- matrix(predict(fit), nrow=length(rows))
y.imp[outlier] <- y[outlier] - y.res[outlier]

show(y.imp, length(rows), length(cols), main="Imputed")
show(y.real, length(rows), length(cols), main="Real")

引用自:https://stats.stackexchange.com/questions/100282

comments powered by Disqus