次のようなデータを考えてみます。
```{r}
n <- 1000
t <- seq(from=-1,to=1,length=n) * 2 * pi * 10
d <- 6
x <- 1.1^t * cos(t)
y <- 1.1^t * sin(t)
z <- 1.1^t
X <- cbind(x,y,z,matrix(0,n,d-3)) + rnorm(n*d,0,sd(c(x,y,z))*0.15)
library(GPArotation)
R <- Random.Start(d)
X. <- X%*%R
plot(as.data.frame(X.),asp=1,cex=0.1)
```
６行６列で色々なプロットが並んでいます。
これは、1000個のサンプルについて、6種類の検査結果があるときに、1000個のサンプルがどんな具合なのかを絵にしてみたものです。
6種類の検査があるので「6次元空間」上に1000個の点が散らばっており、その散らばり具合が知りたいのです。
6次元空間は次元が高いので、視覚的に捉えるためには、6種類の検査から2種類ずつをペアにして、そのすべてをプロットしています。
$$
\begin{equation}
_6 C_2 = \frac{6!}{2!4!} = 15
\end{equation}
$$
全部で15ペアできます。それらを縦軸・横軸を入れ替えた対称なプロットも含めることで30通りのプロットになっています。
なんとなく形がわかりますが、もう少しすっきりとわかる工夫はないでしょうか。
これまでに何度か行列の利用に触れてきましたが、こういうときに、行列は「この6次元のオブジェクトをうまい具合に回転して、見やすい軸を取り出す」という計算を素早く確実にやってくれるので、多次元データの解析では行列を頻繁に利用します。
その結果が以下の図になります。
```{r}
plot(as.data.frame(X),asp=1,cex=0.1)
```
オオマカニ言って、左上の3つの軸に集中して、「らせん」「円錐」といった「意味のある形」があります。
それ以外のパネルには、ランダムに広がった線と円とがあるだけで、「意味の少ない広がり」が現れています。
この主要な3軸だけを取り出して残りを捨ててしまっても、それなりにデータの全体を把握していることになります。
こんな波形を描く現象があるとします。
```{r}
t <- seq(from=0,to=1,length=1000)*2*pi*4
d <- 4
rs <- 1:4
ks <- c(1,sqrt(2),exp(1),pi)
phis <- 1:4*pi/5
x <- rep(0,length(t))
for(i in 1:d){
x <- x + rs[i]*cos(ks[i]*t+phis[i])
}
x. <- x + rnorm(length(t),0,sd(x)*2)
ylim <- range(c(x,x.))
plot(t,x,type="l",col=2,ylim=ylim)
```
しかしながら、実際の観測データとしては、次のような雑音が入ったデータが得られてしまうことはよくあることです。
```{r}
plot(t,x.,type="l",ylim=ylim)
points(t,x,type="l",col=2)
```
こんなとき、この雑音だらけの観測データは、時間軸に沿った観測値としては雑音だらけですが、時間に対して、見る立場を変えてみます。
周期を$1,\frac{1}{2},\frac{1}{4},...$のように半分にしながらだんだんに細かい三角関数の波を考えて、それが足し合わさってこの波が見えているものと想定します。
このようにして異なる周期の波の足し合わせに分解してやって、どの周期の波の成分が多いのかを判定し、成分の多い波と、成分が少ない波とにより分けます。
フーリエ変換という方法です。
さきほど、6次元空間のオブジェクトの見る方向を変えてやり、意味のある方向はとっておき、意味のなさそうな方向は捨てたのと同じ発想です。
このようにして、上記6つの周期の波だけを取り出してやると、次のような波が得られます。
```{r}
fft.out <- fft(x.)
mod.fft <- Mod(fft.out)
ord <- order(mod.fft)
fft.out. <- fft.out
fft.out.[ord[1:(length(t)-6)]] <- 0
x.. <- fft(fft.out.,inverse=TRUE)/length(t)
plot(t,x..,type="l",col=3,ylim=ylim)
```
滑らかなカーブが得られました。
これが、雑音に隠されていた本当の波とどういう関係かを示すと次のようになります。
```{r}
plot(t,x,type="l",col=2,ylim=ylim)
points(t,Re(x..),type="l",col=3)
```
ぴったり一致はしませんが、雑音が邪魔していたことを考えれば随分と当てはまりのよい推定ができたと言えるでしょう。
この技術は雑音の入る情報処理の基本ルーチンの一つですから、医療画像などの分野では避けて通れない知識です。
このように、データはそのまま見ていてもなかなか背後に隠れているルールを見つけることは難しいのですが、見る方向を変えると、集約されることは多いです。
それは、現象を形造っているルールは比較的単純なのですが、そこに乱雑な項が入りこんでわかりにくくしている上に、私たちが普通にアプローチするときのものの見方は、背景のルールを見出しにくい方向であることが重なっているからかもしれません。
そんなときに数学は視点変更という強力な武器を提供します。
行列を使った視点変換、フーリエ変換と呼ばれる周期への視点変換を紹介しました。
どちらも強力なツールですが、計算機で一発計算できる、ある意味で単純な仕組みですし、「こうだったら、こうなるべし」という、予測しやすい現象が対象となっていました。
次は、予期せぬことが起きるときのアプローチです。
「予期しないことが起きる」というのはどういうことなのでしょうか。
「Xだったら、Yになる」と思っていたところが、「YではなくZにな」り、「ZがYと全然違う」ときに、びっくりします。
もし「Xだったら、Yになると思っていたけれど、Zになってびっくりした。でも、2回目もXだったらZになったので、もうびっくりしない」というのなら、これはもう「予想範囲内」です。
「XだったらZになってびっくりしたけれど、2回目にXだったら、YでもZでもなくWになってびっくりした。何が起こるか予想がつかない」というのが、「予期せぬこと」でしょう。
このような予期せぬことには大きく2つの原因があると考えることが多いです。
1つは、本当にランダムに起きている場合。サイコロの目の予想がつかないようなこと。
もう一つは、Xが正確にわかれば、Y,Z,Wのどれになるかは決まっているのだけれど、XのちょっとのずれがY,Z,Wへの道を分けるので、予想しがたい、という場合。
後者は結果が、初期値に敏感な系と呼びます。
そんな例を絵で見てみましょう。
以下の例では、3要素が簡単な連立微分方程式の関係を作っています。微分方程式(相互作用のルール)は簡単ですが、3要素の量の変化を3次元の状態空間で描くと突然値が変わったりしてびっくりします。
$$
\begin{equation}
\frac{dx}{dt} = -10x + 10y\\
\frac{dy}{dt} = 28x - y - xz
\frac{dz}{dt} 0 \frac{8}{3}z + xy
\end{equation}
$$
```{r,warning=FALSE}
Niter<-100000
dt<-0.001
xs<-matrix(0,Niter,3)
xs[1,]<-c(0.3,0.4,0.5)
for(i in 2:Niter){
xs[i,1]<-xs[i-1,1]+(-10*xs[i-1,1]+10*xs[i-1,2])*dt
xs[i,2]<-xs[i-1,2]+(28*xs[i-1,1]-xs[i-1,2]-xs[i-1,1]*xs[i-1,3])*dt
xs[i,3]<-xs[i-1,3]+(-8/3*xs[i-1,3]+xs[i-1,1]*xs[i-1,2])*dt
}
```
```{r,echo=FALSE,warning=FALSE}
library(rgl)
```
```{r setup,warning=FALSE}
knit_hooks$set(rgl = function(before, options, envir) {
if (before && rgl.cur() > 0) rgl.close()
hook_rgl(before, options, envir)
})
```
3要素の時間変化のグラフと、それを3次元空間での軌跡としたものとを示します。
```{r}
obs <- 1:1000
matplot(obs,xs[obs,],type="l")
```
```{r, rgl=TRUE,echo=FALSE}
plot3d(xs[obs,],col=rainbow(1000),cex=0.1)
```
これは比較的、早期の様子ですが、しばらく時間が経った後に、再度観察すると次のようになります。
落ち着いた状況に見えます。
```{r}
obs <- 10001:11000
matplot(obs,xs[obs,],type="l")
```
```{r, rgl=TRUE,echo=FALSE}
plot3d(xs[obs,],col=rainbow(1000),cex=0.1)
```
落ち着いているので、しばらく目を離していて、再度、観察を始めたところ
```{r}
obs <- 70001:80000
matplot(obs,xs[obs,],type="l")
```
```{r, rgl=TRUE,echo=FALSE}
plot3d(xs[obs,],col=rainbow(1000),cex=0.1)
```
途中で目を離さずにずっと観察すると次のようになります。
```{r}
obs <- 1:length(xs[,1])
matplot(obs,xs[obs,],type="l")
```
```{r, rgl=TRUE,echo=FALSE}
plot3d(xs[obs,],col=rainbow(1000),cex=0.1)
```
3次元状態空間に2つの渦巻きがあって、そのどちらかをぐるぐる回っているかと思えば、2つの円を行ったり来たりする局面も登場します。
この類の動きは、よくよくその背景がわかれば現象に不思議なことはないのですが、「XだったらYとなるだろう」という判断や推定をするときに、「直線的〜線形」な関係があるものとして考え始める習性がある人間は、このような現象を「不思議だ、カオスだ」と考えます。
しかしながら、このような変化は複数の要素が集まるとごく普通に発生する現象ですし、生命現象の魅力はこのような一筋縄ではいかないものであることも多いので、非線形と呼ばれるこのような現象を掴まえるための基礎を持つことも有用でしょう。
腫瘍という増えなくてもよいところで細胞がもりもりと増えるタイプの病気があります。
悪性の腫瘍を癌と呼んだりもします。
そんな腫瘍ですが、クルリとしたきれいなボールのような場合もあれば、レタスのようにひだひだがたくさんある場合もあります。その中間はひらべったい腫瘍でしょうか。
クルリと丸い腫瘍とレタスのような場合とでは、総じてレタス型の方がタチが悪いことになっています。
では、こんな形状の違いはどのように説明すればよいでしょうか。
細胞増殖の速度の違い(速い方がタチが悪いことが多いでしょう)にその説明を求めることもできます。
幾何学的に「球状」と「レタス状」と「平坦状」との違いは、曲率として説明できます。
平坦なのは、曲率がゼロです。
「球状」なのは、常に曲がろう、曲がろうとしていて、その曲がろうという傾向を「正の曲率」と呼び、球(や円)は曲率が正であって、どこでも一定な図形のことです。
曲率がゼロの場合と正の場合がありましたが、曲率が負の場合というのもあるのでしょうか。
それが、「レタス状」の場合です。
曲率が正の場合は「面積」が不足するので、そのつじつまを合わせるために「閉じて」しまいますが、
曲率が負の場合は「面積」が常に余ってくるので、2次元空間に納まりきらずひらひら・びらびらと3次元空間にはみ出してきます。
レタス状の広がりは双曲平面と呼ばれますが、描いてみればこんな具合です。
```{r}
library(igraph)
x <- c(1)
k <- 11
current <- 1
n.kid <- 1:2
P <- 0.3
p <- c(P,1-P)
el <- matrix(0,0,2)
for(i in 1:k){
n <- sample(n.kid,length(current),replace=TRUE,prob=p)
new.kid <- sum(n)
cumsum.kid <- cumsum(n)
cumsum.kid <- c(0,cumsum.kid)
max.x <- max(x)
if(new.kid > 1){
tmp <- (max.x+1):(max.x+new.kid-1)
el <- rbind(el,cbind(tmp,tmp+1))
}
for(j in 1:length(n)){
el <- rbind(el,cbind(rep(current[j],n[j]),(max.x + cumsum.kid[j]+1):(max.x+cumsum.kid[j+1])))
}
current <- (max.x+1):(max.x+new.kid)
x <- c(x,current)
}
g <- graph.edgelist(el,directed=FALSE)
```
```{r, rgl=TRUE,echo=FALSE}
coords <- layout.fruchterman.reingold(g, dim=3)
rglplot(g, layout=coords,vertex.size=0.1,vertex.label=NA)
```
医学・生命科学は生き物の仕組みを解明し理解するために、どんな道具を使っても構いません。
最近の医学・生命科学は数学・情報学・計算機科学の最先端ですが、医学科のカリキュラムでは、残念ながら、医学生命科学のための数学を提供できていません。また、一般教養科目としての数学も、医学・生命科学のための数学とはなっていません。
ですが、医学・生命科学への数学的なアプローチは非常に多様で奥が深く、今後、ますます重要性を増すと思われます。
本日、紹介しきれなかった内容で、本講義の教員が興味を持っている内容が、断片的に、
> http://www.genome.med.kyoto-u.ac.jp/StatGenet/lectures/MyBook/mybookcurrent.svg
> http://www.genome.med.kyoto-u.ac.jp/wiki_tokyo/index.php/%E3%83%A2%E3%82%B8%E3%83%A5%E3%83%BC%E3%83%AB%E3%82%92%E4%BD%BF%E3%81%A3%E3%81%A6%E5%AD%A6%E3%81%B6
にあります。
興味を持った学生さんは覗いてみてください。
