Hatena::ブログ(Diary)

魚料理好きの日記

2016-09-01

 1.2 複雑な変化

周期関数の振幅をα倍、周期数をβ倍していき、その級数和をとると、ワイエルシュトラス関数とよばれ、
x(t)=¥sum^{N}_{j=0} ¥alpha^{j+1} ¥cos (¥beta ^j t)  , | ¥alpha | < 1
と表せる。
ホワイトノイズはこれで表現できる。
x( ¥beta t)= ¥alpha^{-1} x(t)- ¥cos(t)
で、時間軸と関数値軸を異なるスケール変換しているので、自己アフィン相似性とよぶ。

s<-0.1
T<-1:100*s
N<-50
alpha<-2/3
beta<-2
x<-c()
ax<-c()
#x(t)の計算
for(i in 1:length(T)){f<-function(x,i){(alpha^(x+1))*cos((beta^(x))*i*s)}
                      x[i]<-sum(f(0:N,i))
                      }
plot(T,x,type="l",ylim=c(-1.5,1.5),ylab=(""))
#αx(βt)の計算
for(i in 1:length(T)){f<-function(x,i){(alpha^(x+1))*cos((beta^(x+1))*i*s)}
                      ax[i]<-sum(f(0:N,i))
                      }
par(new=TRUE)
plot(T,ax,type="l",col=2,ylim=c(-1.5,1.5),ylab=(""))
par(new=TRUE)
#2関数の差
plot(T,(alpha*ax-x),type="l",col=3,ylim=c(-1.5,1.5),ylab=(""))

f:id:AOdragon:20150313184926j:image

べき関数x(t)=A(t)t^{¥gamma}¥alpha x ( ¥beta t)=x(t)を満たすときの条件は
¥alpha ¥beta ^{¥gamma}=1  ,  A(¥beta t)=A(t)なので
x(t)=t^{¥frac{¥log 1/¥alpha }{¥log ¥beta}} ¥sum^{¥infty}_{j=-¥infty} A_j  ¥exp (i 2 ¥pi ¥frac{¥log t}{¥log ¥beta} j )と表される。
グラフとしてのフラクタル次元Dは
D=2-¥frac{¥log ( 1/¥alpha )}{¥log ¥beta}
と表される。
H=2-Dはハースト数とよばれる。

べき関数x(t)=¥sum^{¥infty}_{j=-¥infty} A_j t^{H_j}と書き直すと
べき指数H_j=¥frac{¥log 1/¥alpha }{¥log ¥beta}-¥frac{i2¥pi j}{¥log ¥beta}複素数で、人の肺・地震金融に応用されるらしい。

 1.1 複雑な形

複雑系と呼ばれる事象数学的にどのように認識されるのかについての本。
複雑系の数理

長さ尺度を例に。
単位尺度rを用いて、求められる長さをL(r)とする。
このとき、L(r)はrについて減少関数
マンデルブロ海岸線について
L(r)¥approx r^{1-D}
が成立することを発見した。
現象がべき関数であらわせるとき、べき則に従うといい、成り立つrの領域をスケーリング領域という。
上式のDをフラクタル次元という。
L(r)=rN(r)
とおくと、N(r)は尺度で測った回数となる。
今、rを¥frac{r}{2}としてみると
L(¥frac{r}{2})=cL(r)...(c ¥geq 1)
となる。c=1のときは直線などの滑らかな直線のとき成り立つ。
上式は関数方程式なので、べき関数が解となる。
すると、
D=1+¥frac{¥log c}{¥log 2}
と表せる。

c>1(D>1)のとき、フラクタル図形という。
1<D<2では平均値は存在しない(発散する)。
べき則が成り立つ場合、rについてスケール変換不変性が成り立つ。

