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)
数学・コンピュータ関連の姉妹ブログ『ryamadaのコンピュータ・数学メモ』
京都大学大学院医学研究科ゲノム医学センター統計遺伝学分野のWiki
講義・スライド
医学生物学と数学とプログラミングの三重学習を狙う学習ツール
駆け足で読む○○シリーズ
ぱらぱらめくるシリーズ
カシオの計算機
オンライン整数列大辞典