Блог

очень простая проблема, но я все еще не могу найти решение:

У меня есть матрица / data.frame с измерениями (температуры поверхности), содержащими NAs.

matrix(c(3,4,NA,NA,5,5,6,7,6,NA,NA,NA),ncol = 3,byrow = T)

Итак, у меня есть NAs внутри «поверхности» и снаружи.

 > matrix(c(3,4,NA,NA,5,5,6,7,6,NA,NA,NA),ncol = 3,byrow = T)
     [,1] [,2] [,3]
[1,]    3    4   NA
[2,]   NA    5    5
[3,]    6    7    6
[4,]   NA   NA   NA
 

Мне нужно интерполировать / экстраполировать недостающие значения, например, с помощью билинейной интерполяции.

Спасибо за вашу помощь, Маттиас

Поскольку у вас есть матрица, вы можете заполнить пробелы по строкам, используя любой подходящий тип регрессии. Для каждой строки вы настраиваете регрессию, где переменная x — это номер столбца, а переменная y — значение ячейки. После этого вы можете заполнить пробелы, просто используя predict вектор номеров столбцов. Затем вы проделываете то же самое по столбцам.

Это оставляет вас с двумя матрицами, в каждой из которых заполнено большинство NA значений. Затем вы берете средний результат для каждой ячейки (отбрасывая NA значения там, где они все еще существуют). Это даст вам интерполированные значения для каждой ячейки (если только у вас нет целых строк и целых столбцов, которые есть NA , и в этом случае NA остается точка пересечения (обратите внимание, что вы можете исправить это, выполнив алгоритм дважды).

Какой метод регрессии вы используете, зависит от базовых данных. Для сетки размером 3 х 3, an lm — это нормально:

 mat <- matrix(c(3,4,NA,NA,5,5,6,7,6,NA,NA,NA),ncol = 3,byrow = T)

rowwise <- t(apply(mat, 1, function(x) {
  if(sum(!is.na(x)) < 2) rep(NA, length(x))
  else predict(lm(x ~ seq_along(x)), newdata = list(x = seq_along(x)))
}))

colwise <- apply(mat, 2, function(x) {
  if(sum(!is.na(x)) < 2) rep(NA, length(x))
  else predict(lm(x ~ seq_along(x)), newdata = list(x = seq_along(x)))
})

colwise[is.na(colwise)] <- rowwise[is.na(colwise)]
rowwise[is.na(rowwise)] <- colwise[is.na(rowwise)]

mat_fixed <- mat
mat_fixed[is.na(mat)] <- ((colwise   rowwise)/2)[is.na(mat)]

mat
#>      [,1] [,2] [,3]
#> [1,]    3    4   NA
#> [2,]   NA    5    5
#> [3,]    6    7    6
#> [4,]   NA   NA   NA

mat_fixed
#>      [,1]     [,2] [,3]
#> [1,] 3.00 4.000000  4.5
#> [2,] 4.75 5.000000  5.0
#> [3,] 6.00 7.000000  6.0
#> [4,] 7.50 8.333333  7.0
 

Вы заметите, что исходные значения из mat не изменились, и мы проделали довольно хорошую работу по интерполяции недостающих значений. Я предпочитаю получить визуальное представление об интерполяции, чтобы мы могли сравнить исходную матрицу:

 library(ggplot2)

ggplot(reshape2::melt(t(mat)), aes(Var1, Var2, fill = value))   
  geom_tile()   
  scale_y_reverse()   
  coord_equal()
 

К интерполированному:

 ggplot(reshape2::melt(t(mat_fixed)), aes(Var1, Var2, fill = value))   
  geom_tile()   
  scale_y_reverse()   
  coord_equal()
 

Для меня это выглядит довольно хорошо.

В зависимости от размера вашей матрицы и т. Д., Возможно, Лучше попробовать интерполирующую функцию, такую как approx , или использовать loess также работает хорошо, хотя они не будут экстраполироваться, поэтому они не гарантируют заполнения всех пропущенных значений по краям вашей матрицы. Если бы у меня была очень большая матрица (размером с изображение), я мог бы использовать approx сначала для интерполяции, а затем lm для экстраполяции, но для матрицы 3 x 3 lm само по себе нормально.

Более полная реализация для интерполяции будет выглядеть примерно так:

 interpolate <- function(y) 
{
  complete_cases <- is.finite(y)
  if(sum(complete_cases) == 0) return(rep(NA, length(y)))
  if(sum(complete_cases) == 1) return(rep(y[complete_cases], length(y)))
  approx(seq_along(y), y, seq_along(y))$y
}
 

И для экстраполяции:

 extrapolate_left <- function(y)
{
  x <- seq_along(y)
  ssOK <- which(!is.na(y))[1:2]
  ssNA <- which(is.na(y))
  ssNA <- ssNA[ssNA < ssOK[1]]
  
  y[ssNA] <- predict(lm(y ~ x, data = data.frame(x, y)[ssOK,]), list(x = x[ssNA]))
  round(y, 6)
}

extrapolate_right <- function(y)
{
  x <- seq_along(y)
  ssOK <- tail(which(!is.na(y)), 2)
  ssNA <- which(is.na(y))
  ssNA <- ssNA[ssNA > ssOK[2]]
  
  y[ssNA] <- predict(lm(y ~ x, data = data.frame(x, y)[ssOK,]), list(x = x[ssNA]))
  round(y, 6)
}
 

Который вы могли бы использовать для заполнения векторов:

 fill_vector <- function(y) {
  y <- interpolate(y)
  x <- seq_along(y)
  if(all(is.na(y))) return(y)
  if(is.na(y[1])) y <- extrapolate_left(y)
  if(is.na(tail(y, 1))) y <- extrapolate_right(y)
  y
}
 

и матрицы:

 fill_matrix <- function(mat) {
  rowwise <- t(apply(mat, 1, fill_vector))
  colwise <- apply(mat, 2, fill_vector)
  colwise[is.na(colwise)] <- rowwise[is.na(colwise)]
  rowwise[is.na(rowwise)] <- colwise[is.na(rowwise)]
  mat[is.na(mat)] <- ((colwise   rowwise)/2)[is.na(mat)]
  mat
}
 

Вы можете видеть, что это работает довольно хорошо, если мы берем матрицу 50 * 50 и заполняем половину ее NA значениями:

 x <- seq(0, pi, length.out = 50)
m <- outer(sin(x), sin(x))
m[sample(length(m), 1250)] <- NA
 

Постройте график со NA значениями:

 ggplot(reshape2::melt(t(m)), aes(Var1, Var2, fill = value))   
  geom_tile()   
  scale_y_reverse()   
  coord_equal()
 

И нанести на график значения, «заполненные» нашим методом fill_matrix :

 ggplot(reshape2::melt(t(fill_matrix(m))), aes(Var1, Var2, fill = value))   
  geom_tile()   
  scale_y_reverse()   
  coord_equal()