超一様分布列(LDS)と積分
超一様分布列というのは、凄い一様に分布してる列のことです(は？)。
これは英語で言うとLow-Descrepancy Sequence なので、略してLDSと呼ばれることもあるようです。
例えば、1/2 , 1/4 , 3/4 , 1/8 , 5/8 , 3/8 , ……
というような数列を第100項まで図示すると、以下のようになります。
これに対して、ただ単に一様乱数を100個とっただけの数列を図示すると、以下のようになります。
比べてみれば一目瞭然ですね、超一様分布列(LDS)の方が、ただの一様分布の列よりも、わざとらしくもっと"一様に"分布しているのがわかります。
こういった超一様分布列を利用することで、積分計算を通常のモンテカルロよりも早く終わらせることができる、可能性があります。便利ですねえ～
LDSを利用したモンテカルロ計算(いやまあ、もはやモンテカルロ計算ではないんですけど)のことを、Quasi -モンテカルロと呼ぶのですが、QMCと略すことにしましょう。
Quasiとは、"疑似"という意味らしいので、日本語訳すると疑似モンテカルロということになりますけど……超一様分布列は「決定論的なただの数列」なので、これを使った積分をモンテカルロ計算と呼ぶとややこしいですね。モンテカルロ法は通常は(疑似)乱数からの積分計算を指しますからね。まあしょうがないですね、そういう用語なんで。
(追記：Quasi-モンテカルロの日本語訳は準モンテカルロらしいです。)
というわけで、今回はLDSを利用した積分をちょっとやってみましょう。
積分をやってみよう！
以下のような関数hを考えます。
関数hは、0<x<1 , 0<y<1 の範囲について、以下のような関数になっています。
(なんかうまくgif作れなくて、長時間眺めてると気持ち悪くなってきそうな図になりました。申し訳ない)
これについて、0<x<1 , 0<y<1 での積分をやってみましょう。
超一様分布列には、Halton列を利用します。
計算結果は以下のようになりました。
一様分布を利用した通常のモンテカルロよりも、Quasiモンテカルロの方が収束が速そうに見えますね。
今回の記事を再現するためのコードを、以下に置いておくので、試しにRで動かしてみてくださ～い。(おわり)
コード
## 設定
f <- function(x){
a <- sin(10*x) * atan(10*x)^2
return(a)
}
curve(f , xlim=c(-1.2,1.2),col="red",lwd=2)
abline(h=0,lwd=2)
g <- function(y){
b <- -cos(2*pi*y) * abs(lgamma(3*y))^asin(y)
return(b)
}
curve(g,col="blue",lwd=2)
abline(h=0,lwd=2)
h <- function(x,y){
c <- (f(x)+g(y))*((x+y)/2+1/2)
return(c)
}
x <- y <- seq(0,1,length=50)
z <- outer(x,y,FUN=h)
persp(x=x,y=y,z=z,
theta = -40,
phi = 20)
## 超一様分布列の利用
# 2,3進数変換
trans2 <- function(n){
if(n <= 0){
return(NULL)
}else{
return(append(Recall(n%/%2), n%%2))
}
}
trans2(1)
trans3 <- function(n){
if(n <= 0){
return(NULL)
}else{
return(append(Recall(n%/%3), n%%3))
}
}
## Halton列の作成
H2 <- function(N){
a <- numeric()
for(i in 1:N){
b <- rev(trans2(i))
a[i] <- sum(2^(-1*(1:length(b))) *b)
}
return(a)
}
plot(H2(100),main="超一様分布列")
plot(runif(100) , main="ただの一様分布列")
H3 <- function(N){
a <- numeric()
for(i in 1:N){
b <- rev(trans3(i))
a[i] <- sum(3^(-1*(1:length(b))) *b)
}
return(a)
}
plot(runif(1000),runif(1000),main="(疑似)乱数列")
plot(H2(1000) , H3(1000) , main="Halton列")
N <- 10000
res1 <- numeric()
res2 <- numeric()
X <- runif(N)
Y <- runif(N)
for(i in 1:N){
res1[i] <- mean(h(X[1:i] , Y[1:i]))
}
plot(res1,type="l")
X <- H2(N)
Y <- H3(N)
for(i in 1:N){
res2[i] <- mean(h(X[1:i] , Y[1:i]))
}
plot(res2,type="l")
plot(res1,type="l",col="red",ylim=c(-0.4,0.4),main="通常のモンテカルロとQuasiモンテカルロの比較")
lines(res2,col=adjustcolor(alpha=0.7,col="blue"))
legend("topright",col=c("red","blue"),legend=c("MC","QMC"),lty=1,lwd=2,cex=1.5)
