分布の分布を推定したい

• 昨日と同じ設定。今日は、各観測点に同じ寄与力を授けたうえで、その寄与力を最大にする１パラメタの最適値推定をする方法でやってみる

# 説明変数次元
n.x <- 2
# 被説明変数に影響を及ぼす変数数
n.x.a <- 2
# 被説明変数レベル数
n.y <- 2
Y.level <- 1:n.y
# 標準的な分布を定める
# 説明変数空間分布を定める
# 構成分布の数
n.d <- 3
# 分布の構成比
library(MCMCpack)
r <- rdirichlet(1,rep(1,n.d))
# 分布は正規分布
k.ctr <- 10
k.sd <- k.ctr/2
ctrs <- matrix(runif(n.d*n.x)*k.ctr,ncol=n.x)
sds <- matrix(runif(n.d*n.x)*k.sd,ncol=n.x)

# 確率密度分布関数を定める
# 原点からの距離Lに応じた三角関数 k1 * (sin(k2 * L)+1)/2 + k3 ; k1 <=1, k1+k3<=1
k1 <- runif(1)
k1 <- 1
k2 <- 1
k3 <- runif(1) * (1-k1)

my.calc.p <- function(x,k1,k2,k3){
	k1 * (sin(k2*sqrt(sum(x^2)))+1)/2+k3
}

# サンプリングする
# サンプル総数
n.s <- 50

X <- matrix(0,0,n.x)
Y <- rep(0,n.s)
r.d <- c(rmultinom(1,size=n.s,prob=r))
for(i in 1:n.d){
	tmp <- matrix(0,r.d[i],0)
	for(j in 1:n.x){
		tmp <- cbind(tmp,rnorm(r.d[i],ctrs[i,j],sds[i,j]))
	}
	X <- rbind(X,tmp)
}
P <- apply(X[,1:n.x.a],1,my.calc.p,k1,k2,k3)
P.level <- matrix(0,length(P),n.y)
for(i in 1:n.s){
	tmp.prob <- c(P[i],rep((1-P[i])/(n.y-1),n.y-1))
	tmp.prob <- tmp.prob/sum(tmp.prob)
	P.level[i,] <- tmp.prob
	Y[i] <- sample(1:n.y,1,prob=tmp.prob)
}
plot(X[,1:2],col=Y)

###
# optim推定
x <- X
y <- Y
y.level <- Y.level

my.optim.out <- optim(1,my.optim.r,method="Brent",lower=0,upper=range(x)[2])
est.cnt <- my.r.count(x,y,y.level,my.optim.out$par)
#plot(x[,1:2],col=rgb(P.level[,1],P.level[,2],P.level[,3]))
est.P <- (est.cnt[,1]+1)/(est.cnt[,1]+est.cnt[,2]+2)
est.P.level <- (est.cnt+1)/apply(est.cnt+1,1,sum)
par(mfcol=c(2,2))
plot(P.level,est.P.level,xlim=c(0,1),ylim=c(0,1))
plot(x[,1:2],col=y+1)
rgb1 <- rgb2 <- matrix(0,length(x[,1]),3)
for(i in 1:n.y){
	rgb1[,i] <- 1-P.level[,i]
	rgb2[,i] <- 1-est.P.level[,i]
}

plot(x[,1:2],col=rgb(rgb1[,1],rgb1[,2],rgb1[,3]),pch=20)
plot(x[,1:2],col=rgb(rgb2[,1],rgb2[,2],rgb2[,3]),pch=20)

par(mfcol=c(1,1))


###
# 予測
n.pt <- 500
X.new <- matrix(0,0,n.x)
r.d.new <- c(rmultinom(1,size=n.pt,prob=r))

for(i in 1:n.d){
	tmp <- matrix(0,r.d.new[i],0)
	for(j in 1:n.x){
		tmp <- cbind(tmp,rnorm(r.d.new[i],ctrs[i,j],sds[i,j]))
	}
	X.new <- rbind(X.new,tmp)
}
P.new <- apply(X.new[,1:n.x.a],1,my.calc.p,k1,k2,k3)
P.level.new <- matrix(0,length(P.new),n.y)
for(i in 1:n.pt){
	tmp.prob <- c(P.new[i],rep((1-P.new[i])/(n.y-1),n.y-1))
	tmp.prob <- tmp.prob/sum(tmp.prob)
	P.level.new[i,] <- tmp.prob
}

est.cnt.new <- my.r.count.2(X.new,x,y,y.level,my.optim.out$par)

est.P.new <- (est.cnt.new[,1]+1)/(est.cnt.new[,1]+est.cnt.new[,2]+2)
est.P.level.new <- (est.cnt.new+1)/apply(est.cnt.new+1,1,sum)
par(mfcol=c(2,2))
xlim <- range(c(x[,1],X.new[,1]))
ylim <- range(c(x[,2],X.new[,2]))
plot(P.level.new,est.P.level.new,xlim=c(0,1),ylim=c(0,1))
plot(x[,1:2],col=y,pch=20,xlim=xlim,ylim=ylim)
rgb1.new <- rgb2.new <- matrix(0,length(X.new[,1]),3)
for(i in 1:n.y){
	rgb1.new[,i] <- 1-P.level.new[,i]
	rgb2.new[,i] <- 1-est.P.level.new[,i]
}

plot(X.new[,1:2],col=rgb(rgb1.new[,1],rgb1.new[,2],rgb1.new[,3]),pch=20,xlim=xlim,ylim=ylim)
plot(X.new[,1:2],col=rgb(rgb2.new[,1],rgb2.new[,2],rgb2.new[,3]),pch=20,xlim=xlim,ylim=ylim)

par(mfcol=c(1,1))