在 R 中使用嵌套循环进行计算

Question

我有一个包含三列的地理编码数据集：纬度、经度和集群。我计算了集群的平均中心并将结果存储在两个列表 Center_lat 和 Center_lon.

中

现在我想用 Haversine 公式计算从每个观测值 (3000+) 到每个聚类中心 (30) 的距离。得到一个 3000 x 30 的矩阵。

我尝试使用嵌套 for 循环，但我得到的所有集群的距离都相同。这是代码。

for (i in 1:n){
    for (k in 1:c){
    lat1=radians(Geocode[i,1])
    lon1=radians(Geocode[i,2])
    lat2=radians(Center_lat[k,2])
    lon2=radians(Center_lon[k,2])
    }
    R <- 3958.756 # Earth mean radius [miles]
    dist_mat[i,] <- acos(sin(lat1)*sin(lat2) + cos(lat1)*cos(lat2) * cos(lon2-lon1)) * R
  }

我也在考虑使用 lapply 来替代嵌套循环。但是我不确定如何使用该功能...感谢您的帮助。

# Convert to radian
radians = function(theta=0){return(theta * pi / 180)}

# Calculates the geodesic distance from each property to the center of it's current cluster using the
# Spherical Law of Cosines (slc)
get_dist <- function(lat1, lon1, lat2, lon2) {
  R <- 3958.756 # Earth mean radius [miles]
  d <- acos(sin(radians(lat1))*sin(radians(lat2)) + 
              cos(radians(lat1))*cos(radians(lat2)) * cos(radians(lon2)-radians(lon1))) * R
  return(d) # Distance in miles
}

dist_mat<-lapply()

Answer 1

您似乎想每行每列写入一次矩阵，因此您想要在两个 for 循环中更改矩阵，如下所示：

for (i in 1:n){
    for (k in 1:c){
        lat1=radians(Geocode[i,1])
        lon1=radians(Geocode[i,2])
        lat2=radians(Center_lat[k,2])
        lon2=radians(Center_lon[k,2])
        R <- 3958.756 # Earth mean radius [miles]
        dist_mat[i,k] <- acos(sin(lat1)*sin(lat2) + cos(lat1)*cos(lat2) * cos(lon2-lon1)) * R
    }
}

Answer 2

这是您要在 R 中矢量化的计算类型。这里我们使用 outer 从您的 Geocode 数据和 Center_x 数据生成所有可能的行索引组合，然后一举应用距离函数。

首先，以更易于使用的形式获取数据（一个矩阵表示位置，另一个矩阵表示中心，第一列纬度，第二列经度）：

# See Data section below for actual data used

# G <- radians(Geocode)
# C <- radians(cbind(Center_lat[, 2], Center_lon[, 2]))    
R <- 3958.756 # Earth mean radius [miles]

定义函数，注意我们如何使用索引在 G 和 C 中查找实际坐标，以及函数是如何向量化的（即我们只需要调用它一次所有数据):

my_dist <- function(xind, yind)
  acos(
    sin(G[xind, 1]) * sin(C[yind, 1]) + 
    cos(G[xind, 1]) * cos(C[yind, 1]) * cos(C[yind, 2] - G[xind, 2])
  ) * R

并应用 outer:

DISTS <- outer(seq.int(nrow(G)), seq.int(nrow(C)), my_dist)

str(DISTS)
# num [1:3000, 1:30] 4208 6500 8623 7303 3864 ...

quantile(DISTS) # to make sure stuff is reasonable:
#     0%       25%       50%       75%      100% 
#  0.000  4107.574  6204.799  8333.155 12422.059

这在我的系统上运行大约 30 毫秒。

数据：

set.seed(1)
lats <- runif(10000, -60, 60) * pi / 180
lons <- runif(10000, -179, 180) * pi / 180

G.ind <- sample(10000, 3000)
C.ind <- sample(10000, 30)

G <- cbind(lats[G.ind], lons[G.ind])
C <- cbind(lats[C.ind], lons[C.ind])

在 R 中使用嵌套循环进行计算

Calculation with Nested loops in R

r

nested-loops