2013-03-10

乱数で多項式計算

決断分布論理ゲート Stochastic logic gate Biocomputing

昨日の続き
$f(t) = \sum_{i=0}^n a_i t^i$ という関数が $0\le t \le 1$ の範囲で大まかには $[0,1]$ に入っているような場合に、乱数を使って $f(t)$ の値を求めることができる、という話
$f(t) = \sum_{i=0}^n a_i t^i = \sum_{i=0}^n b_i^n B_i^n(t)$ という置き換えができる
$B_i^n(t)=\begin{pmatrix} n\\i \end{pmatrix}t^i(1-t)^{n-i}$
$b_i^n$ はその値がわかっていれば、それを使う。わかっていなくても、 $b_i^n$ の確率でランダムに0/1の値を発生するような状況にあれば、0/1の列を発生させたうえで、それを数え上げてもよい。以下のソースではR.bsとして作っている
他方 $B_i^n(t)$ は確率tで0/1の値を取り分けられるような乱数を発生すると、確率的に(n.iter回発生させて、その比率)得ることができる。

# 確率的に値をつくるときの試行(セット)回数
n.iter <- 100000
# 多項式を適当に作る
# 多項式の次数-1 (n-1次多項式 )
n <- 4
as <- sort(rnorm(n))
# 台は0-1の範囲
x <- seq(from=0,to=1,length=100)
# ひとまずこの多項式の値を作る
y <- rep(0,length(x))
for(i in 1:n){
	y <- y + x^(i-1)*as[i]
}
plot(x,y)
# が、多項式の値が0-1の範囲でないといけないので関数の値を補正する
# 補正に合わせて係数も補正する
min.y <- min(y)-0.1
as.1 <- as
as.1[1] <- as.1[1]-min.y
y.1 <- y-min.y
as.1 <- as.1/max(y.1)
y.1 <- y.1/max(y.1)
plot(x,y.1)
# できた多項式の係数からBernstein 形式の係数に変換する
bs <- Power2Bernstein(as.1)
# 全部0-1の範囲に収まるはずなのでそれを確認
bs
# 修正した関数値を計算する
y.p <- rep(0,length(x))
for(i in 1:length(as.1)){
	y.p <- y.p + x^(i-1) * as.1[i]
}
plot(x,y.p,type="l")
# Berntein 形式での関数値を計算する
y.b <- rep(0,length(x))
n <- length(bs)-1
for(i in 1:length(bs)){
	y.b <- y.b + bs[i] * Bernstein(n,(i-1),x)
}
# ２つの形式の値が一致することを確認する
plot(y.b,y.p)
# b_i^nとB_i^nとを計算する
# 乱数を行列に収めて、それを基にする
R.xs <- R.bs <- matrix(0,n.iter,length(bs))
for(i in 1:length(bs)){
	R.bs[,i] <- sample(0:1,n.iter,replace=TRUE,prob=c(1-bs[i],bs[i]))
	
}
# b_i^n
R.bs <- apply(R.bs,2,sum)/n.iter
print(R.bs)
print(bs)
# ある値tのf(t)の値を確率的に出したいので、そのtの値を決める
z <- sample(1:length(x),1)
t <- x[z]
# B_i^nを計算
R.xs <- matrix(sample(0:1,n.iter*(length(bs)-1),replace=TRUE,prob=c(1-t,t)),nrow=n.iter)

R.xs <- apply(R.xs,1,sum)
R.xs <- tabulate(R.xs+1,length(bs))/n.iter
print(R.xs)
# 確率的計算と、ふつうの多項式計算の値がほぼ等しくなることの確認
sum(R.bs*R.xs)
y.1[z]

上記を実施するための諸関数

# Bernstein polynomial
Bernstein <- function(n,m,x){
	if(m<0){
		return(rep(0,length(x)))
	}else if(m > n){
		return(rep(0,length(x)))
	}else{
		#tmp <-lgamma(n+1)-lgamma(m+1)-lgamma(n-m+1) +m * log(1-x) + (n-m) * log(x)
		#return(exp(tmp))
		return(gamma(n+1)/gamma(m+1)/gamma(n-m+1)*x^m*(1-x)^(n-m))
	}
}
# power formの係数ををBernstein formの係数へ
Power2Bernstein <- function(a){
	b <- rep(0,length(a))
	n <- length(a)-1
	for(i in 1:length(b)){
		for(j in 1:i){
			tmp.i <- i-1
			tmp.j <- j-1
			b[i] <- b[i] + factorial(tmp.i)*factorial(n-tmp.j)/factorial(tmp.i-tmp.j)/factorial(n)*a[j]
		}
	}
	b
}
# Bernstein formの係数セットの最高次数を１つ上げて変更する
Bernstein.1plus <- function(b){
	m <- length(b)-1
	ret <- rep(0,length(b)+1)
	ret[1] <- b[1]
	ret[length(ret)] <- b[length(b)]
	for(i in 2:length(b)){
		ret[i] <- (1-(i-1)/(m+1))*b[i]+(i-1)/(m+1)*b[i-1]
	}
	ret
}
# Bernstein formの係数セットを最高次数mに上げて変更する
Bernstein.m <- function(b,m){
	ret <- b
	n <- length(b)-1
	for(i in 1:(m-n)){
		ret <- Bernstein.1plus(ret)
	}
	ret
}