メモ

CategoryVector<-function (nc = 3){
	df <- nc - 1
	d <- df + 1
	diagval <- 1:d
	diagval <- sqrt((df + 1)/df) * sqrt((df - diagval + 1)/(df - diagval + 2))
	others <- -diagval/(df - (0:(d - 1)))
	m <- matrix(rep(others, df + 1), nrow = df + 1, byrow = TRUE)
	diag(m) <- diagval
	m[upper.tri(m)] <- 0
	as.matrix(m[, 1:df])
}

make.simplex <- function(k){
	cv <- CategoryVector(k)
	rbind(t(cv*sqrt(1-1/k)),rep(1/sqrt(k),k))
}
make.simplex2 <- function(k){
	cv <- CategoryVector(k)
	t(cv*sqrt(1-1/k))
}

make.simplex.multi <- function(r){
	X <- make.simplex(r[1])
	k <- length(r)
	if(k > 1){
		for(i in 2:k){
			X <- X %x% make.simplex(r[i])
		}
	}
	X
}
make.simplex.multi2 <- function(r){
	X <- make.simplex2(r[1])
	k <- length(r)
	if(k > 1){
		for(i in 2:k){
			X <- X %x% make.simplex2(r[i])
		}
	}
	X
}
n <- 5
r <- c(rep(3,n),2)
r.2 <- r[length(r):1]
X <- make.simplex.multi(r.2)

H <- 2^n
h <- sample(2:(H/2),1)
hap <- sample(0:(H-1),h)
library(MCMCpack)
af <- rdirichlet(2,rep(1,h))
af.m <- apply(af,2,mean)
t <- 0.1/max(abs(af[1,]-af[2,]))
af.ca <- af.m + (af[1,]-af.m)*t
af.co <- af.m + (af[2,]-af.m)*t
max(abs(af.ca-af.co))


ncaco <- c(200,180)
geno <- matrix(sample(1:h,2*sum(ncaco),replace=TRUE),sum(ncaco),2)

from_int.ry <- function(k, x,m=2){
	if(k <= 0){
		return(c())
	} else {
		return(append(Recall(k-1, x%/%m), x%%m))
	}
}
to_int.ry <- function(x,m){
	sum(m^((length(x)-1):0)*x)
}
hap.bin <- matrix(NA,h,n)
for(i in 1:h){
	hap.bin[i,] <- from_int.ry(n,hap[i],2)
}
hap.bin
to_int.ry(hap.bin[1,],2)

hap2geno <- function(h1,h2){
	h.sum <- h1+h2
	tmp <- to_int.ry(h.sum,3)+1
}
geno.2 <- rep(NA,length(geno[,1]))
for(i in 1:length(geno.2)){
	geno.2[i] <- hap2geno(hap.bin[geno[i,1],],hap.bin[geno[i,2],])
}

geno.3 <- matrix(NA,length(geno[,1]),n)
for(i in 1:length(geno.3[,1])){
	geno.3[i,] <- hap.bin[geno[i,1],] + hap.bin[geno[i,2],]
}
nbins <- prod(r)/2
case.table <- tabulate(geno.2[1:ncaco[1]],nbins=nbins)
cont.table <- tabulate(geno.2[(ncaco[1]+1):sum(ncaco)],nbins=nbins)

a <- c(case.table,cont.table)
A <- array(a,r)

non.zero <- which(case.table+cont.table!=0)

case.table2 <- case.table[non.zero]
cont.table2 <- cont.table[non.zero]

#X2 <- make.simplex.multi(c(length(non.zero),2))
X <- make.simplex.multi(c(2,length(non.zero)))

ds <- sample((-2):2,length(non.zero),replace=TRUE)
ds[1] <- -(sum(ds[-1]))

tmp.case.table2 <- case.table2
tmp.cont.table2 <- cont.table2
tmp.case.table2 <- tmp.case.table2 + ds
tmp.cont.table2 <- tmp.cont.table2 - ds


rot.ds2 <- X %*% c(ds,-ds)

plot(rot.ds2)

a <- c(tmp.case.table2,tmp.cont.table2)
r <- c(length(non.zero),2)
A <- array(a,r)
N <- sum(A)
E <- apply(A,1,sum)
for(i in 2:length(r)){
	E <- E %x% t(apply(A,i,sum))/N
}

#A.col <- apply(A,1,sum)
#A.row <- apply(A,2,sum)
#E <- A.col %x% t(A.row)/sum(A.col)
e.vec <- c(E)
apply(E,1,sum)
apply(E,2,sum)

############
df <- length(non.zero)-1

to.use <- 1:df

Xr <- X[to.use,]
Xt <- t(X)
Xtr <- Xt[,to.use]
W <- t(Xtr)%*%diag(1/e.vec)%*%Xtr
W

eigen.W <- eigen(W)
sq.e.val <- sqrt(eigen.W$values)
S <- diag(sq.e.val)
Sinv <- diag(1/sq.e.val)
V <- eigen.W$vectors
V%*%S%*%S%*%t(V) - W

PtoQ <- S%*%t(V)%*%Xr
QtoP <- Xtr%*%V%*%Sinv

###
# table.O

table.O <- A
simplex.conversion <- function(A){
	r <- dim(A)
	r <- r[length(r):1]
	N <-sum(A)
	E <- apply(A,1,sum)
	for(i in 2:length(r)){
		E <- E %x% t(apply(A,i,sum))/N
	}
	e.vec <- c(E)
	Xr <- make.simplex.multi2(r)
	Xtr <- t(Xr)
	W <- t(Xtr)%*%diag(1/e.vec)%*%Xtr
	eigen.W <- eigen(W)
	sq.e.val <- sqrt(eigen.W$values)
	S <- diag(sq.e.val)
	Sinv <- diag(1/sq.e.val)
	V <- eigen.W$vectors

	PtoQ <- S%*%t(V)%*%Xr
	QtoP <- Xtr%*%V%*%Sinv

	return(list(r=r,N=N,E=E,e.vec=e.vec,Xr=Xr,W=W,eigen.W=eigen.W,PtoQ=PtoQ,QtoP=QtoP))
}
sc.out <- simplex.conversion(A)