Rのpurrr::mapをmiraiで並列化 – 伊東生態統計研究室

purrrパッケージのヘルプを読んでいると、miraiパッケージと、purrrのin_parallel関数で並列化する方法がのっていました。バージョン1.1.0から追加された機能のようです。まだexperimentalということですが、ためしてみました。

準備

purrrパッケージと、実行例の作成のためにAHMbookパッケージとunmarkedパッケージを、グラフ作成のためggplot2パッケージを読み込みます。

また、set.seed関数で擬似乱数を固定しておきます。

library(purrr)
library(AHMbook)
library(unmarked)
library(ggplot2)

set.seed(123)

データ

AHMbookパッケージのsimNmix関数で、N混合モデルのデータをつくります。サイト数267、調査回数10回、ゼロ過剰なし、1サイトの平均個体数は3、平均発見率0.6のデータを4反復作成します。

nmix_data <- purrr::map(
  1:4,
  \(i) {
    simNmix(nsites = 267, nvisits = 10,
            mean.theta = 1, mean.lam = 3, mean.p = 0.6,
            Neg.Bin = TRUE,
            show.plots = FALSE)
})
## ***** New simulation *****
## 
## No. sites visited:       267 
## No. rep. visits:           10 
## Total no. visits:        2670 
## 
## 
## Naive overdispersion measure (var/mean) for true abundance (N): 1.22 
## Naive overdispersion measure (var/mean) for observed counts (C): 1.14 
## ***** New simulation *****
## 
## No. sites visited:       267 
## No. rep. visits:           10 
## Total no. visits:        2670 
## 
## 
## Naive overdispersion measure (var/mean) for true abundance (N): 1.39 
## Naive overdispersion measure (var/mean) for observed counts (C): 1.22 
## ***** New simulation *****
## 
## No. sites visited:       267 
## No. rep. visits:           10 
## Total no. visits:        2670 
## 
## 
## Naive overdispersion measure (var/mean) for true abundance (N): 1.31 
## Naive overdispersion measure (var/mean) for observed counts (C): 1.19 
## ***** New simulation *****
## 
## No. sites visited:       267 
## No. rep. visits:           10 
## Total no. visits:        2670 
## 
## 
## Naive overdispersion measure (var/mean) for true abundance (N): 1.6 
## Naive overdispersion measure (var/mean) for observed counts (C): 1.35

発見された個体数のヒストグラムを各反復について表示します。観測で発見されなかった個体も含めた個体数をタイトルに入れています。

for (i in 1:4) {
  p <- data.frame(Num = c(nmix_data[[i]]$C)) |>
    ggplot(aes(x = Num)) +
    geom_bar() +
    scale_x_continuous(name = "Observed number of individuals",
                       breaks = seq(0, 30, 5),
                       minor_breaks = 1:30) +
    labs(title = paste("Replication", i, ", Ntotal = ", nmix_data[[i]]$Ntotal))
  plot(p)
}

並列化してモデルあてはめ実行

unmarkedパッケージのunmarkedFramePCount関数で、各反復のデータにN混合効果モデルをあてはめます。

map関数で、繰り返して実行する関数をin_parallel関数の引数に入れます。in_parallelの中で使用する関数は名前空間を明示的に指定し、使用するオブジェクトは引数として渡します。

並列実行の前にmirai::daemons(n = 4)で、4系列を並行実行するようにします。処理が終わったら、mirai::daemons(n = 0)でリセットしておきます。

macOSで実行していますが、アクティビティモニタでプロセスを確認すると、実行中は確かに4つのプロセスが走っていました。

mirai::daemons(n = 4)

fit <- purrr::map(
  1:4,
  in_parallel(
    \(i) {
      d <- unmarked::unmarkedFramePCount(nmix_data[[i]]$C)
      unmarked::pcount(~1 ~1, data = d,
                       K = 200, mixture = "NB")
    }, nmix_data = nmix_data
  )
)

mirai::daemons(n = 0)

結果

purrr::mapの並行実行が本題なのですが、せっかくですので、実行結果を表示します。ここでは、総個体数の推定値を表示します。

for (i in 1:4) {
  print(sum(bup(ranef(fit[[i]]))))
}
## [1] 775.4051
## [1] 742.4653
## [1] 774.5735
## [1] 830.5136

元の値に近い値が推定されたようです。