ryamadaの遺伝学・遺伝統計学メモ このページをアンテナに追加 RSSフィード

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

2017-08-11 多項観測ベクトルの異同

[][]多項観測ベクトルの異同

  • (たとえば)スキンフローラを観測して、2サンプル間でその異同を数値にしたいとする
  • Bray Curtis法というのがある。0-1におさまる
  • こちら
  • こちらの論文で使っている

2017-08-07 ぱらぱらめくる『Information Geometry and Population Genetics』

[][][]ぱらぱらめくる『Information Geometry and Population Genetics』

  • この本は、基本的に、集団のハプロタイプ頻度の変化をWright-Fisherモデルの下で扱うにあたって、多項分布を単体に対応付け、そこでの確率質量ベクトルの時間変化を多項分布の作る情報幾何空間での表現で書き下すことをしている
  • そういう意味では、ほとんどがいわゆる集団遺伝学の記載で、3 Geometric Structures and Information Geometryで情報幾何が導入されているので、そこの情報幾何と集団遺伝学との記法の対応性を確認すればよい
  • また、多項分布が単体に対応するところを、平方根を使って(多次元)球面に対応付けたり、そこでの酔歩を論じることも大事な要素
  • また、有限・離散から連続極限へと話を飛ばすこともしている
  • したがって、(1)集団遺伝学における多型の扱い、と、(2)多項分布の情報幾何表現(指数分布族表現)、との2つがわかれば、読める(し、読み終えたも同然)
  • 目次
  • 1 Introduction
  • 2 The Wright-Fisher Model
  • 3 Geometric Structures and Information Geometry
  • 4 Continuous Approximations
  • 5 Recombination
  • 6 Moment Generating and Free Energy Functions
  • 7 Large Deviation Theory
  • 8 The Forward Equation
  • 9 The Backward Equation
  • 10 Applications

2017-08-03 ぱらぱらめくる『言葉と物ー人文科学の考古学ー』

[][][]ぱらぱらめくる『言葉と物ー人文科学の考古学ー』

言葉と物―人文科学の考古学

言葉と物―人文科学の考古学

  • まずはあとがきから
    • キリスト教・マルクス主義を超える思想が第二次大戦後に様々な領域から出てきた。レヴィストロースとかもそう
    • 構造主義とまとめよう
    • それらの共通項の説明になっている
    • そのような解釈が可能(比較的容易)であったのは、西洋の諸学問領域が同根でそれらを説明する用語のルーツが同じであることから、語彙上の関係を言葉ハンドリングの上で行うことで、著者とその理解者(読者?)とに共通認識ができやすかったから(?)
      • 日本語は、西洋諸学問を輸入して、それぞれの学問領域で用語を新造・借用したので、その点での共感がしにくい
  • 目次
    • 第一部
      • 第一章 侍女たち (言葉以前?)
      • 第二章 世界という散文 (韻文ではなくて散文。叫び・心情としての韻文に対する、人文科学としての散文??)
      • 第三章 表象すること (音節記号としての言語(の部品))
      • 第四章 語ること (意味・論理)
      • 第五章 分類すること (離散記述統計)
      • 第六章 交換すること (集団としての人行動)
    • 第二部
      • 第七章 表象の限界
      • 第八章 労働、生命、言語(ランガール)
      • 第九章 人間とその分身
      • 第十章 人文諸科学
  • 各章の( )内コメントは、3分で読んで書いたコメント…
  • 上記にないのは、「連続・量的・微分・力学・量子・確率」
  • 「視覚が主」で「聴覚が従」で、その聴覚体系としての言語、ということなのか?

2017-07-28 L1ノルムで中央値、L2ノルムで平均値

[][][][][][][]L1ノルムで中央値、L2ノルムで平均値

  • 標本平均値は、標本の値の和を標本の数で割ることでえらる
  • それは、標本から推定された期待値
  • それは、L2ノルムの和を最小にするような値
  • 標本の中央値は、ソートして真ん中に来る値
  • それは、L1ノルムの和を最小にするような値
  • こちらも参考