多次元に拡張すると、
今単位容積に対して、尺度rを用いるとN(r)回充てる必要があったとすると、
N(r)=(¥frac{1}{r}^{D})
なので、
D=¥frac{¥log N(r)}{¥log (¥frac{1}{r})}
と表せる。
それぞれの次元で異なる場合
D=¥frac{¥log N_{x}(r)}{¥log (1/r_{x})}+¥frac{¥log N_{y}(r)}{¥log (1/r_{y})}+...
と表せる。
さらに、非整数となるフラクタル次元を計算するために
D=¥lim_{r ¥to 0}¥frac{¥log N(r)}{¥log (¥frac{1}{r})}
と拡張して定義される。

具体例として
カントール集合ではD=¥frac{¥log 2}{¥log 3}=0.63093...となる。
シェルピンスキー・ガスケットではd次元のとき、D=¥frac{¥log (d+1)}{¥log 2}となる。

2015-12-12

4章 Evolutionary games

4.2 Nash equilibrium
2種類の生き物(もしくは戦略)A,Bがいて、AがそれぞれA,Bと出会った時の利得をa,bとし、Bの場合c,dとする。
これを↓のようなpayoff行列で表すとする。
¥matrix{ & A & B ¥cr A & a & b ¥cr B & c & d}

このとき、
1.Aが狭義のNash均衡⇔a>c
2.AがNash均衡⇔a≧c
3.Bが狭義のNash均衡⇔d>b
4.BがNash均衡⇔d≧b

4.3 Evolutionary stable strategy(ESS)
Aが多数いる中で変異種Bがはびこらないような条件
AがESS⇔a>c or (a=c and b>d)

4.4 More than 2 strategies
strategyS_{j}に対するstrategyS_{i}のpayoffをE(S_{i},S_{j})とする。

1.strategyS_{k}が狭義のNash均衡⇔E(S_{k},S_{k})>E(S_{i},S_{k}), ¥forall i ¥neq k
2.strategyS_{k}がNash均衡⇔E(S_{k},S_{k}) ¥succeq E(S_{i},S_{k}), ¥forall i ¥neq k
3.strategyS_{k}ESSE(S_{k},S_{k})>E(S_{i},S_{k}), ¥forall i ¥neq kもしくはE(S_{k},S_{k})=E(S_{i},S_{k}), E(S_{k},S_{i})>E(S_{i},S_{i}),  ¥forall i ¥neq k
4.strategyS_{k}が広義のESSE(S_{k},S_{k})>E(S_{i},S_{k}), ¥forall i ¥neq kもしくはE(S_{k},S_{k})=E(S_{i},S_{k}), E(S_{k},S_{i}) ¥succeq E(S_{i},S_{i}),  ¥forall i ¥neq k
条件の強さとして1>3>4>2

無敵なstrategyはE(S_{k},S_{k})>E(S_{i},S_{k}),E(S_{k},S_{i})>E(S_{i},S_{i})で定義される

4.5 Replicator dynamics
戦略i(もしくは変異種)に対する戦略jのpayoffをa_{ij}とするとn×nのpayoff行列Aができる。
また、fitnessをf_{i}(¥vec{x})=¥sum^{n}_{j=1}a_{ij}x_{j}とすれば
再生産の方程式¥dot{x_{i}}=x_{i}¥[ f_{i}-¥phi¥] として表せる。

4.6 Hawk or Dove?
タカ派ハト派をモチーフにしたゲーム。
¥matrix{ & H & D ¥cr H & ¥frac{b-c}{2} & b ¥cr D & 0 & ¥frac{b}{2} }
自分と敵のタカ派ハト派をとる確率をp_{1},p_{2}とすると
期待所得E(p_{1},p_{2})E(p_{1},p_{2})=¥frac{b}{2}(1+p_{1}-p_{2}-¥frac{c}{b}p_{1}p_{2})で表される。
安定点はp^{¥ast}=¥frac{b}{c}のとき。

