関連を絵で見る

• トーラス上の乱点を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)
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>
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)

*/
```