f:id:ryamada22:20170727094348j:image

x <- c(rnorm(100,0,1),rnorm(200,3,1))
plot(x)
hist(x)

mean(x)
median(x)

X <- seq(from=min(x),to=max(x),length=1000)

L1 <- L2 <- rep(0,length(X))

for(i in 1:length(X)){
	L1[i] <- sum(abs(x-X[i]))
	L2[i] <- sum((x-X[i])^2)
}

matplot(X,cbind(L1,L2),type="l")

abline(v=mean(x),col=3)
abline(v=median(x),col=4)
abline(h=min(L1),col=3)
abline(h=min(L2),col=4)

2017-07-27 ぱらぱらめくる『Nonparametric Inference on Manifolds: With Appli

[][][][]ぱらぱらめくる『Nonparametric Inference on Manifolds: With Applications to Shape Spaces』

  • pdf
  • Preface
    • パラメトリックな方法があった
      • Landmarksbased shape spaces
      • Directional statistics
    • この本はノンパラメトリック
    • 大きな分類:Extrinsic vs. Intrinsic
      • Extrinsic
        • 多様体を高次元ユークリッド空間に埋め込んで、ユークリッド座標系を使う
        • パラメトリック
        • ノンパラメトリック
      • Intrinsic
        • 多様体を埋め込まずに扱う
    • Frechet mean
      • 多様体の上に確率分布があったときに、それのFrechet meanを検討するのが主題
      • 平均を扱うには、距離の定義が必要で、Extrinsicに(埋め込んだ高次空間のEuclid距離を使う)やるか、多様体上の距離(測地線・測地距離を使う)を使ってIntrinsicにやるかに分かれる
    • 構成
      • Chapter 1 Examples
        • 簡単な例
      • Chapter 2 Location and Spread on Metric Spaces
        • Frechet mean and dispersionのnotation
      • Chapter 3 Extrinsic Analysis on Manifolds
        • Extrinsic なFrechet mean and dispersionの推定
      • Chapter 4 Intrinsic Analysis on Manifolds
        • Intrinsicにやる。Chapter 4と対比的
      • Chapter 5 Landmark Based Shape Spaces
        • Landmark-based shape manifold
        • Shape spaces概要
      • Chapter 6 Kendall's (Direct) Similarity Shape Spaces ¥Sigma_m^k
        • Kendall's similarity shape spaces
      • Chapter 7
        • Planar shape spaces
      • Chapter 8
        • Reflection (similarity) shape spaces
      • Chapter 9
        • Stiefel manifolds
      • Chapter 10
        • Affine shape spaces
      • Chapter 11
        • Projective shape space
      • Chapter 12
        • ノンパラベイズでの分布推定(を使った形解析)の難しさ
      • Chapter 13
        • 形データで何をするか
        • 回帰、分類、検定

[][][][]1 Examples ぱらぱらめくる『Nonparametric Inference on Manifolds: With Applications to Shape Spaces』

  • 円周上の観測点セットから、分布の中心(平均)の点推定・区間推定
  • 球面上の観測点セットから、分布の中心(平均)の点推定・区間推定
  • ゴリラの頭蓋形の男女差比較・検定、分類器
    • 特徴点の位置の平均を求めるのに、Extrinsicにもできるし、Intrinsicにもできる
  • 健病の脳の形比較で検定・推定。ただしサンプル数が特徴点の数に比べて少ない場合
  • 手書き文字判定:Affine shape
  • 緑内障判定:神経乳頭の形解析

[][][][]2 Location and Spread on Metric Spaces ぱらぱらめくる『Nonparametric Inference on Manifolds: With Applications to Shape Spaces』

  • 空間があってそこに確率分布があるとき、その平均(期待値)とばらつきとを扱うのがこの本
  • ある点があったときに、その点と分布に従う点からのコスト(これの定義の仕方が色々ある)の期待値が最小になるとき、それを「平均」とする、これが基本
    • コストとしては「二乗距離」をとることが多いが、距離についても取り方があって、Extrinsic か Intrinsic かはこの本での2大別になる
    • 二乗距離にすると、mean、一乗距離にすると、median
  • 多様体の上に滑らかなコスト関数(と確率分布と)が定義されるべき