4.7 Thre is always a Nash equilibrium
自分も相手もn通りの戦略をとれるとし、それぞれの戦略をとる自分の分布を¥vec{p}=(p_{1},p_{2},...,p_{n})、相手の分布を¥vec{q}=(q_{1},q_{2},...,q_{n})とする。
payoff matrixを A={a_{ij}} とすれば、期待所得E(¥vec{p},¥vec{q})=¥sum^{n}_{i=1}¥sum^{n}_{j=1}a_{ij}p_{i}q_{j}=¥vec{p}A¥vec{q}と表せる。
このとき、¥vec{p}A¥vec{p}¥geq¥vec{q}A¥vec{p} , ¥forall ¥vec{q}なる¥vec{p}が必ず存在する。
ペロン―フロベニウスの定理使えば証明できそう(たぶん)。

4.9 Game theory and ecology
ロトカ・ヴォルテラの式について

a<-1
b<-1
c<-1
d<-1
f<-function(x,y){
x<-x+x*(a-b*y)
y<-y+y*(-c+d*x)
return(c(x,y))
}
n<-50
LV<-matrix(0,2,n+1)
LV[,1]<-c(1.2,1.2)
for(i in 1:n){
LV[,i+1]<-f(LV[1,i],LV[2,i])
}
matplot(t(LV),type="l")

f:id:AOdragon:20151212113232j:image

2015-09-09

3章 Fitness landscape and Sequence spaces

3.1 sequence space
長さLのsequenceがあったときに、有り得るsequenceは4^{L}通り。(A,T,C,G)
このL次元空間をsequence spaceとする。
これを2進数(0,1)で考えると、L次元の超立方体の格子点があり得る組み合わせとなる。
各格子点の距離はマンハッタン距離で考えられる。

3.2 Fitness landscape
Fitness landscapeはsequence space状の格子点に対し、数値(reproduction rate)を与えたもの。
各格子点の評価みたいなもの。

3.3 the equasispecies equation
2進数でのsequence spaceを考え、それぞれの格子点を10進数に変換すると0から2^{L}-1まで1対1で配列番号として与える。
i番目の配列に対して、x_{i},f_{i},Q_{ji}をそれぞれiの割合、適応度、配列jからiへの変異行列とすると、
¥dot{x_{i}}=¥sum^{n}_{j=1}x_{j}f_{j}Q_{ji}-¥phi x_{i}
ただし、¥phi = ¥sum^{n}_{i=1} f_{i} x_{i}
をequasispecies equationという。
Qが単位行列、つまり変異が起きないなら最も高いf_{i} をもつ配列iが生き残る。
Qが単位行列でないなら、単体内のある点x^{¥ast} に収束する。このときφは必ずしも最大値とはならない。
しかし、総個体数はφの指数関数で増加する。

#試行回数
t<-30
#変異しうる配列数
n<-10
#遺伝子iを持つ個体の割合ベクトル
x<-runif(n)
x<-x/sum(x)
#遺伝子iを持つ個体のfittnessベクトル
f<-runif(n)
#mutation matrix
Q<-matrix(runif(n*n),n,n)
q<-apply(Q,1,sum)
Q<-Q/q
#average fittness
phi<-x%*%f

a<-matrix(0,n,t+1)
a[,1]<-x
for(i in 1:t){
dx<-(x*f)%*%Q-phi%*%x
x<-x+dx
a[,i+1]<-x
phi<-x%*%f
}
matplot(t(a),type="l")

f:id:AOdragon:20150814012054j:image

3.4 a mutation matrix for point mutations
配列i,jのハミング距離をh_{ij}として、点変異の入る確率をuとすると
変異行列Qの各成分(iからjと変異する確率)は
q_{ij}=u^{h_{ij}}(1-u)^{L-h_{ij}}

3.5 a adaptation is localization in sequence space
変異する確率が高すぎると適応していた個体(fittnessの値が高い)の割合がいまいち増えない。

f:id:AOdragon:20150831234254j:image
f:id:AOdragon:20150831234253j:image
上の方は変異確率を1/100にしたもので、下の方は1/10にしたもので、上だともっとも高いfittnessをもつものに収束する。
単純なモデルでこの変異の閾値を計算するとu<¥frac{1}{L} という条件が求められる。
つまりuL<1で、RNAウイルス大腸菌酵母菌、マウス、人などでuLの実測値が1より低いという論文が出ているらしい(Drake(1991,1998))。
一方で、QβウイルスとVZVウイルスは1よりはるかに大きく、その理由は不明。

