Halton列の作成をするためのRコード

Disclaimer

・この記事は筆者の所属する組織・団体とは一切関係がありません。
・この記事の内容を利用することによって生じた損失等について、筆者は一切責任を負いません。

Halton列の生成をするRの関数 Halton(N,r)

Quasiモンテカルロ（準モンテカルロ）で利用する機会のある超一様分布列(Low-Discrepancy Sequence; LDS)、の一種であるHalton列の生成コードをここに置いておきます。

## r進数に基づくHalton列の生成コード
trans <- function(n, r){
  if(n <= 0){
    return(NULL)
  }else{
    return(append(Recall(n%/%r,r), n%%r))
  }
}
Halton <- function(N,r){
  a <- numeric()
  for(i in 1:N){
    b <- rev(trans(n = i,r = r))
    a[i] <- sum(2^(-1*(1:length(b))) *b) 
  }
  return(a)
}
### Halton(N,r)はr進数に基づき、Halton列をN個まで生成する。
### 出力はベクトル。

使い方：例えば2進数に基づくHalton列を横軸、3進数に基づくHalton列を縦軸にして、256個の点列をプロットすると以下のようになります。

N <- 256
plot(Halton(N=N , r = 2),Halton(N=N , r = 3))

これは、単に一様乱数を二次元にマップさせた場合（下図）と比べると、より「幅広く」分布する点列であることが視覚的に確認できます。

set.seed(1)
plot(runif(n=N),runif(n=N))

参考文献

Dick, Josef, Frances Y. Kuo, and Ian H. Sloan. "High-dimensional integration: the quasi-Monte Carlo way." Acta Numerica 22 (2013): 133-288.