2013-07-14

関連を絵で見る

トーラス上の乱点をvon Mises Fisher乱点を使って作る
そのトーラス上乱点の方向ベクトルを期待ベクトルとしたvon Mises Fisher乱点方向（だがノルムはばらばら）であるような点があるとする
それらの相互関係を図に示してみる

library(FNN)

my.dist.KL <- function(X,method="euclidean",k=10){
	if(is.matrix(X)){
		ret <- as.matrix(dist(X,method=method))
		return(ret)
	}else{
		n <- length(X)
		num.sample <- sapply(X,nrow)
		K <- min(k,min(num.sample)-1)
		ret <- array(0,c(n,n,K))
		for(i in 1:n){
			for(j in 1:n){
				ret[i,j,] <- KL.divergence(X[[i]],X[[j]],k=K)
			}
		}
		return(ret)
	}
}


my.perm.KL <- function(X,Y,n.iter=100,method="euclidean",k1=10,k2=10){
	D.X <- my.dist.KL(X,k=k1)
	log.D.X <- log(D.X)
	D.Y <- my.dist.KL(Y,k=k1)
	log.D.Y <- log(D.Y)
	#k-NN.X
	k.NN.X <- t(apply(D.X,1,order))
	k.NN.Y <- t(apply(D.Y,1,order))
	K <- k2
	n <- length(D.X[,1])
	St.K <- rep(0,K)
	for(i in 1:K){
		St.K[i] <- sum(log.D.Y[cbind(1:n,k.NN.X[,i+1])])
		St.K[i] <- St.K[i] + sum(log.D.X[cbind(1:n,k.NN.Y[,i+1])])
	}
	
	St.Perm <- matrix(0,n.iter,K)
	for(i in 1:n.iter){
		sh <- sample(1:n)
		sh.log.D.Y <- log.D.Y[sh,]
		sh.log.D.Y <- sh.log.D.Y[,sh]
		for(j in 1:K){
			St.Perm[i,j] <- sum(sh.log.D.Y[cbind(1:n,k.NN.X[,j+1])])
		}
		sh <- order(sh)
		sh.log.D.X <- log.D.X[sh,]
		sh.log.D.X <- sh.log.D.X[,sh]
		for(j in 1:K){
			St.Perm[i,j] <- St.Perm[i,j] + sum(sh.log.D.X[cbind(1:n,k.NN.Y[,j+1])])
		}

	}

	plot(c(sort(St.Perm[,1]),St.K[1]))
	abline(h = St.K[1])
	ps <- rep(0,K)
	for(i in 1:K){
		ps[i] <- length(which(St.Perm[,i]< St.K[i]))/n.iter
	}
	
	return(list(St.K=St.K,St.Perm=St.Perm,p.value=ps))
}

library(Rcpp)
library(RcppArmadillo)
sourceCpp("RvonMisesFisher.cpp")

d1 <- 2
d2 <- 2
r1 <- 1
r2 <- 0.3
mu1 <- c(1,rep(0,d1-1))
mu2 <- c(1,rep(0,d2-1))
k1 <- 1.2
k2 <- 1.8
n <- 300
out1 <- RvonMisesFisher(n,mu1,k1)

out2 <- RvonMisesFisher(n,mu2,k2)


torus.coords <- function(x1,x2){
	r1 <- sqrt(apply(x1^2,1,sum))
	r2 <- sqrt(apply(x2^2,1,sum))
	
	X <- matrix(0,length(x1[,1]),length(x1[1,])+length(x2[1,])-1)
	X[,1:length(x1[1,])] <- x1
	X <- X * (1+x2[,1])
	X[,(length(x1[1,])+1):length(X[1,])] <- x2[,2:length(x2[1,])]
	X
}

x1 <- out1[[1]]*r1
x2 <- out2[[1]]*r2

theta <- pi/2
M <- matrix(c(cos(theta),sin(theta),-sin(theta),cos(theta)),byrow=TRUE,2,2)

x2 <- x2 %*% M

my.torus <- torus.coords(x1,x2)
rg <- range(my.torus)
my.torus <- rbind(my.torus,rep(rg[1],length(my.torus[1,])))
my.torus <- rbind(my.torus,rep(rg[2],length(my.torus[1,])))
library(rgl)
plot3d(my.torus)

X <- my.torus
Y <- matrix(0,length(X[,1]),length(X[1,]))
for(i in 1:length(X[,1])){
	Y[i,] <- RvonMisesFisher(1,X[i,]/sqrt(sum(X[i,]^2)),10)[[1]][1,]
	Y[i,] <- Y[i,]*runif(1)
}

Y0 <- Y+2
XY <- rbind(X,Y0)
plot3d(X)
open3d()
plot3d(Y0)
open3d()

plot3d(XY)
for(i in 1:length(X[,1])){
	segments3d(c(X[i,1],Y0[i,1]),c(X[i,2],Y0[i,2]),c(X[i,3],Y0[i,3]))
}
Y0 <- Y*3
XY <- rbind(X,Y0)
plot3d(X)
open3d()
plot3d(Y0)
open3d()

plot3d(XY)
for(i in 1:length(X[,1])){
	segments3d(c(X[i,1],Y0[i,1]),c(X[i,2],Y0[i,2]),c(X[i,3],Y0[i,3]))
}

#include <RcppArmadilloExtensions/sample.h>
// [[Rcpp::depends(RcppArmadillo)]]
using namespace Rcpp ;

// [[Rcpp::export]]

List RvonMisesFisher(int n,NumericVector mu,double k) {
	int d = mu.size();
	NumericMatrix ret(n,d);
	int cnt = 0;
	RNGScope scope;
	
	while(cnt<n){
		NumericVector tmp(d);
		double rsq = 0;
		for(int i=0;i<d;++i){
			tmp[i] = ::Rf_rnorm(0,1);
			rsq += tmp[i]*tmp[i];
		}
		rsq = sqrt(rsq);
		for(int i=0;i<d;++i){
			tmp[i] /= rsq;
		}
		double tmp2 = 0;
		for(int i=0;i<d;++i){
			tmp2 += mu[i]*tmp[i];
		}
		NumericVector tmpr = runif(1);
		if(exp(tmp2*k)/exp(k) > tmpr[0]){
			for(int i=0;i<d;++i){
				ret(cnt,i) = tmp[i];
			}
			cnt ++;
		}
		
	}


	return List::create(
		Named("") = ret
	);

}

/*** R
library(MCMCpack)
d <-2
mu <- c(1,rep(0,d-1))
k <- 10
n <- 1000
out <- RvonMisesFisher(n,mu,k)
tmp <- out[[1]]
tmp <- rbind(tmp,rep(1,d))
tmp <- rbind(tmp,rep(-1,d))
plot(tmp)

d <-3
mu <- c(1,rep(0,d-1))
k <- 10
n <- 1000
out <- RvonMisesFisher(n,mu,k)
library(rgl)
tmp <- out[[1]]
tmp <- rbind(tmp,rep(1,d))
tmp <- rbind(tmp,rep(-1,d))
plot3d(tmp)

*/