## 2015年10月 2日(金曜日) # -*- mode: R -*- ##EmacsとESS使いの方に、モードラインの設定。 ##plot(sin) ##setwd(choose.dir()) ## EIが677個、Nが1478個。全部で2155個。 source("pat.txt")#これは、dump()コマンドで作成した物、 ## #ニューラルネットワークのプログラム source("snn.txt") ##300回学習済みの重みパラメータをセットする。 source("snnWeight.txt") ##エラーの様子。 plot(sapply(1:2155, activate)) abline(h=0.1, col=2) hist(sapply(1:2155, activate),nc=300) hist(log(sapply(1:2155, activate)),nc=300) abline(v=log(0.1), col=2) ##これは、Histgramを重ねるときの技 hist(log(sapply(1:2155, activate)),nc=300) -> b ##BINの位置を教えてもらうだけ。 hist(log(sapply(1:677, activate)), breaks=b$breaks, col=1,main="") par(new=T) hist(log(sapply(678:2155, activate)), breaks=b$breaks, col=2,main="") abline(v=log(0.1), col=3) ##失敗例について、 (1:2155)[sapply(1:2155, activate)>0.1] (dimnames(BufInput)[[2]])[231] ##[1] "CTGGTTCTC" (dimnames(BufInput)[[2]])[sapply(1:2155, activate)>0.1] ## "CTGGTTCTC" ## "AGGGGCCAC" ## "ACAGTGCTG" ## "CTACAACGG" ## "AAGGTGCCA" ## "CAGGTATAT" ## "AAGGTGCCA" ##593はEI,1944はN。 (dimnames(BufInput)[[2]])[593]==(dimnames(BufInput)[[2]])[1944] ##で、中間層の出力は、どうなっているの?? hOut <- c() for(iSet in 1:2155){ activate(iSet); hOut <- cbind(hOut, Hidden$out) } ##実際に使うには、SlidingWindow ##部分配列の切り出し。 substring("ATGC",1,1) ##{開始を1、終了位置を1}とすると、"A"が取れる。 substring("ATGC",1:4,1:4) ##{開始を1、終了位置を1}、{開始を2、終了位置を2}、{開始を3、終了位置を3}、{開始を4、終了位置を4} ##[1] "A" "T" "C" "C" ##1塩基の変換(塩基を実数に) zz<-function(x){switch(EXPR = x, "A"=c(1,0,0,0), "T"=c(0,1,0,0), "G"=c(0,0,1,0), "C"=c(0,0,0,1))} ##長い配列を1つ1つに分解できた。のに、zz関数に適用(apply)する。 sapply(substring("ATGC",1:4,1:4),zz) ## A T C C ## [1,] 1 0 0 0 ## [2,] 0 1 0 0 ## [3,] 0 0 0 0 ## [4,] 0 0 1 1 #で欲しいのは、36個の{0,1}のながーいベクトルへの型なので、 #行列になってしまったのを、ベクトルに(無理やり)変換 as.vector(sapply(substring("ATCC",1:4,1:4),zz)) ## [1] 1 0 0 0 0 1 0 0 0 0 0 1 0 0 0 1 ##ひっくるめて関数化。 ##nchar(x)は、文字列の長さ yy<-function(str){ as.vector(sapply(substring(str,1:nchar(str),1:nchar(str)),zz))}; ##yy("TTCAGCGGCCTCAGCCTGCCTGTCTCC") yy("TTCAGCGG") ##長いテスト文字列が有った。 ##EI, ATRINS-DONOR-521, ##CCAGCTGCATCACAGGAGGCCAGCGAGCAGGTCTGTTCCAAGGGCCTTCGAGCCAGTCTG testSeq <- "CCAGCTGCATCACAGGAGGCCAGCGAGCAGGTCTGTTCCAAGGGCCTTCGAGCCAGTCTG" nchar(testSeq)## 全部で60文字 ##9文字ずつに切って、{0,1}数字列に変換 sapply(substring(testSeq,1:(nchar(testSeq)-8), 9:nchar(testSeq)), yy) rm(BufInput)##一度、クリアにする。 BufInput <- sapply(substring(testSeq,1:(nchar(testSeq)-8), 9:nchar(testSeq)), yy) resp <- c() for(iSet in 1:(nchar(testSeq)-8)){ activate(iSet) Output$out resp <- rbind(resp, Output$out) } matplot(resp, type="b") 1; #中間層の出力を行列の形式で。行方向に15個の中間層の出力、列方向にデータ。 dim(hOut) par(mfrow=c(1,1)) ##相関係数について、おさらい。 ##100個の(-1,1)の一様乱数。 x <- runif(n=100,min=-1,max=1) hist(x) y <- x+runif(100,-0.1,0.1) plot(x,y) ##相関が高い。 cor(x,y) title(sprintf("Cor %5.3f", cor(x,y))) y <- runif(n=100,min=-1,max=1) plot(x,y) ##相関が低い cor(x,y) title(sprintf("Cor %5.3f", cor(x,y))) y <- -x+0.1*runif(n=100,min=-1) plot(x,y) ##相関が高い。 cor(x,y) title(sprintf("Cor %5.3f", cor(x,y))) ##相関係数の計算cor cor(hOut[1,],hOut[2,]) ##ヒストグラム hist(hOut) #中間層の出力のデータに対する、レスポンス。平均と分散。 plot(apply(hOut,2,"mean")) plot(apply(hOut,2,"var")) cor(t(hOut))->a a[row(a)-col(a)>0] hist(abs(a[row(a)-col(a)>0]),nc=20) #####「相関係数の辺りの説明用グラフ」書き専用。 dat1 <- read.table("dat1.txt") plot(dat1,xlim=c(-8.5,2.5),ylim=c(-5,6)) ##主成分分析のグラフ用。角度をいれると射影したヒストグラムが出来る。 poi <- function(angle){ par(mfrow=c(1,2)) plot(dat1,xlim=c(-8.5,2.5),ylim=c(-5,6)) abline(a=mean(dat1[,2])-tan(angle*pi/180)*mean(dat1[,1]),b=tan(angle*pi/180) ,col="blue") d <- matrix(c(cos(angle*pi/180),sin(angle*pi/180)),2,1)##重み付け hist( as.matrix(dat1) %*% d, breaks=seq(-10,6,0.5), main="") title(paste("angel",angle)) } for(i in -90:90 ){ poi(i) dev.copy(jpeg, sprintf("D:/tmp/aa%03d.jpg",i+90), width=640, height=480) dev.off(); } sapply(seq(-90,90,1),poi) par(mfrow=c(1,1)) ##主成分分析 par(mfrow=c(1,1)) prcomp(t(hOut))->b plot((b$sdev)^2,type="b") abline(h=0) #累積寄与率 plot(cumsum((b$sdev)^2/sum((b$sdev)^2)),type="b",ylim=c(0,1)) #クラスター解析 dist(hOut) ##2次元平面にプロット ##library(MASS) ##dat.cmd<-cmdscale(as.dist(1.0-cor(t(hOut)))) #15個の中間層の出力の"距離"を総当たりで求める。 hclust(dist(hOut), method="average") #それに基づいて、階層的クラスタリングを行う。 plot(hclust(dist(hOut), method="average")) ##相関係数の計算cor cor(hOut[3,],hOut[7,]) ##で、相関係数を"距離"と見立てて行うためには、 plot(hclust(as.dist(1.0-cor(t(hOut))), method="average")) ##15次元で分布しているデータを2次元に押しつぶす。 ##押し花。 library(MASS) c <- cmdscale(as.dist(1.0-cor(t(hOut)))) plot(c) ###中間層を1つにしてみる。 #学習結果 source("snnWeight1.txt") Hidden ## $net ## [1] 49.89437 ## $out ## [1] 1 drop(InputHidden$w) (1:2155)[sapply(1:2155, activate)>0.1] ##重みの様子を見てみる。 barplot(drop(InputHidden$w), col=rep(c("red","green","blue","yellow"),9) ) dim(InputHidden$w) matrix(InputHidden$w, nrow=4) barplot(as.vector(rbind(matrix(InputHidden$w, nrow=4),0)), col=c("red","green","blue","yellow","white") ) ##オマケ、lda()でもやってみよう。 library(MASS) i <- as.vector(t(BufInput)) im <- matrix(i,nr=2155) t <- dimnames(BufTarget)[[2]] b <- lda(im,t) plot(b) matplot(predict(b,im)$posterior) hist(predict(b,im)$x) b$scaling barplot(-drop(b$scaling), col=rep(c("red","green","blue","yellow"),9) ) barplot(as.vector(rbind(matrix(b$scaling, nrow=4),0)), col=c("red","green","blue","yellow","white") ) ##数だけから推定。 ##EI境界のもの apply(im[tm[,1]==1,],2,sum) ##そうじゃないのもの apply(im[tm[,1]==0,],2,sum) barplot(apply(im[tm[,1]==1,],2,sum)-apply(im[tm[,1]==0,],2,sum), col=rep(c("red","green","blue","yellow"),9)) barplot(as.vector(rbind(-1*matrix(apply(im[tm[,1]==1,],2,sum)-apply(im[tm[,1]==0,],2,sum), nrow=4),0)), col=c("red","green","blue","yellow","white") ) ##強引に、全ての入力データをクラスタリングしてみる。 ##1つめのデータ im[1,] ## [1] 0 0 0 1 1 0 0 0 0 0 1 0 0 0 1 0 0 1 0 0 0 0 0 1 0 1 0 0 0 0 1 0 0 1 0 0 ##2つめのデータ im[2,] ## [1] 1 0 0 0 0 0 1 0 0 0 1 0 0 0 1 0 0 1 0 0 0 0 1 0 1 0 0 0 0 0 1 0 0 0 0 1 ##差の2乗は、 (im[1,]-im[2,])^2 ##その総和、 sum((im[1,]-im[2,])^2) ##総和の最小値は、0.0。これが、ニューラルネットワークでは「2乗誤差」 sum((im[1,]-im[1,])^2) ##つまり、まったく同じ塩基配列だったら、距離は0。2つの塩基配列中に、異 ##なる塩基が有れば、距離は増える。このとき、塩基の(違いの)種類によらな ##い。異なる塩基数だけが距離の元。 ##全部だと大変 sss <- runif(2155)<0.5 #sss <- runif(2155)<0.2 ##plot(hclust(dist(im[sss,]))) hclust(dist(im[sss,]))->hc ##plot(hc) ##ラベルでEIかそうでないかを知る。 ##hc$orderで、デンドログラム中の葉が判る。それが==1だと、 plot(hc, labels=sapply((tm[sss,1])[hc$order]==1, function(x){ if(x){ return("*");##EI境界にある } else { return(""); } } ) )