加速 R 循环

Question

我是这个论坛的新手，也是 R 的初学者，如果我的脚本/问题以令人困惑的方式编写，我深表歉意。 我有 3 年内来自大约 130 个不同站点的天气数据，我想填补这些空白。现在，我只为全球辐射这样做，但我还有四个变量。 我的桌子是这样的：

tbl <- read.table(text = 
"    Date.and.Time I glbRad I precipitation.mm.day I rel.hum... I wind.speed..m.s I temperature.
    1 I 2010-01-01-01 00:00:00 I 0.6 I 0.1 I 99.6 I 1 I 2.3 
    2 I 2010-01-01-01 01:00:00 I 0.6 I 0 I 99.5 I 1 I 2.2 ", 
sep = "I", header = TRUE)

如果间隔只持续一两个小时，我会取之前和之后的测量值的平均值。 如果间隔持续时间超过两个小时，我将使用最近的气象站的值，该气象站在所需时间段内拥有完整数据。我有一个表 distanzen.csv，其中包含第一列中气象站的 name_i 和相邻气象站的 name_j。相邻站点按距离排序。

neighbors <- read.table(header = TRUE, sep = "I",stringsAsFactors = FALSE,text = 
"name_j I name_i I distance
1 I Ainersthofen I Edelshausen I 16.303
2 I Ainersthofen I Gablingen I 19.684")

一般来说，该脚本有效。但这太慢了。你知道我怎样才能加快速度吗？我知道我应该以某种方式摆脱循环，但我不太知道该怎么做。 此外，如果相邻站的日期完全丢失（整行丢失），我会收到错误“参数长度为零”。在这种情况下，我想选择第二近的邻居。

#reading data
file_path="F:/SkriptAktion/wetter_csv_spalten_richtig_Ortsnamen/"
setwd(file_path)
names <-list.files()
d =1
for (n in names){
  table<-read.csv(paste(file_path,n, sep=""), sep=",", header=TRUE, stringsAsFactors=FALSE)

  #change date format
  date <- as.POSIXlt(table$Date.and.Time, tz="utc", format="%d.%m.%Y %H:%M")
  table$Date.and.Time<-date

  #add a column "gaps_radiation" where A) it says „ok“ if the value is not missing  B) it says „MW“ if one or two subsequent values are missing C) it says the name of the neighbouring station  if data of the neighbouring station has been used
 # write „MW“ for all missing values
  table$gaps_radiation <- character(nrow(table))
  table$gaps_radiation<-lapply(table[,"glbRad"],function(x) ifelse (x!=".", "ok", "MW"))

  #change global.radiation from character to numeric
    table$glbRad <- as.numeric(table$glbRad)

  # If the gap lasts only one or two hours, I take the average of the previous and the subsequent measurements.
  #1h gap  
  for (i in 2:(length(table$glbRad)-1)){
    if (table$gaps_radiation[i] == "MW" & table$gaps_radiation[i-1]=="ok" & table$gaps_radiation[i+1]=="ok"){
  table$glbRad[i] <- (table$glbRad[i-1]+table$[i+1])/2
}else {
  #if ((table$gaps_radiation[i] == "MW"){(table$gaps_radiation[i] == "MW"}
  table$glbRad[i] <- table$glbRad[i]
}
  }

  #2h gap

  for (i in 3:(length(table$glbRad)-1)){
    if (table$gaps_radiation[i] == "MW" 
        & table$gaps_radiation[i-1] == "MW"
        & table$gaps_radiation[i-2] == "ok"
        & table$gaps_radiation[i+1]=="ok"){
      table$glbRad[i] <- (table$glbRad[i-2]+table$glbRad[i+1])/2
  table$glbRad[i-1] <- (table$glbRad[i-2]+table$glbRad[i+1])/2
}else {table$glbRad[i] <- table$glbRad[i]
}  
  }


   # gaps in the beginning/end of table
  # 1h gap

  if (table$gaps_radiation[length(table$glbRad)]== "MW" & table$gaps_radiation[length(table$glbRad)-1]=="ok"){
table$glbRad[length(table$glbRad)] <- table$glbRad[length(table$glbRad)-1]
  }else {table$glbRad[length(table$glbRad)] <- table$glbRad[length(table$glbRad)]
  }

  if (table$gaps_radiation[1]== "MW" & table$gaps_radiation[2]=="ok"){
table$glbRad[1] <- table$glbRad[2]
  }else {table$glbRad[1] <- table$glbRad[1]
  }

  # 2h gap

  if (table$gaps_radiation[length(table$glbRad)]== "MW" & table$gaps_radiation[length(table$glbRad)-1] == "MW" & table$gaps_radiation[length(table$glbRad)-2]=="ok"){
table$glbRad[length(table$glbRad)] <- table$glbRad[length(table$glbRad)-2]
table$glbRad[length(table$glbRad)-1] <- table$glbRad[length(table$glbRad)-2]
  }else {table$glbRad[length(table$glbRad)] <- table$glbRad[length(table$glbRad)]
     table$glbRad[length(table$glbRad)-1] <- table$glbRad[length(table$glbRad)-1]
  }

  if (table$gaps_radiation[1]== "MW" & table$gaps_radiation[2] == "MW"& table$gaps_radiation[3]=="ok"){
table$glbRad[1] <- table$glbRad[3]
table$glbRad[2] <- table$glbRad[3]
  }else {table$glbRad[1] <- table$glbRad[1]
     table$glbRad[2] <- table$glbRad[2]
  }


  #gaps > 2h

  mis_dates <- table[(is.na(table$glbRad)),"Date.and.Time"]
  if (length(mis_dates)>=1){

neighbours <- read.csv(file="F:/SkriptAktion/distanzen.csv", header=TRUE, sep=",", dec=".", fill=TRUE, stringsAsFactors=FALSE)


tab1 <- read.csv(file=paste(file_path, neighbours$name_j[d*130+1], ".csv", sep=""), sep=",", header=TRUE, stringsAsFactors=FALSE)
tab1$Date.and.Time <- as.POSIXlt(tab1$Date.and.Time, tz="utc",format="%d.%m.%Y %H:%M")
tab1$glbRad <- as.numeric(tab1$glbRad)

for (i in 1:length(mis_dates)){
  table[table$Date.and.Time == mis_dates[i], "glbRad"] <- tab1[tab1$Date.and.Time == mis_dates[i], "glbRad"]
  table[table$Date.and.Time == mis_dates[i],"gaps_radiation"] <- neighbours$name_j[d*130+1]}

if (nrow(table[is.na(table$glbRad),])>0) {  
  tab1 <- read.csv(file=paste(file_path, neighbours$name_j[d*130+2], ".csv", sep=""), sep=",", header=TRUE, stringsAsFactors=FALSE)
  tab1$Date.and.Time <- as.POSIXlt(tab1$Date.and.Time, tz="utc",format="%d.%m.%Y %H:%M:%S")
  for (i in 1:length(mis_dates)){
    table[table$Date.and.Time == mis_dates[i], "glbRad"] <- as.numeric(tab1[tab1$Date.and.Time == mis_dates[i], "glbRad"])
    table[table$Date.and.Time == mis_dates[i],"gaps_radiation"] <- neighbours$name_j[d*130+2]}
}else {table <- table}

if (nrow(table[is.na(table$glbRad),])>0) {  
  tab1 <- read.csv(file=paste(file_path, neighbours$name_j[d*130+3], ".csv", sep=""), sep=",", header=TRUE, stringsAsFactors=FALSE)
  tab1$Date.and.Time <- as.POSIXlt(tab1$Date.and.Time, tz="utc",format="%d.%m.%Y %H:%M:%S")
  for (i in 1:length(mis_dates)){
    table[table$Date.and.Time == mis_dates[i], "glbRad"] <- tab1[tab1$Date.and.Time == mis_dates[i], "glbRad"]
    table[table$Date.and.Time == mis_dates[i],"gaps_radiation"] <- neighbours$name_j[d*130+3]}
}else {write.table(table,paste("F:/SkriptAktion/Lueckenfueller_radiation/", n, sep=""),sep=",", row.names=FALSE, col.names=TRUE, na="")}

if (nrow(table[is.na(table$glbRad),])>0) {
  write.table(table,paste("F:/SkriptAktion/Lueckenfueller_radiation/", "lueckig", n, sep=""),sep=",", row.names=FALSE, col.names=TRUE, na="")
}else {table <- table}

  }else {write.table(table,paste("F:/SkriptAktion/Lueckenfueller_radiation/", n, sep=""),sep=",", row.names=FALSE, col.names=TRUE, na="")}
  d<- d+1   
}

Answer 1

我认为内部循环可以很容易地矢量化，您只需要注意索引，因为您不想使用第一个和最后一个元素。

i <- 2:(length(table$global.radiation..W.qm.) -1)
i <- 1 + which(table$gaps_radiation[i] == "MW" & table$gaps_radiation[i-1]=="ok" & table$gaps_radiation[i+1]=="ok")
table$global.radiation..W.qm.[i] <- (table$global.radiation..W.qm.[i-1]+table$global.radiation..W.qm.[i+1])/2

Answer 2

这是问题的一个版本的解决方案，符号更简洁。

您可能想要做的第一件事是将所有表同时加载到内存中，而不是在执行过程中重新加载相邻表。接下来我要做的是创建一个数据框或矩阵，其中每列包含给定站点的观测值，每一行对应于给定时间，并且在没有数据的地方缺少值（即NA）。这使得数据更容易操作，而不是放在不同的表中。很容易做到，而且你把这种格式的观察数据清洗干净后，也很容易把清洗后的数据放回原来的表中。

我将生成此格式的示例数据集，如下所示：

> set.seed(100)
> numID<-5
> data<-list()
> for(i in 1:numID){
+     data[[letters[i]]]<-rnorm(1000)
+     data[[i]][sample(998,200)+1]<-NA
+ }
> data<-data.frame(data)
> head(data,10)
            a          b            c           d          e
1  -0.5021924  0.1832545  0.465130835 -0.41210403 -0.8573020
2   0.1315312 -1.4173952  1.301940661          NA  0.9045634
3          NA  0.7547373 -0.427443347  0.02168948  0.8159008
4   0.8867848  0.8888487           NA -1.01383931 -1.1543267
5   0.1169713 -0.6939272 -0.540616369  0.42388204  0.4156978
6          NA -1.8599799  1.038092588 -0.75247680 -1.0199797
7          NA  0.3463114  0.714788709  2.00850576  0.2821374
8   0.7145327         NA           NA  0.81969681         NA
9          NA         NA           NA -1.14063105 -1.0967526
10 -0.3598621         NA -0.009403063          NA  0.9392961

所以这里的电台名称是abcde。为了便于说明，我将邻居列表数据表示为一个矩阵，其中包含每个站点的邻居作为行，按距离排序：

< nn<-c()
> for(i in 1:numID){
+   nn<-rbind(nn,sample(letters[(1:numID)[-i]]))
+ }
> rownames(nn)<-names(data)
> nn
  [,1] [,2] [,3] [,4]
a "c"  "b"  "e"  "d" 
b "c"  "d"  "e"  "a" 
c "a"  "d"  "e"  "b" 
d "c"  "e"  "b"  "a" 
e "c"  "b"  "d"  "a" 

所以这意味着离a最近的车站是c，其次是b等

好的，所以首先要用两个小时或更短的间隔的平均值来填充缺失值。这可以通过调用apply 和辅助函数来完成，无需任何循环：

> fillAve<-function(x){
+   w<-which(!is.na(x))
+   d1<-w[which(diff(w)==2)]
+   d2<-w[which(diff(w)==3)]
+   x[d1+1]<-(x[d1]+x[d1+2])/2
+   x[d2+2]<-x[d2+1]<-(x[d2]+x[d2+3])/2
+   x  
+ }
> data2<-data.frame(lapply(data,fillAve))
> head(data2,10)
            a          b            c           d          e
1  -0.5021924  0.1832545  0.465130835 -0.41210403 -0.8573020
2   0.1315312 -1.4173952  1.301940661 -0.19520727  0.9045634
3   0.5091580  0.7547373 -0.427443347  0.02168948  0.8159008
4   0.8867848  0.8888487 -0.484029858 -1.01383931 -1.1543267
5   0.1169713 -0.6939272 -0.540616369  0.42388204  0.4156978
6   0.4157520 -1.8599799  1.038092588 -0.75247680 -1.0199797
7   0.4157520  0.3463114  0.714788709  2.00850576  0.2821374
8   0.7145327         NA  0.352692823  0.81969681 -0.4073076
9   0.1773353         NA  0.352692823 -1.14063105 -1.0967526
10 -0.3598621         NA -0.009403063 -0.67154451  0.9392961

最后，对于仍然缺失的值，我们需要搜索最近的邻居来填充它们。下面的实现使用了两个循环，但是这些循环只遍历一组位置，而不是整个数据集，所以性能打击不会太差。尝试从 R 代码中删除循环时要记住的一件事是，最重要的是避免在 R 中执行最内层的循环，在这种情况下，循环将遍历各个观察值。只要适当地对内部循环进行矢量化，在 R 中在外部级别使用显式循环通常不会对性能造成太大影响，因为每次外部迭代通常都会做足够多的工作，相比之下，R 的解释机制的开销会很小。

所以这段代码将生成一个新的数据框data3，其中包含干净的数据：

> data3<-data2
> for(i in 1:numID){  
+   x<-data3[[i]]
+   w<-which(is.na(x))
+   j<-1
+   while(length(w)>0 && j<=(numID-1)){
+     y<-data2[[nn[i,j]]]
+     x[w]<-y[w]
+       w<-which(is.na(x))
+       j<-j+1
+   }
+   data3[[i]]<-x
+ }
> head(data3,10)
            a            b            c           d          e
1  -0.5021924  0.183254452  0.465130835 -0.41210403 -0.8573020
2   0.1315312 -1.417395156  1.301940661 -0.19520727  0.9045634
3   0.5091580  0.754737319 -0.427443347  0.02168948  0.8159008
4   0.8867848  0.888848672 -0.484029858 -1.01383931 -1.1543267
5   0.1169713 -0.693927195 -0.540616369  0.42388204  0.4156978
6   0.4157520 -1.859979946  1.038092588 -0.75247680 -1.0199797
7   0.4157520  0.346311411  0.714788709  2.00850576  0.2821374
8   0.7145327  0.352692823  0.352692823  0.81969681 -0.4073076
9   0.1773353  0.352692823  0.352692823 -1.14063105 -1.0967526
10 -0.3598621 -0.009403063 -0.009403063 -0.67154451  0.9392961