[][][][]3 Extrinsic Analysis on Manifolds ぱらぱらめくる『Nonparametric Inference on Manifolds: With Applications to Shape Spaces』

  • 多様体を高次元ユークリッド空間に埋め込む
  • 多様体上の2点の「二乗距離」を、埋め込んだユークリッド空間でのユークリッド二乗距離にしてしまう、というのがExtrinsicなアプローチの基本
  • 式はたくさん出てくるけれど、普通の推定の話
  • 高次元ユークリッド空間にどう埋め込むか、という問題はある
  • 多次元球面という多様体をExtrinsicに扱うという問題は、基本問題の1つ
  • 球をもう少し拡張するとStiefel manifoldになる
    • Stiefel manifoldが重要なのは、ある特定の長方形行列が作る多様体で、行列の期待値やばらつきの推定に使えるから
  • 地球の磁極推定
  • 地球上の火山の分布

[][][][]4 Intrinsic Analysis on Manifolds ぱらぱらめくる『Nonparametric Inference on Manifolds: With Applications to Shape Spaces』

  • 多様体上の距離を測地線に沿ってリーマン計量にのっとって測ることにするのが、すなおな距離のとり方
  • 二乗ノルムをコストにしてFrechetを定めるのが基本。それ以外のコストノルム関数も使える
  • このようにしたとき、IntrinsicなFrechet mean & dispersionが推定される
  • 用語・道具立て
    • Geodesic 測地線
    • Exponential map
      • 測地線が決まったら、測地線方向に進む初速度が接空間内ベクトルとして取れる。決まる。単位時間後の到達地点をexp_p(v) = ¥gamma(1)と定める。ただしpは多様体上の点、¥gamma(t)は測地線、vは初速度
    • Cut Locus
      • ある点pから、ある測地線上に到達するのに、かかる時間は単位速度で進むと決めることで定まる。測地線以外の経路でもその点に到達することができるが、測地線で行くのが一番短時間で行ける場合がある(そのような場合も多い)。測地線経路が最短時間経路であるような、測地線上の点であって、かつ、もっとも時間がかかる点の集合をCut Locusと言う(球面の場合だったら、測地線を逆に辿って到達することとの競争になるので、「ぐるり周回測地線」における「最遠点」がそれになる。球面なら対蹠点
    • Sectional Curvature
      • 曲率を表すリーマン測度はテンソル。1つのスカラーな曲率指標としてガウス曲率がある。それはその点における全方位の曲率のスカラー化だが、それを、四分扇について値化したもの
    • Injectivity Radius
      • 多様体上の点から、そのCut Locusまでの測地線距離を考える。その値を多様体上の点すべてについて考慮したとき、その最小値が定まる。それのこと
    • Convex Set
      • 多様体上のある部分(set)がConvexであるとは、そのsetの任意の2点間に唯一の最短な測地線が存在すること
  • 上記、道具立てを使って、多様体上のある点における観測を多様体上のどこの点にどのくらいの重みとするべきか、というのを定める、重みの計算には、「距離」に応じたコスト関数が用いられる。

[][][][]5 Landmark Based Shape Spaces ぱらぱらめくる『Nonparametric Inference on Manifolds: With Applications to Shape Spaces』

  • m次元空間(多様体)にk > m 個の点があるとき
  • 変換不変量を問題にすることで、k個の点座標の異同を定量して、検定・推定・分類等をする
  • 形空間(Shape spaces)の点として形を扱う
    • Kendall's similarity shape spaces (平行移動と回転を無視。スケール変換で標準化)
    • Reflection shape spaces (直交変換(回転と反転)を無視)
    • Affine shapes (Affine transformationsを無視)
    • Projective shapes (Projective transformationsを無視)
  • Shape spacesは、ユークリッド空間の部分空間ではない
  • Shape spacesはリーマン多様体Nを変換群Gを法とした剰余群[M=N/G])
    • NがS^dであることも多い
  • 変換は多様体上の軌道を作る(群作用がfreeなとき)
  • 点には接空間があって、その接空間には軌道とそうでないものとに分かれる。軌道とそうでないものは相互に直交
  • Extrinsic analysisでは、Nをユークリッド空間に埋め込むから、Nも埋め込まれ、軌道も埋め込まれる
  • Intrinsic analysisでは、軌道同一視によって、NをMに沈めこむ(Riemannian submersion)
    • 沈めこんだ減次元多様体では、リーマン測度(内積)の変換も出来る(内積が保たれる?)
    • これを使って、沈めこんだ後の多様体での推定等を、オリジナルの多様体での距離・測度を使って考えることが出来る
  • いくつかのShape spaces
    • (Direct) similarity Shape Spaces ¥Sigma_m^k
      • m=2(2次元)の場合には、群作用はfree(単位元のみが、不変作用を持つ)なので、変換同一視の空間¥Sigma_2^kはコンパクト微分可能多様体である
    • Reflection Similarity Shape Spaces R¥Sigma_m^k
      • m>2のときsimilarity transformationsの群作用はfreeとは限らない。それは¥Sigma_m^kが必ずしも多様体にならないことを意味する
      • Reflectionも同一視すると、m>2でのDirct similarity Shape Spacesの問題が解消する(こともある)ので、これを使ってExtrinsic解析をする手もある
    • Affine Shape Spaces A¥Sigma_m^k
      • 二次元泳動像の試行ごとのずれをAffine変換とみな巣ことで、この枠組みに入れることが出来るらしい
    • Projective Shape Spaces P¥Sigma_m^k
      • コンピュータビジョンでは透視図法的なことをやるが、それに使えるのがProjective shape spaces

