メモ

• ２項分布の連続化：β分布
• 多項分布の連続化：ディリクレ分布
  • 確率ベクトル、観測度数,
  • 確率
• に制約をつけつつ、確率を最大にするようなが知りたいとする
• Rのoptim()関数で最適化することにする
• 少し普通ではない変数の置き方をする
  • 長さのベクトルを考える
  • これをある回転ベクトルXで回転したとしてなるベクトルが得られる
  • の個の要素のうち、個は可変変数とし、残りの個は不変とするような制約を考える

# 単位ベクトルfを作る
k <- 6 # 次元
f <- runif(k)
f <- f/sum(f)
# すべての要素が非負のk次元ベクトルを色々作る
library(MCMCpack)
V <- c(rdirichlet(1,rep(1,k)))*100
# 回転行列を作る
library(GPArotation)
X <- Random.Start(k)
# 回転する
U <- X%*%V
# k'個を決めて、可変要素を決める
#k_pr <- sample(1:(k-1),1)
#frees <- sample(1:k,k_pr)
k_pr <- k-1
frees <- 2:k
non.frees <- (1:k)[-frees]

# 以下でlgamma()関数を使いたいのだがlgamma()の引数が負になると困る…
# それを回避する「正しいかどうか不明な」工夫
my.lgamma <- function(x){
	ret <- rep(0,length(x))
	for(i in 1:length(x)){
		if(x[i] >= 0){
			ret[i] <- lgamma(x[i])
		}else{
			ret[i] = exp(-x[i])
		}
	}
	ret
}

# 確率を最大にする(確率の対数を最大にする)可変変数は長さk_prのベクトルinit.var
# 上で決めてきた変数を使う
# optim()は最小化なので、確率の対数の負数を最小化することにする
my.fx <- function(init.var){
 new.U <- U
 new.U[frees] <- init.var
 # 回転行列Xの逆行列はt(X)
 new.V <- t(X)%*%new.U
 # lgamma()に計算不能範囲があるので、うまく行かないのだが…
 #if(min(new.V) < 0){
  #return(.Machine$double.eps)
 #}else{
  #non.zeros <- which(new.V >= 0)
  #non.zero.new.V <- new.V[non.zeros]
  #return(-(my.lgamma(sum(new.V)+1)-sum(my.lgamma(new.V+1)) + sum(new.V * log(f))))
  return(-(my.lgamma(sum(new.V)+1)-sum(my.lgamma(new.V+1)) + sum(new.V * log(f))))
  #return(-(lgamma(sum(non.zero.new.V)+1)-sum(lgamma(non.zero.new.V+1)) + sum(non.zero.new.V * log(f[non.zeros]))))
 #}
 
}
# 初期値を変えて何回かやってみる
n.iter <- 10
for(i in 1:n.iter){
	init.var <- U[frees]+rnorm(k_pr)*0.0001
	#optim.out <- optim(init.var,my.fx,method = "L-BFGS-B",lower =0)
	optim.out <- optim(init.var,my.fx)

	#print(optim.out$par)
	print(optim.out$value)
	opt.U <- U
	opt.U[frees] <- optim.out$par
	opt.V <- t(X) %*% opt.U
	print(c(opt.V))
}

• 追加のメモ

R <- 3

N <- 3

f <- runif(R)
f <- f/sum(f)
n.iter <- 10000
X <- rdirichlet(n.iter,rep(1,R)) * N 

v <- lgamma(apply(X,1,sum)+1) - apply(lgamma(X+1),1,sum) + apply(t(X) * log(f),2,sum)
v.2 <- apply((t(X)-f*N)^2,2,sum)


d <- apply((t(X)-f*N)^2,2,sum)

plot(d,v.2)
cols <- (v * 0.1) %%3

plot3d(X,col=cols)

library(rgl)
library(MCMCpack)
R <- 3

N <- 1000

f <- runif(R)
f <- f/sum(f)
#f <- rep(1/R,R)
n.iter <- 10000
X <- rdirichlet(n.iter,rep(1,R)) * N 
X <- X + matrix(rnorm(n.iter*R),ncol=R)*N/10
v <- lgamma(apply(X,1,sum)+1) - apply(lgamma(X+1),1,sum) + apply(t(X) * log(f),2,sum)
v.2 <- apply((t(X)-f*N)^2/(f*N),2,sum)
v.2 <- v.2 - min(v.2)+0.01

d <- apply((t(X)-f*N)^2,2,sum)

plot(v,v.2)
cols <- (v * 0.5) %%3
#cols <- (d * 1.3) %%3

plot3d(X,col=cols)
open3d()
plot3d(X,col=(d * 0.0001) %%3)

plot3d(X,col=(v.2/(d)) %% 3)

plot3d(X,col=gray((max(d/v.2)-d/v.2)/(max(d/v.2)-min(d/v.2))))