読者です 読者をやめる 読者になる 読者になる

廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

(2次元)ポアソン過程のシミュレーション

R 確率過程

アルゴリズム

Jochen Voss によればポアソン過程の乱数を生成するアルゴリズムは以下の通りである.

(p.60 参照)

入力(input):

  • 強度関数(intensity function) \lambda : \mathbb{R}^d \to   \mathbb{R}
  • 集合  D \subseteq \mathbb{R}^d とその関数  \Lambda(D) < \infty

使用する乱数(randomness used):

  •  N \sim \rm{Pois}(\Lambda(D))
  • 独立同分布のサンプル  X_i \sim \frac{1}{\Lambda(D)}1\!\!1_D \lambda(\cdot)i = 1,…, N

出力(output):

  • 集合 D 上のポアソン過程. その強度(intensity)は λ

擬似コード

  1. 乱数  N \sim \rm{Pois}(\Lambda(D))
  2.   \Pi \leftarrow \phi空集合
  3. for i = 1, 2, … N. do
  4.  X_i \sim \frac{1}{\Lambda(D)}1\!\!1_D \lambda(\cdot) を生成.
  5.   \Pi \leftarrow \Pi \cup \{X_i\}
  6. end for
  7. return   \Pi

ここで,

\displaystyle \Lambda(D) =\int_D \lambda(x) \, dx,

\displaystyle \frac{1}{\Lambda(D)}1\!\!1_D \lambda(x) =\left\{ \begin{array}{ll}\frac{\lambda(x)}{\Lambda(D)} & (x \in D)\\0 & (\mbox{その他})\end{array}\right.,

である.

\displaystyle D =[0,1]\times[0,1] = \subseteq \mathbb{R}^2,

\displaystyle \lambda(x) =1000,

とすると,

\Lambda(x) =\int_D 1000 \,dx =\int^{1}_{0}\int^{1}_{0} 1000 \,dx_1 dx_2 =1000

より,

\displaystyle \frac{1}{\Lambda(D)}1\!\!1_D \lambda(x) =\left\{ \begin{array}{ll} 1 & (x \in D)\\0 & (\mbox{その他})\end{array}\right.,

となる.

f:id:abrahamcow:20150108154113p:plain

よって一様乱数で2次元のポアソン過程が作れる.

N <- rpois(1,1000)
x<- runif(N)
y <- runif(N)
plot(x,y)

f:id:abrahamcow:20150108154641p:plain

library(spatstat) にはこのような関数が用意されている.

library(spatstat) 
N <- rpois(1,1000)
X <- runifpoint(N)
plot(X)

f:id:abrahamcow:20150108155329p:plain

http://web.sfc.keio.ac.jp/~maunz/ESTRELA/200912_77.pdf