[][][][]6 Kendall's (Direct) Similarity Shape Spaces ¥Sigma_m^k ぱらぱらめくる『Nonparametric Inference on Manifolds: With Applications to Shape Spaces』

  • d次元球面上の点集合があったときに、m=2,3次元回転を同一視して、商群を取ったもの
  • S^d/G
  • G = SO(m): 特殊直交群 (mxm直交行列であって、行列式が +1)
  • m=2,3が大事な例
  • 特徴点がm=2次元座標として得られているとする
  • 特徴点の数k=30の例
  • まず、m=2次元の両軸について原点中心にする(自由度が2減る)
  • さらに、mxk次元ベクトルとして、「長さ」を1に標準化する
  • 自由度はmkから(2+1)下がってmk-3

f:id:ryamada22:20170727163717j:image

m <- 2
k <- 30

z <- matrix(rnorm(m*k),nrow=m)
z[1,] <- z[1,]*0.3 + 1
z[2,] <- z[2,]*0.2 + 2

ang <- Arg(z[1,] + i * z[2,])
z <- z[,order(ang)]


plot(t(z),pch=20,asp=1)

z.normed <- z - apply(z,1,mean)
plot(t(z.normed),pch=20,asp=1)

z.normed. <- z.normed/sqrt(sum(z.normed^2))

