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

ｋ次元球のk-1次元面上での平等なk方向 by RY

多次元球 立体角
  • こちらからのつづき
  • ユークリッド空間の場合
    • ユークリッド空間で正規直交基底を取る(\mathbf{e_1}=(1,0,0,...),\mathbf{e_2}=(0,1,0,...)のようなベクトルで張られたもの）
    • このベクトルの先端を結んだ多次元多角形は\sum_{i=1}^k a_i \mathbf{e_i}, \sum_{i=1}^k a_i=1で表される面上にある
    • この面はディリクレ面
    • 多次元多角形の重心をこの面上の原点に取り直すと、k個の頂点座標は、k-1次元空間中の正単体の頂点に相当する（こちらを参照）
    • これは、k-1次元ユークリッド空間にk本のベクトルを平等に配置する方法である
  • 球面座標の場合
    • 3次元の球面座標のオイラー表示を思いだそう（こちら）
      • x-y平面を三角関数で表して、そこにz軸を積み上げる形をしている
      • これを一般化した多次元球の角座標形もあってそれは、この本の冒頭などにも書かれている
    • ユークリッド空間のときにディリクレ面を取って、正規直交基底の先端が作る多角形の重心を原点として取り直したのと同じことを球面座標でもやるとする
    • そのときの重心は(\frac{1}{\sqrt{k}},\frac{1}{\sqrt{k}},...,\frac{1}{\sqrt{k}})なる点である
    • この原点から見て、k個の正規直交基底の先端はk-1次元亜空間中のk個の平等な点（正単体の頂点）方向にあるとみなせる
    • また、この頂点からの移動は、必ず、このk-1亜空間中の単位ベクトル（自由度k-2：k-1個の要素があり、その間に１つの制約式(ノルムが１)であるからk-2)とその球面上の距離（角度・弧の長さ）として表せる
    • 今、球面上の任意の点と、この原点(\frac{1}{\sqrt{k}},\frac{1}{\sqrt{k}},...,\frac{1}{\sqrt{k}})とを通る大円が作る、原点での接線方向と、この大円上の弧の距離とで、この点の座標を定めることとする
    • この大円と、正規直交基底に相当する正単体の頂点方向との関係を角度で表すことにする
    • この方法では、k個の正規直交基底のベクトルには順序や特別扱い（オイラー極座標でのz軸のような）はなくなっている
    • 平等化ができた
    • 大円方向と距離が定まると、そのときに、ある基底の頂点方向の座標がいくつになるかは計算可能である
    • この大円と基底頂点方向との関係（角度）によって、この基底の頂点方向の座標を定める関係が、『三角関数』の多次元版なのかもしれない
      • ２次元の場合には、方向は正か負か、しかなく、ある方向への距離は「角度」そのもの。その角度を使って、x,yの座標を定めているのがいわゆる\cos, \sinで表した三角関数
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])
}
RandomSphere <-
function (df = 3, r = 1, n = 100) 
{
    rs <- matrix(rnorm(df * n), nrow = n)
    rs/sqrt(apply(rs^2, 1, sum)) * r
}

k<-20 # 次元

# k-1の空間にk本のベクトルを均等に配置する(正単体)
Xs<-CategoryVector(k)
# それに、第k次元要素を0として与える
Xs2<-cbind(Xs,rep(0,k))

# k-次元空間の原点に、第k番目方向の成分が0の平面を取る
# (0,0,...,1)を中心とする単位球が、この平面に接しているとして
# 原点から、この平面上(k-1次元面)のベクトルxの方向への移動を始め、
# 球に沿って、長さtだけ移動したときの座標を返す関数
CoordsSphere<-function(k,x,t){
 O<-rep(0,k)
 O[k]<-1
 B<-tan(t)*x
#A<-(1-cos(t))*O+cos(t)*B
 Bp<-B-O
 Ap<-cos(t)*Bp
 Ap
}

# 原点を中心とする単位球を移動して、(0,0,...,1)を中心とするようにし、
# 原点を中心としたときに(1/sqrt(k),1/sqrt(k),...,1/sqrt(k))という座標にある点が原点に来るようにしたものとする
# このとき、単位ベクトル(1,0,0,...,0),(0,1,0,...,0),...,(0,0,...,0,1)の座標(の１例)をXspに格納する
Xsp<-Xs2
evec<-rep(0,k)
evec[1]<-1
cent<-rep(1/sqrt(k),k)
tunit<-acos(sum(evec*cent))
for(i in 1:k){
 Xsp[i,]<-CoordsSphere(k,Xs2[i,],tunit)
}

