特異値分解 - ryamadaの遺伝学・遺伝統計学メモ

特異値分解 singular value decomposition
$M=UDV^{*}$ というMの分解のこと。ただし、 $D$ は対角成分以外がゼロで、対角成分は非負。 $U,V^{*}$ はユニタリ行列。
特異値分解をするには、Ｒには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ケースコントロールデータを６つの数値のリストとして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表の総サンプルサイズがＮとすると、特異値分解の $U,V\{*}$ は変わらず、 $D$ がＮ倍になる
$V^{*}$ は２ｘ２のユニタリ行列であるが、これは、ユニタリである条件(構成する２つの転置ベクトル、 $v_1{*},v_2{*}$ のそれぞれのノルムが１で、その内積がゼロであるという条件)から、 $|v_{11}v_{22}-v_{12}v_{21}|=1$ という性質がある。
ケースとコントロールとで分布に差がないとき、 $D$ の対角成分２つのうちいずれかがゼロである
観測テーブルの行列と期待度数の行列との差の行列の第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}$ で表される