特異値分解

  • 特異値分解 singular value decomposition
  • M=UDV^{*}というMの分解のこと。ただし、Dは対角成分以外がゼロで、対角成分は非負。U,V^{*}はユニタリ行列。
  • 特異値分解をするには、Rにはsvd関数がある。
> mtest<-matrix(c(1,2,3,4,5,6),ncol=2)
> mtest
     [,1] [,2]
[1,]    1    4
[2,]    2    5
[3,]    3    6
> svd(mtest)
$d
[1] 9.5080320 0.7728696

$u
           [,1]       [,2]
[1,] -0.4286671  0.8059639
[2,] -0.5663069  0.1123824
[3,] -0.7039467 -0.5811991

$v
           [,1]       [,2]
[1,] -0.3863177 -0.9223658
[2,] -0.9223658  0.3863177
  • 2x3SNPケースコントロールデータを6つの数値のリストとしてsvd関数に渡すとする
svd2x3<-function(x){
m<-matrix(x,ncol=2)
return(svd(m))
}

として

> svd2x3(c(1,2,3,4,5,6))
$d
[1] 9.5080320 0.7728696

$u
           [,1]       [,2]
[1,] -0.4286671  0.8059639
[2,] -0.5663069  0.1123824
[3,] -0.7039467 -0.5811991

$v
           [,1]       [,2]
[1,] -0.3863177 -0.9223658
[2,] -0.9223658  0.3863177
  • ちょっと個人的な用途で、ごちゃごちゃと計算するのを組み込んだ関数をメモしておくと・・・
svd2x3testRT<-function(x,y){
p<-x[1]+x[2]+x[3]
q<-x[4]+x[5]+x[6]
a<-x[1]+x[4]
b<-x[2]+x[5]
c<-x[3]+x[6]
N<-p+q
e<-matrix(c(p*a,p*b,p*c,q*a,q*b,q*c),ncol=2)
e<-e/N
d<-c(cos(y[2]*pi+2/3*pi),cos(y[2]*pi),cos(y[2]*pi-2/3*pi),-cos(y[2]*pi+2/3*pi),-cos(y[2]*pi),-cos(y[2]*pi-2/3*pi))
d<-y[1]*d*N
matd<-matrix(d,ncol=2)
m<-e+d
print(matd)
print(e)
print(m)

s<-svd(m)
print("$u")
print(s$u[,1]*s$u[,1])
print(s$u[,2]*s$u[,2])
print(s$u[1,]*s$u[1,])
print(s$u[2,]*s$u[2,])
print(s$u[3,]*s$u[3,])
print("$v")
print(s$v[,1]*s$v[,1])
print(s$v[,2]*s$v[,2])
print(s$v[1,]*s$v[1,])
print(s$v[2,]*s$v[2,])

return(svd(m))
  • 2x3表の総サンプルサイズがNとすると、特異値分解U,V\{*}は変わらず、DがN倍になる
  • V^{*}は2x2のユニタリ行列であるが、これは、ユニタリである条件(構成する2つの転置ベクトル、v_1{*},v_2{*}のそれぞれのノルムが1で、その内積がゼロであるという条件)から、|v_{11}v_{22}-v_{12}v_{21}|=1という性質がある。
  • ケースとコントロールとで分布に差がないとき、Dの対角成分2つのうちいずれかがゼロである
  • 観測テーブルの行列と期待度数の行列との差の行列の第i,j成分はu_{i,j}\times \sum_{k\ne i,l\ne j}u_{k,l})-u_{i'\ne i,j}\times \sum_{k\ne i',l\ne j} u_{k,l}で表される