# k 個の正規直交基底ベクトルは単位ベクトルであって
# 相互に直交するしていることを確認する
r2<-matrix(0,k,k)
for(i in 1:k){
 for(j in 1:k){
  r2[i,j]<-sum(Xsp[i,]*Xsp[j,])
 }
}
r2


# 今、任意の方向に任意の距離(角度)で移動した点の座標をこの方法を用いて計算してみる
x<-c(RandomSphere(df=k-1,r=1,n=1))

x<-c(x,0)
t<-runif(1)*2*pi
# t<-tunit # 移動距離(移動角度)を正規直交基底のベクトルまでの距離にすると、relationとrelationVとは『ほぼ』一致する????
#t<-1.54
Ap<-CoordsSphere(k,x,t)

# この点と正規直交基底のベクトルとの関係を調べて見る
# relationは角度で表した関係
# euclidは正規直交基底での座標になる
relation<-euclid<-rep(0,k)
for(i in 1:k){
 euclid[i]<-sum(Ap*Xsp[i,])
 relation[i]<-acos(euclid[i])

}

relation
euclid
# また、Apの移動方向と、正規直交基底ベクトルの関係もみてみる
relationV<-rep(0,k)

for(i in 1:k){
 relationV[i]<-acos(sum(x*Xs2[i,]))
}

relationV

# relationV と relation って、順序が守られている
plot(relationV,relation)

#######
# 逆に、ユークリッド座標が与えられているときに、この点を「重心」からの方向と
# 角度で表してみる

#euclid<-rnorm(k)
#euclid<-euclid/sqrt(sum(euclid^2))



# Xsp %*% Ap = euclidなので
# Ap = Xsp^(-1) %*% euclid
Ysp<-solve(Xsp)
Ap2<-Ysp %*% euclid

# Ap と Ap2 とが同じことを確認
Ap-Ap2


# CoordsSphereの逆関数
CoordsSphereInv<-function(k,Ap){
 O<-rep(0,k)
 O[k]<--1
 t<-acos(sum(Ap * O))
 Bp<-1/cos(t)*Ap
 B<-O-Bp
 list(x=B/sqrt(sum(B^2)),t=t)
}

xANDt<-CoordsSphereInv(k,Ap)

# きちんと元に戻ることを確認
x2<-xANDt$x
t2<-xANDt$t
x-x2
t-t2


plot(euclid,relationV)
euclid/relationV
euclid/relation



# Function to estimate polynomial coefficients and coordinates estimated
# 多項式近似の係数と、推定座標を返す関数
glmPolynomial<-function(y,x,dg){# y:dependent values (従属値列),x:descriptive values (説明変数),dg:degrees of polynomials (多項式の次数)
 # Xs: Matrix consisted of x^1,x^2,...x^{dg}
 # Xs: xの1,2,...,dg 乗の値でできた行列
 Xs<-matrix(0,length(x),dg) 
 for(i in 1:dg){
  Xs[,i]<-x^i
 }
 # glm() is generalized linear regression function in R
 # glm() はＲの一般化線形回帰関数
 fit<-glm(y~Xs[,1:dg]) 
 # beta: coefficients estimated
 # beta: 推定係数ベクトル
 beta<-fit$coefficients
 # Xs2 has an additional column for x^0
 # Xs2 には、x^0のための列を追加
 Xs2<-cbind(rep(1,length(x)),Xs)
 # predY is a matrix of values for individual degrees of x
 # predY は xの各次数の項
 predY<-t(t(Xs2)*beta)
 # apply(): Sum up all columns corresponding to 0 to dg degrees
 # apply()関数で0:dg次数のすべてを足し合わせる
 sumupPredY<-apply(predY,1,sum)
 list(beta=beta,x=x,predY=sumupPredY)
}

# Generate data
x<-euclid
y<-relationV

xlim<-c(min(x),max(x))
ylim<-c(min(y),max(y))
plot(x,y,xlim=xlim,ylim=ylim,col="red")

# GLM is to be applied to multiple degrees (dgs)
# GLMを複数の次数に適用
dgs<-1:5
for(i in dgs){
 outglm<-glmPolynomial(y,x,i)
 par(new=TRUE)
 # gray() changes contrast 
 # gray() 関数は線の濃淡を変える
 plot(outglm$x,outglm$predY,xlim=xlim,ylim=ylim,col=gray(i/(max(dgs)*1.5)),type="l")
}