3.6 selection of the quasispecies
2峰性のfittness landscapeをもち、1つは高いが狭く、1つは高くはないが幅広とする。
このとき、変異の割合によって生き残りが変わってくる。
変異の割合が高すぎればすべての種が生き残り、低すぎると最も高いものに収束する。
その中間であれば、2つのピークをもつ種が生き残る。

#試行回数
t<-100
#変異しうる配列数
n<-30
#遺伝子iを持つ個体の割合ベクトル
x<-runif(n)
x<-x/sum(x)
#遺伝子iを持つ個体のfittnessベクトル
f<-c(1,rep(0,n-1))+c(rep(0,n-5),0.1,0.2,0.3,0.2,0.1)
#mutationの割合を調整するパラメーター
h<-100
#mutation matrix
Q<-matrix(1,n,n)+h*diag(n)
q<-apply(Q,1,sum)
Q<-Q/q
#average fittness
phi<-x%*%f

a<-matrix(0,n,t+1)
a[,1]<-x
for(i in 1:t){
dx<-(x*f)%*%Q-phi%*%x
x<-x+dx
a[,i+1]<-x
phi<-x%*%f
}
plot(a[,t+1])

f:id:AOdragon:20150909222353j:image:w640
f:id:AOdragon:20150909222352j:image:w640
f:id:AOdragon:20150909222351j:image:w640
f:id:AOdragon:20150909222350j:image:w640

2015-08-31

2章 what evolution is

evolutionary dynamicsの基本概念replication,selection,mutationについて。
replication:自己の(遺伝情報も含めて)複製。
selection:異なる個体間での淘汰。
mutation:異なる(遺伝情報をもつ)個体の生産。

1.reproduction
期間毎にそのときの個体数x(t)とある正の定数rの積で増殖するとき、
¥dot{x}=rx
と表され、解は
x(t)=x_{0}e^{rt}
となる。
これと同時に、ある正の定数dで死亡するとき、
¥dot{x}=(r-d)x
と表され、寿命の平均は¥frac{1}{d}となり、個体の増加率はrとなる。
¥frac{r}{d}>1のときに、個体数は増加していく。

最大値Kの制限付き増加モデルは、上式に制限項をかけたもので
¥dot{x}=rx(1-¥frac{x}{K})
と表され、解は
x(t)=¥frac{Kx_{0}e^{rt}}{K+x_{0}(e^{rt}-1)}
と表される。
この解はKに収束する。

今、K=1として差分方程式を考える。
x_{t+1}=ax_{t}(1-x_{t})
負の値を取らない条件は0 ¥leq a ¥leq 4
収束x^{¥ast}はaの値によって変わる。
0 ¥leq a <1のとき、x^{¥ast}=0
1<a<3のとき、x^{¥ast}=¥frac{a-1}{a}
3<a ¥leq 4のとき、x^{¥ast}は不安定。

f<-function(a,x){y<-a*x*(1-x)
                 return(y)}
t<-50
a<-1:40*0.1
x<-0.5
X<-matrix(0,t,length(a))
X[1,]<-x
for(i in 2:t){
X[i,]<-f(a,X[i-1,])
}
matplot(X[,1:10],type='l')
matplot(X[,11:30],type='l')
matplot(X[,31:40],type='l')

f:id:AOdragon:20150831213351j:image
f:id:AOdragon:20150831213352j:image
f:id:AOdragon:20150831213353j:image

このモデルはさらに決定論カオスで初期値鋭敏性をもつ。

f<-function(a,x){y<-a*x*(1-x)
                 return(y)}
t<-50
a<-4
x<-c(0.3156,0.3157)
X<-matrix(0,t,length(x))
X[1,]<-x
for(i in 2:t){
X[i,]<-f(a,X[i-1,])
}
matplot(X[,1:2],type='l')

