2014-09-28

数式を知らない生物のパターン認識

R 尤度ベイズ

生物は、別に直線とは $y=ax+b$ であることを知っていたりするわけではなく、「ありふれたインプット」を蓄積して、新たなインプットをそれと照合していることだろう
たとえば二次元視覚入力の場合、2次元平面に光刺激のパターンを受け取り、それが頻回ならば、それをパターンとして保持することはできるだろう。その「記録」は視野細胞の格子に換算すると二次元多項分布
二次元多項分布からディリクレ分布を推定するわけではないけれど、それっぽいことは乱数を使ったりして行うことは「プロセスとしては可能」
なので、新たな網膜入力が、「保持している多項分布由来のディリクレ分布から得られる多項分布の尤度」もプロセスとして作れる
そんな具合に「観測データ」が「蓄積としてのパターン」に対する尤度を計算してみよう
これまでに溜めたすべてのパターンで尤度がイマイチなら、それを「新たなパターン」として蓄えておくこともよい作戦

# n次元格子空間にzなる離散観測多項分布をもたらすような
# 実体Zがあるとする
# 今、別の観察xなる離散観測多項分布が、この実体Zによるものであるかどうかを
# 尤度で算出することとする

# zからはZの各項の割合はディリクレ分布と考える
# このとき、各項の観察確率の期待値は(z[i]+1)/(sum(z)+length(z))である
# これを基に、観察多項分布が観察される確率を計算する

my.dirichlet.likelihood <- function(x,z,log=TRUE){
	Z.exp <- (z+1)/(sum(z)+length(z))
	ret <- sum(x * Z.exp)
	if(!log){
		ret <- exp(ret)
	}
	ret
}

# Zが平均m、分散sのn=1次元正規分布の場合
# 実体情報の精度が高いときは、zの点の数が多い

n.z <- 10000
m <- 3.4
s <- 2.1
z <- matrix(rnorm(n.z,m,s),ncol=1)
#zz <- cbind(z,z)
# n次元領域を区域に分け方を格納する
# 1次元の定幅グリッド
my.grid <- function(x,k=100,alpha=0.1){
	rg <- range(x)
	seq(from=rg[1]-(rg[2]-rg[1])*alpha,to=rg[2]+(rg[2]-rg[1])*alpha,length=k)
}
K <- list()
K[[1]] <- my.grid(zz[,1])
#K[[2]] <- my.grid(zz[,2])

# n次元座標行列と次元別グリッドベクトルから、番地を出す
my.grid.hist <- function(x,K){
	y <- x
	for(i in 1:length(K)){
		tmp <- outer(K[[i]],x[,i],"<")
		y[,i] <- apply(tmp,2,sum)
	}
	as.data.frame(apply(y,2,as.factor))
}
h.z <- my.grid.hist(z,K)

# 観測データ
n.x <- 100
m.x <- 3.1
s.x <- 2.5
x <- matrix(rnorm(n.x,m.x,s.x),ncol=1)
h.x <- my.grid.hist(x,K)

# 観測と理論とを合わせて、そのうえで
h.zx <- rbind(h.z,h.x)
zx <- c(rep("z",n.z),rep("x",n.x))

h.zx2 <- cbind(h.zx,zx)

tab <- table(h.zx2)

total.cell.n <- prod(sapply(K,length)+1)

obs <- rep(0,total.cell.n)
obs[1:length(tab[,1])] <- tab[,1]
dens <- rep(0,total.cell.n)
dens[1:length(tab[,2])] <- tab[,2]
dens.st <- (dens+1)/(sum(dens)+length(dens))


dmultinom(obs,prob=dens.st,log=TRUE)

ms <- seq(from=0,to=4,length=50)
ss <- seq(from=0,to=4,length=50)
ss <- ss[-1]

like <- matrix(0,length(ms),length(ss))
n.z <- 1000
for(i in 1:length(ms)){
	for(j in 1:length(ss)){
		tmp.z <- matrix(rnorm(n.z,ms[i],ss[j]),ncol=1)
		tmp.h.z <- my.grid.hist(tmp.z,K)
# 観測と理論とを合わせて、そのうえで
h.zx <- rbind(tmp.h.z,h.x)
zx <- c(rep("z",n.z),rep("x",n.x))

h.zx2 <- cbind(h.zx,zx)

tab <- table(h.zx2)

total.cell.n <- prod(sapply(K,length)+1)

obs <- rep(0,total.cell.n)
obs[1:length(tab[,1])] <- tab[,1]
dens <- rep(0,total.cell.n)
dens[1:length(tab[,2])] <- tab[,2]
dens.st <- (dens+1)/(sum(dens)+length(dens))


like[i,j] <- dmultinom(obs,prob=dens.st,log=TRUE)


	}
}
image(ms,ss,like)