Md ばらばらに起きること

提供: StatGenetKyotoU
移動: 案内検索
ばらばらに起きること
カテゴリ 名称 説明
医学生物学トピック
数学トピック ベルヌーイ過程 起きるか起きないか
数学トピック ポアッソン過程 ぽつりぽつりと起きること
R関数・Rスキル 確率分布関係の諸関数 dxxxx(),pxxxx(),qxxxx(),rxxxx()

目次

このモジュールで登場する分布などの相互関係

  • この節でWikipedia記事へリンクが張ってある項目は、その項目を読むことが必要である。Wikipediaの記事のすべてが理解できないかもしれないが、それでも、通読することは必要である。Wikipediaの記事のうち、「確率分布に共通した内容」が、記事の右カラムの「確率分布の特徴一覧表」の内容になる。この「特徴一覧表」の項目は確率分布を相互に比較して理解するための基本項目である。それについては、md_確率密度関数・累積分布関数・ハザード関数・累積ハザード関数で扱うので、そちらも適宜参照すること。
  • ベルヌーイ過程(Wikiの記事参照)と二項分布
  • ポアッソン過程とポアッソン分布から、つながる分布
    • 1次元空間でのポアッソン事象
      • 生起間隔(最近接点間距離)の分布として、指数分布
      • i番近接点間距離の分布として、アーラン分布(i版のiを連続にしたのがガンマ分布
        • 複数回起きることで達成される現象は、アーラン分布・ガンマ分布を取る
    • 多次元空間のポアッソン事象
      • 最近接点間距離の分布として、ワイブル分布(したがって、ワイブル分布の1次元版が指数分布)
        • 時間経過とともに、生起確率が増加する(老朽化による故障頻度の上昇)などを説明する分布
        • 複数の要因のどれかが故障すると全体の故障につながる場合を説明する分布

ばらばらに起きることの整理整頓

どこで・いつ、ばらばらに起きるか

どのくらい起きるか

  • 起きる回数だけ、起きる回数の分布だけを考慮する場合には、「0次元時空間」での現象を観察している。もしくは、「0次元空間」において「一定時間」かけて起きる現象を観察し、その回数を集計して(時間のことは忘れて)いる。
  • 起きる回数のとらえ方による違いについては、#二項分布と正規分布とポアッソン分布を参照

起きていることの『何を』測定して、『何を』問題にするか

  • 時空間における「発生」を観測する
    • 時空間の次元
        • 1次元
          • 離散的な場合
          • 連続的な場合
        • 2以上次元
          • (離散的な場合)
          • 連続的な場合
  • 「発生点」の分布を問題にする
    • 「発生点」の数による大別
      • たくさん発生する場合
      • まばらに発生する場合
  • 問題にすること
    • たくさん発生する場合
      • 密度(単位時空間体積あたりの回数)を問題にする
    • まばらに発生する場合
      • 単位時空間体積あたりの回数(これが離散的)を問題にする
      • 発生間隔(発生しない時空間の広さの分布)を問題にする
      • 特定回数の発生点を持つ時空間の広さを問題にする
      • ポアッソン過程で何を問題にするかについては#ポアッソン過程で何を問題にするかを参照

ばらばらに起きることが集まって、初めて、観察できるときのこと

  • ばらばらに起きていること、それ自体を観察できることもあれば、そうでないときもある
  • そのようなときの「観察できること」の発生の様子は、「ばらばらに起きていることそれ自体」の発生の様子とは異なってくる
  • ばらばらに起きていることが、何個か集まって初めて、何かしらの事件が起き、その事件が観察できる場合もある
    • たとえば、複数の安全装置に守られている現象(癌化など)が該当する
    • 個々の安全装置の破綻がばらばらに起きることであり、すべての安全装置が破綻するまでは、何も見えないが、すべてが破綻すると癌が見えてくる
  • ばらばらに起きることが場所がいくつかあって、そのうちの一つでも、起きれば、何かしらの事件が起き、その事件が観察できる場合もある
    • 精密な機構があって、複数のステップの組み合わせによって達成されている現象(代謝サイクルなど)が該当する
    • 個々のステップの破綻がばらばらに起きることであり、どれか一つでも破綻すると、現象全体に異常が生じて、それが観察される
  • 「ばらばらに起きること」同士の組み合わせと関係とがあり、その結果として、初めて、観察できる事件が起きる場合については#ばらばらな出来事の組み合わせがもたらすこと

複数種類のばらばらな事象

  • 前の節でも、複数のばらばらなことがある場合を扱った
  • 前の節の場合には、ばらばらなこと自体は観察できず、ばらばらなことが集まって起こす現象を観察することができた。そして、その現象の起きる様子を問題にした
  • 本節では、2(多)種類のばらばらなことがあって、それらを観察できる場合が対象である
  • 例えば、ある現象がばらばらに起きていて、別の種類の現象もばらばらに起きているとする。ここで、ある現象ともう一つの種類の現象の発生が誘発関係にあるのではないか、と疑ったときには、「ある現象ともう一つの種類の現象との最近接距離」の分布が「ばらばらな仮定」より短めに偏っているかを知る必要がある。2種類のばらばらに起きる現象の因果関係を問題にするときに気になる分布である
  • これについては#2(多)種類のばらばらな出来事の相互関係を参照

ばらばらに起きること、と、ばらばらでなく起きることを比べる

  • 「ばらばらに」起きるとは、次の2つの要素からなる
    • どことも知れず、いつとも知れず、「平等に」起きる
    • 「ばらばらに」~「勝手気ままに」~「相互に無関係に」起きる
  • 2要素の組み合わせが作る4パターンを以下の表のようにパターン1,2,3,4とする
平等に 不平等に
相互に無関係に パターン1 パターン2
相互に関係を持って パターン3 パターン4
  • 以下、パターン1、パターン2、パターン3、パターン4の順にプロットを描いてみる。「平等」と「無関係」の意味するところを理解すること
  • パターン1
# パターン1 「平等に」「無関係に」
# 次元
k<-2
# 1軸の領域長(L^kの領域に「ポツリポツリ」を起こします)
L<-1
# 生起点の数
N<-1000
# 生起点の座標
# ポアッソン過程
# 縦軸・横軸、ともに均一分布からの乱数としてとる
X<-matrix(runif(k*N)*L,ncol=k)

plot(X,cex=0.01,pch=19)

Pattern1.jpeg

  • パターン2
# パターン2 「不平等に」「無関係に」
# 次元
k<-2
# 1軸の領域長(L^kの領域に「ポツリポツリ」を起こします)
L<-1
# 生起点の数
N<-1000
# 生起点の座標
# 正規分布で発生させ、LxL領域に納まるように、座標を調整する
X<-matrix(rnorm(k*N),ncol=k)
minX<-min(X)
maxX<-max(X)
X<-(X-minX)/(maxX-minX)*L
plot(X,cex=0.01,pch=19,xlim=c(0,L),ylim=c(0,L))

Pattern2.jpeg

  • パターン3-1
# パターン3-1 「平等に」「関係を持って」
# 次元
k<-2
# 1軸の領域長(L^kの領域に「ポツリポツリ」を起こします)
L<-1
# 生起点の数
N<-10000
g<-10
# 10分の1をパターン1で発生させ
# 残りの10分の9は、10分の1の点の周辺に集中するように発生させる
X<-matrix(runif(k*N/g)*L,ncol=k)
for(i in 1:g){
	tmpX<-X+matrix(rnorm(k*N/g,sd=0.01),ncol=k)
	plot(tmpX,cex=0.01,pch=19,xlim=c(0,L),ylim=c(0,L))
	par(new=TRUE)
}
par(new=FALSE)

Pattern3-1.jpeg

  • パターン3-2
# パターン3-2 「平等に」「関係を持って」
# 次元
k<-2
# 1軸の領域長(L^kの領域に「ポツリポツリ」を起こします)
L<-1
# 生起点の数
N<-1000
# パターン1で発生させ
# その後、各点を最近傍3点の重心になるように再配置
X<-matrix(runif(k*N)*L,ncol=k)
for(i in 1:N){
	D<-apply((t(t(X)-X[i,]))^2,1,sum)
	Neighbor3<-order(D)[2:4]
	X[i,]<-apply(X[Neighbor3,],2,sum)/3
	plot(X,cex=0.01,pch=19,xlim=c(0,L),ylim=c(0,L))

}
plot(X,cex=0.01,pch=19,xlim=c(0,L),ylim=c(0,L))

Pattern3-2.jpeg

  • パターン4-1
# パターン4-1 「不平等に」「無関係に」
# 次元
k<-2
# 1軸の領域長(L^kの領域に「ポツリポツリ」を起こします)
L<-1
# 生起点の数
N<-10000
g<-10
# ポアッソン過程
X<-matrix(rnorm(k*N/g),ncol=k)
minX<-min(X)
maxX<-max(X)
X<-(X-minX)/(maxX-minX)*L

for(i in 1:g){
	tmpX<-X+matrix(rnorm(k*N/g,sd=0.01),ncol=k)
	plot(tmpX,cex=0.01,pch=19,xlim=c(0,L),ylim=c(0,L))
	par(new=TRUE)
}
par(new=FALSE)

ファイル:Pattern3.jpeg

# パターン4-1 「不平等に」「無関係に」
# 次元
k<-2
# 1軸の領域長(L^kの領域に「ポツリポツリ」を起こします)
L<-1
# 生起点の数
N<-10000
g<-10
# 10分の1をパターン2で発生させ
X<-matrix(rnorm(k*N/g),ncol=k)
minX<-min(X)
maxX<-max(X)
X<-(X-minX)/(maxX-minX)*L
# 残りの10分の9は、パターン3-1と同様に集中発生させる
for(i in 1:g){
	tmpX<-X+matrix(rnorm(k*N/g,sd=0.01),ncol=k)
	plot(tmpX,cex=0.01,pch=19,xlim=c(0,L),ylim=c(0,L))
	par(new=TRUE)
}
par(new=FALSE)

Pattern4-1.jpeg

  • パターン4-2
# パターン4-2 「不平等に」「無関係に」
# 次元
k<-2
# 1軸の領域長(L^kの領域に「ポツリポツリ」を起こします)
L<-1
# 生起点の数
N<-1000
# パターン2で発生させ
N<-1000
# 生起点の座標
X<-matrix(rnorm(k*N),ncol=k)
minX<-min(X)
maxX<-max(X)
X<-(X-minX)/(maxX-minX)*L
# その後、パターン3-2と同様に、最近傍3点の重心に再配置
for(i in 1:N){
	D<-apply((t(t(X)-X[i,]))^2,1,sum)
	Neighbor3<-order(D)[2:4]
	X[i,]<-apply(X[Neighbor3,],2,sum)/3
	plot(X,cex=0.01,pch=19,xlim=c(0,L),ylim=c(0,L))

}
plot(X,cex=0.01,pch=19,xlim=c(0,L),ylim=c(0,L))

Pattern4-2.jpeg

二項分布と正規分布とポアッソン分布

起きる回数はまちまちだが、大きく3つに分けられる

  • 二項分布と正規分布とポアッソン分布
    • 数えられる程度にたくさん起きる(二項分布)
    • 数えきれないくらいたくさん起きる(正規分布)
      • 二項分布の極限
    • ぽつりぽつりと起きる(ポアッソン分布)
      • これも二項分布の極限
  • 二項分布
    • 時空間は1次元
    • 時空間は離散的
    • 時空間全体を観察する
  • 起きる回数は「単位時空間あたりの頻度」x「時空間の広さ」に依存する
  • 二項分布が正規分布とポアッソン分布に収束する様子を以下で確認する
    • 二項分布→正規分布
p<-runif(1) #生起確率
# 正規分布に近づける
# N[i]*pが二項分布の期待値
# N[i]*pを期待値とし、sqrt(p(1-p)N[i])を分散とする正規分布と
# N[i]*pを期待値とする二項分布はN[i]が大きくなるにつれてほぼ一致する様子を
# N[i]を1からmaxNへと増やしながら確認する
maxN<-100
N<-1:maxN # 試行数を増やしていく
for(i in N){
 dn<-dnorm(0:i,p*i,sqrt(p*(1-p)*i)) #正規分布の生起確率
 db<-dbinom(0:i,i,p) #二項分布の生起確率
 ylim<-c(0,max(dn,db))
 plot(dn,ylim=ylim,type="l")
 par(new=TRUE)
 plot(db,ylim=ylim,type="h",col=2)
}

NormBinom.jpeg

    • 二項分布→ポアッソン分布
      • 二項分布の確率
        • 1回ごとに、起きる確率がpとして、n回試行すると、平均してnp回起きると予想されるが、そのような状況でk回起きる確率
        • P_{{binom}}(k)={\begin{pmatrix}n\\k\end{pmatrix}}p^{k}(1-p)^{{n-k}}={\frac  {n!}{k!(n-k)!}}p^{k}(1-p)^{{n-k}}
        • kがk+1になると、確率はどれくらい小さくなるだろうか?
        • {\frac  {P_{{binom}}(k+1)}{P_{{binom}}(k)}}={\frac  {{\frac  {n!}{(k+1)!(n-(k+1))!}}p^{{k+1}}(1-p)^{{n-(k+1)}}}{{\frac  {n!}{k!(n-k)!}}p^{k}(1-p)^{{n-k}}}}={\frac  {n-k}{k+1}}{\frac  {p}{1-p}}
      • ポアッソン分布の確率
        • 1回ごとに、起きる確率がpとして、n回試行すると(ただしn\to \infty )、平均してnp回(ただし\lim _{{n\to \infty }}np=\lambda )起きると予想されるが、そのような状況でk回起きる確率
        • 今、k回起きる確率をP_{{pois}}(k)とすれば
        • {\frac  {P_{{pois}}(k+1)}{P_{{pois}}(k)}}=\lim _{{n\to \infty }}{\frac  {P_{{binom}}(k+1)}{P_{{binom}}(k)}}=\lim _{{n\to \infty }}{\frac  {n-k}{k+1}}{\frac  {p}{1-p}}=\lim _{{n\to \infty }}{\frac  {n-k}{k+1}}{\frac  {{\frac  {n}{\lambda }}}{1-{\frac  {n}{\lambda }}}}
        • {\frac  {P_{{pois}}(k+1)}{P_{{pois}}(k)}}=\lim _{{n\to \infty }}{\frac  {\lambda }{k+1}}{\frac  {n-k}{n-\lambda }}={\frac  {\lambda }{k+1}}
        • したがって、P_{{pois}}(k)=P_{{pois}}(0){\frac  {\lambda ^{k}}{k!}}
        • ここで\sum _{{i=0}}^{{\infty }}P_{{pois}}(i)=1を満足するようにP_{{pois}}(0)を取ると
        • P_{{pois}}(k)={\frac  {\lambda ^{k}}{k!}}e^{{-\lambda }}という
p<-runif(1)
m<-5
pm<-p*m # 期待度数
# ポアッソン分布に近づける
# 期待度数pmを固定して、試行回数N[i]を増やすにつれて生起確率pm/N[i]を小さくして行き
# ポアッソン分布の確率と二項分布の確率とが収束していく様子を観察する
maxN<-100
N<-ceiling(pm):maxN # 試行
for(i in N){
 xlim<-c(0,max(pm*10,30))
 po<-dpois(0:i,pm)
 bi<-dbinom(0:i,i,pm/i) # 生起確率は、期待度数を試行回数で除したもの
 ylim<-c(0,max(po,bi))
 plot(po,xlim=xlim,ylim=ylim,type="h")
 par(new=TRUE)
 plot(bi,xlim=xlim,ylim=ylim,type="l",col=2)
}

PoisBinom.jpeg

1次元と多次元のポアッソン過程

ポアッソン過程1,2,3、次元を眺めてみる

  • ポアッソン過程では、何かしらが、「ぽつり、ぽつり」と起きます
  • どこでも、同じ程度に起きます
  • ある起きること、と、別の起きることとは、相互に独立して起きます
  • 1次元
    • 放射性同位体の崩壊(放射性崩壊)を時間軸で観察すると1次元のポアッソン過程になります
    • 染色体交叉・組み換えも、「染色体が十分に長」く染色体1本当たり、1から数か所でしか起こらないので、DNA分子を1次元空間と見立てた、ポアッソン過程と見ることもできます。ただし、実際の交叉・組み換えは、起きやすい場所と起きにくい場所があること、2つの交叉が近接して起きることには、制約があることなどから、「平等」という点でも、「相互に無関係」という点でも、「ばらばら(ポアッソン過程)」の性質を満足しません

PoisDF=1.jpeg

# 前の節では、k=2 2次元平面での様子を見ましたが、1次元空間や3次元空間でもポアッソン過程は起こせます
# 次元
k<-1
# 1軸の領域長
L<-1
# 生起点の数
N<-1000
# 生起点の座標
# ポアッソン過程
X<-matrix(runif(k*N)*L,ncol=k)

plot(X,rep(0,N),cex=0.01,pch=19)

# 拡大
# 拡大してプロットします
# 半径
# 最大半径
Rmax<-0.01
r<-runif(1)*Rmax

# 中心
Ctr<-runif(k)*L*(1-2*r)+r


plot(X,rep(0,N),cex=1,pch=19,xlim=c(Ctr[1]-r,Ctr[1]+r))
  • 2次元
    • 雨が地面にぱらぱらと降るような感じです。ある時間に雨粒が落ちた点をプロットしてみることにします

Pois2d.jpeg

# 次元
k<-2
# 1軸の領域長(L^kの領域に「ポツリポツリ」を起こします)
L<-1
# 生起点の数
N<-10000
# 生起点の座標
# ポアッソン過程
X<-matrix(runif(k*N)*L,ncol=k)

plot(X,cex=0.01,pch=19)

# 拡大してプロットします
# 拡大してプロットします
# 半径
# 最大半径
Rmax<-0.1
r<-runif(1)*Rmax

# 中心
Ctr<-runif(k)*L*(1-2*r)+r


plot(X,cex=1,pch=19,xlim=c(Ctr[1]-r,Ctr[1]+r),ylim=c(Ctr[2]-r,Ctr[2]+r))
  • 3次元
    • 次元の拡張は同様です
    • 単位k次元体積当たりの発生数が均一であるような生起パターンであればポアッソン過程である

PoisDF3.jpeg


# 次元
k<-3
# 1軸の領域長
L<-1
# 生起点の数
N<-1000
# 生起点の座標
# ポアッソン過程
X<-matrix(runif(k*N)*L,ncol=k)

plot(X,cex=0.01,pch=19)

# 拡大
# 最大半径
Rmax<-0.1
r<-runif(1)*Rmax

# 中心
Ctr<-runif(k)*L*(1-2*r)+r
library(rgl)
plot(X,cex=1,pch=19,xlim=c(Ctr[1]-r,Ctr[1]+r),ylim=c(Ctr[2]-r,Ctr[2]+r))

plot3d(X,cex=1,pch=19,xlim=c(Ctr[1]-r,Ctr[1]+r),ylim=c(Ctr[2]-r,Ctr[2]+r),zlim=c(Ctr[3]-r,Ctr[3]+r))

ポアッソン過程の最近接距離と次元

  • #ポアッソン過程で何を問題にするかにて、ポアッソン過程の何を問題にするかを見ていくが、その中に、ポアッソン過程の事象間の最短距離を問題にする場合が出てくる。その話題の中で指数分布とワイブル分布との関係に触れるが、それに関する関連記事はこちら

ポアッソン過程で何を問題にするか

いろいろな分布の取り出し方

  • ポアッソン過程で、事象を「時空間」上の点として観察したとする
    • そのプロットから、複数の分布を取り出してみる
  • 取り出す分布
    • 単位体積当たりの点の数の分布
      • 「時空間」を格子状に区切る場合(離散的な扱い)
      • 「時空間」に任意の単位体積を切り取る場合(連続的な扱い)
    • 近接距離の分布
      • ある点から、最も近い点までの距離の分布
        • 離散的格子状「時空間」での距離の場合
        • 連続的「時空間」での距離の場合
      • ある点から、第i番目に近い点までの距離の分布
        • 離散的格子状「時空間」での距離の場合
        • 連続的「時空間」での距離の場合
    • 「はじめて起きる」までの回数の分布(離散的、1次元)
      • 「はじめて」起きるまでの回数の分布
      • 第i回目が起きるまでの回数の分布

「時空間」を格子状に区切って、格子1か所あたりの生起回数を数える

  • 格子状に仕切る

Grid.jpeg

# 次元
k<-2
# 1軸の領域長(L^kの領域に「ポツリポツリ」を起こします)
L<-1
# 生起点の数
N<-200
# 生起点の座標
# ポアッソン過程
X<-matrix(runif(k*N)*L,ncol=k)

# グリッド
plot(X,cex=1,pch=19,xlim=c(0,1),ylim=c(0,1))
t<-seq(from=0,to=1,length=11)

for(i in 1:length(t)){
	abline(h=t[i],col=2)
	abline(v=t[i],col=2)
}
  • 仕切りごとに点の数を数えてヒストグラムを描く

GridHist.jpeg

# 正方形内の点の個数を数える
cnt<-c()
for(i in 1:(length(t)-1)){
	for(j in 1:(length(t)-1)){
		cnt<-c(cnt,length(which(X[,1]>=t[i] & X[,1]< t[i+1] & X[,2]>=t[j] & X[,2]<t[j+1])))

	}
}
hist(cnt)

「時空間」に小円を描いて、その内部の点を数える

  • 円を描く・数える

CirclePois.jpeg CircleHist.jpeg

# 円
plot(X,cex=1,pch=19,xlim=c(0,1),ylim=c(0,1))
Ncircle=15

t<-seq(from=0,to=1,length=100)*2*pi
r<-0.05
for(i in 1:Ncircle){
	par(new=TRUE)
	ctr<-runif(2)*0.9+0.05
	par(new=TRUE)
	plot(r*cos(t)+ctr[1],r*sin(t)+ctr[2],col=2,xlim=c(0,1),ylim=c(0,1),type="l")
}
Ncount<-100
cnt<-c()
for(i in 1:Ncount){
	ctr<-runif(2)*0.9+0.05
	D<-apply(t(t(X)-ctr)^2,1,sum)
	cnt<-c(cnt,length(which(D<r^2)))
}
hist(cnt,breaks=0:ceiling(max(cnt)))

近接点を図示、距離の分布を取る

  • 個々の点から最も近い点を図示する

Closest.jpeg

# 最近接点

plot(X,cex=1,pch=19,xlim=c(0,1),ylim=c(0,1))
Npt<-15

for(i in 1:Npt){
	pt<-X[i,]
	D<-apply(t(t(X)-pt)^2,1,sum)
	Dorder<-order(D)
	#print(Dorder[1])
	segments(pt[1],pt[2],X[Dorder[2],1],X[Dorder[2],2],col=2)
	points(pt[1],pt[2],col=2,cex=1.5)
}
  • 個々の点から1,2,3番目に近い点を図示する
  • 1,2,3番目に近い点の距離の分布を描く
    • 1,2,3番をそれぞれ、赤、緑、青とする

Closest3.jpeg Dist123.jpeg

# 第i番近接点

k<-3
plot(X,cex=1,pch=19,xlim=c(0,1),ylim=c(0,1))
Npt<-15
for(i in 1:Npt){
	pt<-X[i,]
	D<-apply(t(t(X)-pt)^2,1,sum)
	Dorder<-order(D)
	#print(Dorder[1])
	for(j in 1:k){
		lw<-0.5
		if(j==k)lw<-2
		segments(pt[1],pt[2],X[Dorder[j+1],1],X[Dorder[j+1],2],col=j+1,lwd=lw)
	}
	
	points(pt[1],pt[2],col=2,cex=1.5)
}
Ncount=100
distmat<-matrix(0,Ncount,k)
for(i in 1:Ncount){
	pt<-X[i,]
	D<-apply(t(t(X)-pt)^2,1,sum)
	Dorder<-order(D)
	distmat[i,]<-sqrt(D[Dorder[2:(k+1)]])
}
par(mfcol=c(1,3))
for(i in 1:k){
	hist(distmat[,i])
}

par(mfcol=c(1,1))

ポアッソン過程にまつわる分布

ポアッソン過程で何を集計するか

  • 前節では、ポアッソン過程のデータから、いくつかの異なった切り口での集計をした
  • 集計結果は分布を見せる
  • 切り口の実例と、切り口がもたらす分布について、列挙する
  • 単位体積あたりの件数を数える~ポアッソン分布
    • 起きうる可能性が0回から\infty まであるときに、何回起きるかの確率の分布
    • 珍しいこと(たとえば突然変異)がゲノム上でおきるとき、ある長さのあたりに起きる突然変異の回数は(長さによる上限はあるが、それほどたくさんの変異が起きる可能性はほぼ0なので)0から\infty を想定しているとみなしてもよいとする。そのような過程では、もし、ゲノムの場所によらずに突然変異が起きる確率が等しいならば、単位長さあたりの回数の分布はポアッソン分布に従うだろう。もし、従わないとすれば、「平等」のルールか「相互に無関係」のルールのどちらか片方、もしくは、両方が守られていないと予想される
    • Wikipediaの記事はポアッソン分布
    • ポアッソン分布の確率密度関数

Dpois.jpeg

xs<-0:10
ks<-1:10
for(i in 1:length(ks)){
	plot(xs,dpois(xs,ks[i]),ylim=c(0,0.5),col=ks[i],type="b")
	par(new=TRUE)
}
par(new=FALSE)

生起間隔を測る(1次元)~指数分布

  • 単位体積(1次元の場合は単位長さ。1次元が「ゲノム」のようにものならば、その「長さ」、1次元が時間であれば、時間の「長さ)当たりの集計をするためには、あらかじめ、「地図」を作って、「範囲」を決めておくことが必要になる
  • 他方、「間隔」を測る場合には、生起を追跡していくことで集められる情報である点が、データ収集上、やや異なる
  • 実際、時間経過を追いかけるときには、「観測」という行為は、時刻とともに流れているので、「生起間隔」を集計するという作業は、不自然な方法ではない
  • 1次元時空間でのポアッソン過程の生起間隔は指数分布である
  • ポツリポツリと起きる現象(奇病の発生や、発作)に、何かしらきっかけがある(平等ではない)のではないか、とか、相互に関係がある(誘発関係がある(伝染する、も含む))ことを疑えば、間隔が指数分布に従うかどうかを調べることで判断材料が得られる可能性がある
  • 『間隔が短い確率は長い確率よりも高い』という性質がある

Dexp.jpeg

xs<-seq(from=0,to=10,length=100)
ks<-2^(seq(from=-3,to=3,length=7))
for(i in 1:length(ks)){
	plot(xs,dexp(xs,ks[i]),ylim=c(0,2),col=i,type="l")
	par(new=TRUE)
}
par(new=FALSE)

起きるまでの回数~幾何分布・負の二項分布

  • 「起きる、起きない」を1回ごとに確認するような場合には、「何回、続けて起きなかった後で、初めて起きるだろう」とか「n回起きるまで待つことにしているが、n回目が起きるのは、何回目だろうか」というような待ち方をすることがある
  • 症例数を一定数にまで到達するノルマがあるときに、「n例目が現れるまでにかかる日数は」とか
  • ここで「日数」を数える話が出たが、このように、「離散的」に連続する「時空間」での生起するかしないか(ベルヌーイ過程)を追いかけている場合に相当する
  • 「初めて起きるまで」の方を幾何分布、「n回目が起きるのは何回目だろうか」の方を負の二項分布という
  • より詳しくは、こちらの記事を参照

1次元ポアッソン過程~指数分布を一般化する

  • 一般化の話はmd_一般化する
  • 1次元ポアッソン過程はポアッソン過程の基本形である
  • 生起間隔は指数分布になる
  • 指数関数の確率密度関数
    • 確率密度関数はP(X=x)=\lambda e^{{-\lambda x}}と表される
    • これを2つの方法で一般化する
      • 一般化(1)~ガンマ分布へ
        • P(X=x)=\lambda e^{{-\lambda x}}=\lambda x^{{1-1}}e^{{-\lambda x}}=C\lambda x^{{k-1}}e^{{-\lambda x}};k=1(Ck,xによらない定数)と変形すると
        • C\lambda x^{{1-1}}e^{{-\lambda x}}kに関する一般式として、指数分布の確率密度関数はk=1の特殊な場合と読み取れる
        • この形であって、確率密度関数としての条件(積分して1)を満足するのがガンマ分布である
        • ガンマ分布はポアッソン過程で事象が複数回起きることによって引き起こされる現象がおきるまでの時間の分布である。より詳しくは#ばらばらな出来事の組み合わせがもたらすこと
# 第i番近接点

k<-1
L<-1
# 生起点の数
N<-10000
# 生起点の座標
# ポアッソン過程
X<-runif(N)*L

ns<-1:15
Ncount=1000
distmat<-matrix(0,Ncount,length(ns))
for(i in 1:Ncount){
	pt<-X[i]
	D<-abs(X-X[i])
	Dorder<-order(D)
	distmat[i,]<-D[Dorder[ns]]
}
par(ask=TRUE)
for(i in 1:length(ns)){
	hist(distmat[,i],breaks=seq(from=0,to=0.01,length=1001),ylim=c(0,100))
}
      • 一般化(2)~ワイブル分布へ
        • 指数分布確率密度関数はP(X=x)=\lambda e^{{-\lambda x}}
        • g(x)=\lambda xと置くと、{\frac  {dg(x)}{dx}}=-\lambda なので、P(X=x)=-{\frac  {g(x)}{dx}}e^{{-g(x)}}で表せる
        • 今、g(x)=(\lambda x)^{k}とすれば、指数分布はk=1の場合に相当する
        • P(X=x)=C(-{\frac  {dg(x)}{dx}})e^{{-g(x)}};g(x)=(\lambda x)^{k}の形の確率分布をワイブル分布と言う
        • k>1の場合には、xが大きくなるにつれ確率が大きくなる。システムが時間とともに故障しやすくなる様子などに当てはまる。生体をシステムと考えてもあてはまる面がある
        • k<1<1の場合には、xが小さいうちに確率が大きく、xが大きくなると確率が小さくなる。システムで言えば、できたてのときには、未検出の問題が顕在化して故障しやすいが、時間とともに、それらの影響がなくなって落ち着く様子に当てはまる。生体をシステムと考えても、当てはまる面がある
        • ワイブル分布は、k次元ポアッソン過程の最近接点距離分布と関係することが、こちらの記事にある

ばらばらな出来事の組み合わせがもたらすこと

  • ガンマ分布はポアッソン過程で、何回か繰り返して起きるまでの間隔であると説明した
    • これをグラフ的に考えてみる
    • 「4回起きるまで」の経過をシミュレーションしてみよう
      • ノード0から出発して、確率pでノード1へ移動し、確率1-pでノード0にとどまる
      • ノード1,2,3でも同様で
      • ノード4はゴールである

Step4.jpeg

        • (この図を描くコードはstep4)
      • シミュレーションソース
# ステップ数に上限を入れる
# 低い確率ではあるが、最終ノードに到達しない確率が0ではないため
Nmax<-10000
# シミュレーション試行回数
Niter<-1000
# 最終ノードに到達するまでのステップ数(時間)の記録用オブジェクト
result<-rep(0,Niter)
# 1つ先のノードに進む確率
p<-0.1
# シミュレーション
for(i in 1:Niter){
	node<-0
	cnt<-0
	while(node<Nnode & cnt<Nmax){
		cnt<-cnt+1
		if(runif(1)<p)node<-node+1
	}
	result[i]<-cnt
}
# 要したステップ数で度数分布を描く
hist(result)

GammaSim.jpeg

  • 単純なポアッソン過程を振り返る
    • 複数ステップのシミュレーションはグラフ上の移動とゴールまでのステップ数でシミュレーションできた
    • 「初めて起きるまで」はグラフで表せば

SimplePois.jpeg

    • 度数分布にすれば

SimplePoisHist.jpeg

  • 逆にグラフを複雑にしても、「スタート」と「ゴール」が設定できれば、同様にシミュレーションで、「到達するまでのステップ数(~時間)」をシミュレーションすることができることがわかる
    • このようなグラフでの状態推移をシミュレーションするには、ノードからノードへの移動を行列を使って計算することが有用で、そのような行列を推移行列・確率行列と言う(英文Wikipedia記事)。これについてはmd_周期的安定で扱う。

CompleGraph.jpeg

2(多)種類のばらばらな出来事の相互関係

  • 2種類のポアッソン過程が起きているときに、タイプAが起きてからタイプBが起きるまでのタイムラグを測定するとどのような分布になるだろう
  • 簡単のためにどちらも同じ確率で起きるものとする
    • A,B2種類の事象が起きていて、AがBよりも後に起きたときを記録し、AがBよりもどれくらい遅く起きるかの分布を調べるものとする
# シミュレーション回数
N<-1000000
A<-runif(N) # Aの生起時刻
B<-runif(N) # Bの生起時刻
# タイムラグ
D<-A-B
# AがBの後におきた場合のみを取り出す
Dpositive<-D[which(D>0)]
# タイムラグを度数分布化
hist(Dpositive)

Dpositive.jpeg

  • Bが起きた直後の方が、Bが起きてしばらく経ってからよりも、Aが起きやすい
  • このような観察は、Aを引き起こすかもしれない介入B(治療など)の後にAが起きるかどうかを観察しているときに相当する
  • 「A,Bが相互にばらばら」なときのタイムラグがこのようになるから、BがAを誘発して、Bの直後にAが起きがちである場合には、この分布よりもタイムラグが「短め」に偏ることがわかる

確率分布を表す関数

確率分布についてのその他の読み物

  • こちらは、本モジュールの元記事。多少、異なった表現をしているので、参考まで

Rの確率分布関数と乱数発生

こちらのサイトで</math>