f:id:AOdragon:20150831213354j:image:初期値0.3156と0.3157の挙動の違い]

2.selection
2種類の生物の個体数x(t),y(t)がそれぞれa,bの割合で初期値増殖するモデル
¥dot{x}=ax,¥dot{y}=by
について
¥rho (t)=x(t)/y(t)
とすることで、a,bの大きさで淘汰される側がわかる。

総個体数一定の条件x+y=1の下では
¥dot{x}=x(1-x)(a-b)
が得られ、均衡はx^{¥ast}=0,1
また、a,bの大きさによって収束先が変わる。
これが適者生存の概念である。
多次元に応用すれば
¥phi = ¥sum{x_{i} f_{i}}かつ¥sum{x_{i}}=1として
¥dot{x_{i}}=x_{i}(f_{i}-¥phi)とすると、
単体上で考えられる。
fが最も大きい個体に収束する。

a<-matrix(0,500,3)
a[1,1]<-runif(1)/2
a[1,2]<-runif(1)/2
a[1,3]<-1-a[1,1]-a[1,2]
f1<-1/2
f2<-1/3
f3<-1/6
S<-function(x,y,z){sigma<-sum(x*f1+y*f2+z*f3)
                    return(c(x+x*(f1-sigma),y+y*(f2-sigma),z+z*(f3-sigma)))}
for(i in 1:499){a[i+1,]<-S(a[i,1],a[i,2],a[i,3])}
library(rgl)
plot3d(a,type="l",col="2")

この3次元verでは初期状態(1/2,1/3,1/6)から(1,0,0)に収束する軌道になる。(f1が最も大きいので)

3.mutation
2タイプの個体X,Yの個体数x,yについてmutationを起こしてもう一方の個体を生む割合をそれぞれu_{1},u_{2}とすると
¥dot{x}=x(1-u_{1})+yu_{2}-¥phi x
¥dot{y}=x u_{1}+y(1-u_{2})-¥phi y
としてx,yの時間変化は表される。
x+y=1とすると、¥phi =1であり、
このときxの均衡値x^{¥ast}=¥frac{u_{2}}{u_{1}+u_{2}}となる。
3つ以上のタイプの場合、mutation matrix,Q=[q_{ij}]を用いて
¥dot{¥vec{x}}=¥vec{x}Q-¥phi ¥vec{x}と表される。
Qはn*nの推移確率行列となる。


4.mating
閉鎖的な環境下、変異も起きない中でramdom matingが行われればallele頻度は保存される(Hardy-Weinbergの法則)

1章 introduction

Nowakの書いた生物の進化の数学モデルの本。
asin:B00J97FFRI:Evolutionary dyanamics]

1章は進化生物学の歴史と各章の紹介。
・進化生物学のだいたいの歴史
ダーウィン:進化の概念を形にした人。
メンデル遺伝の公式の発見者。
ハーディ・ワインバーグ:その名を冠した(法医で聞いた)公式を作ったコンビ。
Fisher,Haldane,Wright:進化、選別、変異のモデル化。
木村資生:中立進化説を提唱。
ウィリアム・ハミルトン:血縁選択説を提唱。

2012-10-08

確率論 シリンダー集合の意義

シリンダー集合(Cとする)とはコイン投げの有限回の試行といったものの結果の全体集合である。
Cは無限回の試行の結果の集合の部分集合(Ωとする)である。
シリンダー集合の集まりC'とするとき、(Ω、C')上の有限加法的測度μがとれる。
Cから生成される加法族(σ(C')とする)というのはCを含む最小のσ加法族のことである。

これらを前提にすると、ホップの拡張定理というものにたどりつく。E.ホップの拡張定理 - Wikipedia
つまり、μがσ加法性を保有してれば、σ(C')上の測度として使えるし、それは一意的に定まる。

シリンダー集合は測度の拡張のための道具ととらえられると思う。