plot(t(cbind(z,z.normed,z.normed.)),pch=20,asp=1,type="p")
points(t(z),pch=20,col=1,type="b")
points(t(z.normed),pch=20,col=2,type="b")
points(t(z.normed.),pch=20,col=3,type="b"
  • 自由度mxkをmxk- m - 1に下げた。そのうち、最後の「ノルム=1」の自由度の減は、線型変換ではないので、以下に示すように、線型変換によってmxk-m-1次元平面に映すことは出来ないが、mxk-m次元平面にある、単位球に移すことは出来る
    • 以下では、多数のm=2,k=3の長さmxk=6のベクトルを作成し、それを重心標準化、ノルム標準化することによって得られる、6次元空間座標を持った多様体上の点が、mxk-m=4次元平坦空間に存在し、それがノルム=1の制約を持っていること、すなわち、mxk-m次元にあるS^{mxk-m-1}という多様体を成していることを示すRのコードである
library(MASS)
library(Matrix)
library(rgl)

n <- 1000 # データセットの数
# そのmxn次元座標を納める行列
Z <- matrix(0,n,m*k)

for(i in 1:n){
	z <- matrix(rnorm(m*k),nrow=m)
	z[1,] <- z[1,]*0.3 + 1
	z[2,] <- z[2,]*0.2 + 2
	z.normed <- z - apply(z,1,mean)
	z.normed. <- z.normed/sqrt(sum(z.normed^2))
	Z[i,] <- z.normed.
}

# ランクがmxk-m
# mxk-m次元空間に収まっていることがわかる
rankMatrix(Z)
# QR分解を用いて、mxk-m次元空間に移す回転行列を計算する
# 発生させたmxkベクトルのうち、ランク本のベクトルを用いる
df <- m*k - (m)
topZ <- Z[1:df,]
# QR分解
qr.out <- qr(t(topZ))
Q <- qr.Q(qr.out)
R <- qr.R(qr.out)

# QR分解により、t(topZ)のペアワイズなベクトル内積が
# R(減次元ベクトルの束)のペアワイズなベクトル内積と一致することを確認する
# 以下はessentially zero
range(t(R) %*% R - topZ %*% t(topZ))

# t(topZ) = QR と出来たから
# Q.inv t(topZ) = R のQ.invを作る
Q.inv <- ginv(Q)
# 検算
round(Q.inv %*% t(topZ) - R)
# Q.invを使って、すべてのmxkベクトルを減次元写像する
rotZ <- Q.inv %*% t(Z)
# ペアワイズ内積が保たれていることを確認する
d1 <- Z %*% t(Z)
d2 <- t(rotZ) %*% rotZ
range(d1-d2)

# 減次元座標のノルムの確認
range(apply(rotZ^2,2,sum))

# mxk -m = 4次元空間の球面はなので、3次元空間にプロットすると
# 球面と球の内部に点が分布する
plot3d(t(rotZ[1:3,]))
plot3d(t(rotZ[2:4,]))
> rankMatrix(Z)
[1] 4
attr(,"method")
[1] "tolNorm2"
attr(,"useGrad")
[1] FALSE
attr(,"tol")
[1] 2.220446e-13
> # QR分解により、t(topZ)のペアワイズなベクトル内積が
> # R(減次元ベクトルの束)のペアワイズなベクトル内積と一致することを確認する
> # 以下はessentially zero
> range(t(R) %*% R - topZ %*% t(topZ))
[1] -1.110223e-16  1.110223e-16
> # 検算
> round(Q.inv %*% t(topZ) - R)
     [,1] [,2] [,3] [,4]
[1,]    0    0    0    0
[2,]    0    0    0    0
[3,]    0    0    0    0
[4,]    0    0    0    0
> range(d1-d2)
[1] -1.221245e-15  1.221245e-15
> 
> # 減次元座標のノルムの確認
> range(apply(rotZ^2,2,sum))
[1] 1 1

f:id:ryamada22:20170727173223p:image

  • このS^{m¥times k - m -1}球面は、Preshape sphereと呼ぶ。なぜなら、まだ回転同一視をしていないから
  • m=2の場合、回転は群 SO(2)であってそれはS^1(円環)なので、S^{m¥times k - m- 1}/G = S^{2k-3}/S^1なる、2k-4次元のコンパクトな多様体であることがわかる
  • これはk-2次元のcomplex projective spaceと(微分)同相
  • これは次のようにも書ける
    • S_m^k = ¥{z ¥in R^{m¥times k} | Trace(z t(z)) = 1, z 1_k =0¥}
      • 行列z、ただし、Traceの制約は、zの全要素の二乗の和が1であること、1_kなる、全成分が1の長さkのベクトルを含む式は、k要素のすべての平均が0でるように標準化した制約に相当する
    • うまく等長変換して、mxk-m次元平面に移せるので
      • S_m^k = ¥{z ¥in R^{m ¥times k -m} | Trace(z t(z)) =1¥}とまで持ち込める
      • こうすることで、線型制約が無くなるので便利
      • 本ではHelmert変換を使って-m次元する、と書いている
      • 上のコードでは、QR分解を使っている
  • さて、回転同一視をする
  • Direct Similarity Shape Spaces
    • ¥Sigma_m^k = S_m^k/SO(m)
    • m>2では¥Sigma_m^kは多様体になるとは限らない
      • m¥times k行列のランクがm-1以上なら、多様体になる
    • そのような場合をNS_m^k = ¥{z ¥in S_m^k|rank(z) ¥ge m-1¥}と書けば
    • ¥Sigma_{0m}^k = NS_m^k/SO(m)で、これは多様体であって、その次元はkm-m-1-¥frac{m(m-1)}{2}だという
    • 微分同相が担保されれば、接空間、接空間を軌道とその直交成分に分けることなどの議論が進められる
    • ※ 行列と転置行列の積のトレースは、行列をベクトル化して、その内積を取るのと同じであることに注意
    • 接空間はT_z S_m^k = ¥{v ¥in H_m^k| Trace (v t(z))=0¥}、ただしH_m^kS_m^kを含んだ超平面。その超平面内の球上の接空間は球面上のベクトルzとの直交条件で示せる
  • 回転同一視をして、その上での形間(測地線)距離はmin_{T ¥in SO(m)} d_{gs}(x,Ty)、ただしd_{gs}(x,y) = arccos(Trace(yt(x))と表される、球面上の測地線
  • したがってarccos(max_{T ¥in SO(m)} Trace(Ty, t(x)))
  • Trace((Ty) x^T)の最大化はy x^Tのsvd分解によって求まる
    • 最大値y x^T = U¥Lambda Vというsvd分解でのTrace(¥Lambda) = ¥sum ¥lambda_iであり
    • それをもたらす回転はT= U^T V^T
library(svd)

# x,y,z座標球面スカラー場に対して
# k個のスペクトルが得られたとする
# それが形情報
# x,y なるm * k 行列

# 3次元ベクトルがk個
m <- 3
k <- 5

# x,y,z座標を平均0にして
# 全要素の二乗和を1に標準化する
x <- matrix(rnorm(m*k),ncol=k)
y <- matrix(rnorm(m*k),ncol=k)

x <- t(t(x)-apply(x,2,mean))
y <- t(t(y)-apply(y,2,mean))
x <- x/sqrt(sum(x^2))
y <- y/sqrt(sum(y^2))

# 検算
apply(x,1,sum)
apply(y,1,sum)
sum(x^2)
sum(y^2)

# svd分解法
# http://www2.stat.duke.edu/~ab216/monograph.pdf の
# p84 の式(6.4)の辺りを使って
# yを回転して、xと近づける最適化を行う

yxt <- y %*% t(x)

svdout <- svd(yxt)
U <- svdout$u
V <- svdout$v

# svd 分解の確認
# テキストの式表現とsvd()関数の式表現のずれに注意する
U %*% diag(svdout$d) %*% t(V)
yxt

# テキストで言うところの最適回転行列 TをWとする
W <- t(t(V)) %*% t(U)

# Trace(Ty t(x)) = sum(eigenvalue)の検算

sum(diag(W %*% y %*% t(x)))
sum(svdout$d)
W %*% y

[][][][]7 The Planar Shape Space ¥Sigma_2^k ぱらぱらめくる『Nonparametric Inference on Manifolds: With Applications to Shape Spaces』

  • Direct Similarity Shape Spacesの中でm=2の場合
  • 平面上の形の解析は形解析の基本であり、需要も多いので、取り立てる
  • 二次元平面上の点を複素平面上の点とみなす
  • 複素数重心で標準化すると、重心を中心としたk次元複素ベクトルが形に対応する
  • これを¥lambda e^{i¥theta}倍することは、伸縮して回転することに相当するので、これによって一致するベクトルは同一視することにする
  • 方向同一視で複素数的に1次元減り、回転同一視で複素数的にさらに1次元減る
  • これは複素射影空間 CP^{k-2}となる
  • k個の二次元座標を長さkの複素ベクトルzとし、u = ¥frac{z-¥bar{z}}{||z-¥bar{z}||}が標準化複素ベクトル。これを回転すると¥pi (u) = ¥{e^{i ¥theta} u | -¥pi ¥le u ¥le ¥pi¥}。この曲線が「形」を表す
  • これは射影平面H^{k-1}=¥{w ¥in C^{k}/¥{0¥}|¥sum_{i=1}^k w_i ~ 0¥}にあるCS^{k-1}複素球面上の点(CS^{k-1} = ¥{u ¥in C^{k}|¥sum_{i=1}^k u_i=0,||u||=1¥}
  • 形間距離はRe(e^{-i¥theta} Conj(w) z)の最小値をコサインとする角度
  • この回転はe^{i¥theta} = ¥frac{Conj(w)^T u}{|Conj(w)^T u|}で与えられるという

f:id:ryamada22:20170729060901j:image

m <- 2
k <- 5

x1 <- matrix(rnorm(m*k),ncol=m)
x2 <- matrix(rnorm(m*k),ncol=m)

x1. <- x1[,1] + 1i * x1[,2]
x2. <- x2[,1] + 1i * x2[,2]

u1 <- (x1.-mean(x1.))/sum(Mod(x1.-mean(x1.)))
u2 <- (x2.-mean(x2.))/sum(Mod(x2.-mean(x2.)))

sum(u1)
sum(Mod(u1))


z <- sum(Conj(u1) * u2)/Mod(sum(Conj(u1) * u2))

thetas <- seq(from=-1,to=1,length=1000)*pi

d <- rep(0,length(thetas))
for(i in 1:length(thetas)){
	tmp <- thetas[i]
	tmp2 <- u1 * exp(1i*tmp)
	d[i] <- abs(Re(sum(Conj(tmp2) * u2)))
}

plot(thetas, d,type="l")

abline(v=Arg(z),col=2)
  • Intrinsic Analysis in ¥Sigma_2^k
    • ¥Sigma_2^k上の点サンプルが複数あったときに、そこにどんな分布があるとみなすか、という解析ができる
      • 球形に配置されているとみなせるときは、方向統計学の意味での、中心と正規分布の分散共分散行列との推定ができる(これはFrechet分布)
      • より一般的な意味でのFrechet分布も推定できる
  • Extrinsinc Analysis in ¥Sigma_2^k
    • all k × k complex Hermitian matricesにマップできる(らしい)
      • Veronese-Whitney embeddingと呼ぶらしい
      • 複素空間単位球面上の複素ベクトルuがあったときに、u_i ¥bar{u_j}(i,j)要素とする行列のこと
    • このようしにして、Hermitian matricesがたくさん得られるので、それの「平均」とか「分布」とかを推定することができる
    • 分布にできれば、検定もできる

[][][][]8 Reflection (Similarity) Shape Spaces R¥Sigma_m^k ぱらぱらめくる『Nonparametric Inference on Manifolds: With Applications to Shape Spaces』

  • m次元k個の座標xを、preshape S_m^k上のzに変換し、¥sigma(z) = ¥{Az | A ¥in O(m)¥}というS_m^kの集合をorbitと呼ぶ。ただしO(m)は直交変換
  • このときR¥Sigma_m^k = ¥{¥sigma(z) ¥in S_m^k, rank(z) = m¥} = NS_m^k/O(m)という同相関係にある
    • ただしS_m^k = ¥{z ¥in R^{m¥times k} | Trace(z z^T) = 1,z 1_k = 0¥}で、NS_m^kはそのnon-singular な部分
    • これは、¥Sigma_m^kのうちのnon-singularな部分(これは回転同一視)であって、さらにReflectionについてだけ同一視条件を付加したものになる(R¥Sigma_m^k=¥Sigma_{0m}^k/G,GはReflection
  • 左右眼球の神経乳頭の形とかは、Reflection同一視の対象

[][][][]9 Stiefel manifolds V_{k,m} ぱらぱらめくる『Nonparametric Inference on Manifolds: With Applications to Shape Spaces』

  • この章がここに入っているのはなぜなんだろう???
  • (Direct) Similarity Shape Spaces ¥Sigma_m^kと、その変形であるRelfection Shape SpacesはR^{m¥times k}行列であって、その成分の2乗和が1という制約を基本としている
    • S_m^k = ¥{z ¥in R^{m ¥times k} | Trace(z z^T) = 1,z 1_k = 0¥}, m ¥le k
    • これに回転同一視をして、商にしたり、線型独立性制約をつけたり、するバリエーションがある
  • 他方、Stiefel manifoldは以下に示すように、行列の形は同じだが、条件は違う
    • テキストとm,kの取り方が異なるが、以下のように書くことにする
    • V_m^k = ¥{z ¥in R^{m ¥times k} | z^T z = I_k¥}
  • 本来はm’ ¥le m次元空間にあるm個の点が、k次元座標を持っているとき、Stiefel manifoldの行列を使うと、m次元空間座標に等長変換できる
    • 以下に示す方法で変換行列と変換後座標を算出できる
      • A=QR (QR分解)。Q ¥in V_m^k
    • m’==mのときは、以下の方法でもできる
      • A=UP;U= A(A^T A)^(-1/2),P=(A^T A)^(1/2)
library(GPArotation)
library(Matrix)

m <- 3
k <- 5

my.rstiefel <- function(m,k){
	tmp <- Random.Start(k)
	return(tmp[,1:m])
}

x <- my.rstiefel(m,k)

round(t(x) %*% x, 5)

# ランクをm未満にする
real.dim <- m-1

# real.dimランクのm次元ベクトルをk個作る
mu <- matrix(rnorm(real.dim*k),ncol=real.dim)
mu <- cbind(mu,matrix(0,k,m-real.dim))
tmp.rot <- Random.Start(m)
mu <- mu %*% tmp.rot


rankMatrix(mu)

U <- mu %*% solve(sqrtm((t(mu) %*% mu)))
P <- sqrtm(t(mu) %*% mu)

qr.out <- qr(mu)
U. <- qr.Q(qr.out)
P. <- qr.R(qr.out)

# 分解はちゃんと出来ている
round(U %*% P - mu,5)
round(U. %*% P. - mu,5)

# 両者は違う
U - U.
P - P.
# 等長変換になっていない
round(t(P) %*% P - t(mu) %*% mu,5)
round(t(P.) %*% P. - t(mu) %*% mu,5)
# m'=real.dim < mなので、UはStiefel manifoldに乗っていない

round(t(U) %*% U,5)
round(t(U.) %*% U.,5)


# ランクをmにする
real.dim <- m

# real.dimランクのm次元ベクトルをk個作る
mu <- matrix(rnorm(real.dim*k),ncol=real.dim)
mu <- cbind(mu,matrix(0,k,m-real.dim))
tmp.rot <- Random.Start(m)
mu <- mu %*% tmp.rot


rankMatrix(mu)

U <- mu %*% solve(sqrtm((t(mu) %*% mu)))
P <- sqrtm(t(mu) %*% mu)

qr.out <- qr(mu)
U. <- qr.Q(qr.out)
P. <- qr.R(qr.out)

# 分解はちゃんと出来ている
round(U %*% P - mu,5)
round(U. %*% P. - mu,5)
# 両者は違う

U - U.
P - P.
# 等長変換になっている
round(t(P) %*% P - t(mu) %*% mu,5)
round(t(P.) %*% P. - t(mu) %*% mu,5)
# m'=real.dim == mなので、UはStiefel manifoldに乗る

round(t(U) %*% U,5)
round(t(U.) %*% U.,5)

[][][][]12 Nonparametric Bayes Inference on Manifolds ぱらぱらめくる『Nonparametric Inference on Manifolds: With Applications to Shape Spaces』

  • unit sphereとplanar shape space について記述する
  • ノンパラなのでDirichlet processを使ったり
  • ベイズなので事前分布、尤度、KLdivergenceを使ったりする