Избегайте нескольких циклов for в R для вычисления матрицы

Итак, в ходе создания поддельных данных для ответа на вопрос о карте я обнаружил, что пишу следующее:

# Generate some fake data
lat <- seq(-90, 90, by = 5)
lon <- seq(-180, 180, by = 10)
phi <- matrix(0, nrow = length(lat), ncol = length(lon))
i <- 1
for (l1 in lat) {
    j <- 1
    for (l2 in lon) {
        phi[i, j] <- (sin(pi * l1 / 180) * cos(pi * l2 / 180))^2
        j <- j+1
    }
    i <- i+1
}
phi <- 1500*phi + 4500  # scale it properly

Теперь очевидно, что эти два центральных цикла for не такие R'ish, как мне хотелось бы. Кажется, я должен иметь возможность получить mapply или что-то в этом роде, но, к сожалению, это возвращает список и на самом деле не делает то, что я хочу. Другие приложения, похоже, тоже не работают.

Что мне здесь не хватает?


r for-loop matrix mapply
person Mike Wise    schedule 24.02.2016    source источник

Ответы (5)


arrow_upward
17
arrow_downward

Вы должны попытаться использовать матричную алгебру. Не нужно использовать какие-либо функции из семейства apply:

lat <- seq(-90, 90, by = 5)
lon <- seq(-180, 180, by = 10)
1500 * tcrossprod(sin(pi * lat / 180), cos(pi * lon / 180))^2 + 4500
person Raad    schedule 24.02.2016
comment
Оригинал был здесь: stackoverflow.com/questions/35592266/ - я дал вам авторство. - person Mike Wise; 24.02.2016
comment
Ваше здоровье! Должен сказать, что приведенные ниже тесты скорости были весьма неожиданными. - person Raad; 24.02.2016
comment
Как и интерес к вопросу. Конечно, сразил меня. - person Mike Wise; 24.02.2016
comment
Хотя мне все еще очень нравится этот ответ, и я использовал его в своей проблеме, теперь я отмечаю, что приведенное ниже решение outer на самом деле более общее. - person Mike Wise; 27.02.2016

arrow_upward
10
arrow_downward

вы можете использовать outer

   x = outer(lat, lon, FUN = function(x,y) {(sin(pi * x/180) * cos(pi * y /180))^2})
    identical(x * 1500 + 4500, phi)
# [1] TRUE

Ответ NBATrends кажется более быстрым, чем другое решение. Вот некоторый бенчмарк

library(microbenchmark) 
microbenchmark(within(df, {
  phi <- (sin(pi * lat / 180) * cos(pi * lon / 180))^2
  phi <- 1500*phi + 4500
}), 1500 * tcrossprod(sin(pi * lat / 180), cos(pi * lon / 180))^2 + 4500, outer(lat, lon, FUN = function(x,y) {(sin(pi * x/180) * cos(pi * y /180))^2}),
((as.matrix(l1)%*%t(as.matrix(l2)))^2) * 1500 + 4500)
Unit: microseconds
                                                                                              expr     min       lq      mean   median       uq     max neval
 within(df, {     phi <- (sin(pi * lat/180) * cos(pi * lon/180))^2     phi <- 1500 * phi + 4500 }) 255.670 262.0095 270.50948 266.6880 277.7060 385.467   100
                                  1500 * tcrossprod(sin(pi * lat/180), cos(pi * lon/180))^2 + 4500  11.471  12.3770  22.30177  12.9805  13.5850 868.130   100
               outer(lat, lon, FUN = function(x, y) {     (sin(pi * x/180) * cos(pi * y/180))^2 }) 137.645 139.7590 144.39520 141.5700 145.1925 179.905   100
                                            ((as.matrix(l1) %*% t(as.matrix(l2)))^2) * 1500 + 4500  16.301  17.6595  20.20390  19.6215  20.5270  80.294   100
person Mamoun Benghezal    schedule 24.02.2016
comment
По многим тоже. Интересно. - person Mike Wise; 24.02.2016
comment
Теперь, когда у меня было время все обдумать, этот ответ (с использованием outer) во многих отношениях является лучшим и более общим ответом, поскольку мы можем поместить произвольную функцию для x и y, а crossprod на самом деле просто выполняет функции, которые являются мультипликативными. товар. crossprod намного быстрее там, где это применимо, и для моей конкретной проблемы crossprod отлично подходит, поэтому я не буду корректировать правильный ответ и оставлю его на этом примечании. - person Mike Wise; 25.02.2016

arrow_upward
7
arrow_downward

Линейная алгебра может быть проще для вашего приложения, потому что вы просто поэлементно умножаете два вектора, что можно сделать через v * u^T. В R умножение матриц равно %*%.

lat <- seq(-90, 90, by = 5)
lon <- seq(-180, 180, by = 10)

l1 <- sin(pi * lat / 180) 
l2 <- s(pi * lon/ 180)

# compute the matrix
phi <- as.matrix(l1)%*%t(as.matrix(l2))
# square each element of the matrix
phi <- phi^2
# scale properly
# square each element of the matrix
phi <- 1500*phi + 4500
person Mathieu B    schedule 24.02.2016

arrow_upward
5
arrow_downward

Зачем привязываться к матричной структуре и применять ее, когда можно векторизовать?

df <- expand.grid(lat = seq(-90, 90, by = 5),
                 lon = seq(-180, 180, by = 10))
df <- within(df, {
  phi <- (sin(pi * lat / 180) * cos(pi * lon / 180))^2
  phi <- 1500*phi + 4500
  })

Вы всегда можете выполнить обратное преобразование, следуя инструкциям здесь.

person sebastian-c    schedule 24.02.2016

arrow_upward
4
arrow_downward

Используя sapply(), но я бы предпочел outer() решение:

#using sapply
phi_1 <- 
  t(
    sapply(lat, function(l1)
      sapply(lon, function(l2)(sin(pi * l1 / 180) * cos(pi * l2 / 180))^2))
  ) * 1500 + 4500

#compare result
identical(phi_1, phi)
# [1] TRUE
person zx8754    schedule 24.02.2016