2014-11-19

シミュレーション〜ぱらぱらめくるRandomFieldsパッケージのvignette

シミュレーションでは、大きく２つの手法を採用している
- コレスキー分解
  - 地点ごとの値の似ている具合が共分散行列になるので、それから多変量正規分布を使うなどするに際してコレスキー分解
  - 発生地点数ｘ変数の数が増えるとどんどん大変になるのは、想像がつきやすい
- Circulant embedding
  - 全地点を個々に扱おうとすると大変なので、周期関数的に扱ってフーリエ変換が使えるようにする、という戦略
推定方法はいろいろ採用している。どれもなかなかに難しい
- Least squares
- Maximal likelihoodとその変法
  - 最適化が難しい、初期値問題もある
- Kriging
モデル
- Cross-covarianceモデルとCross-variogramモデルとがある
- Cross-covarianceモデル
  - Whittle-Matern モデル
    - 2点間の距離のみに依存するモデル $W_{\nu}(h)=\frac{1}{2^{\nu-1}\Gamma(\nu)}|h|^{\nu}K_{\nu}(|h|)$ という式で表される。ごちゃごちゃしているけれど、関数の形を決めるパラメタは $\nu$ だけで、あとは、ガンマ関数とmodifiedベッセル関数(K)のみで決まり、値を決める変数は、2点間の距離 $|h|$ のみで決まるという単純なもの
    - 距離に応じてどんな感じで関係の強さが減衰するかを形状パラメタ $\nu$ だけで変えられている様子がわかる

library(RandomFields)
h <- seq(from=0,to=3,length=100)
nus <- seq(from=0,to=3,length=50)
nus <- nus[-1]
W <- matrix(0,length(nus),length(h))
for(i in 1:length(nus)){
	nu <- nus[i]
	W[i,] <-1/(2^(nu-1)*gamma(nu)) * h^nu * besselK(h,nu)
}
matplot(t(W),type="l")

- - - その特別形
  - Exponential モデル
  - Gaussian モデル
  - Spherical モデル他
使ってみる
- シミュレーション
- 空間次元d=2、変数の数m=1とする
- Whittle-Matern covariance functionは、距離にのみ応じて関連が減少する１変数ランダムフィールドモデルであって、１パラメタ $\nu$ で指定できる

nu <- 0.3
model <- RMwhittle(nu=nu)
# シミュレーション範囲
x <- y <- seq(from=-3,to=3,by=0.3)
# シミュレーション
simu <- RFsimulate(model,x,y)
plot(simu)

- 変数の数を増やす。空間次元d=2、変数の数m=3とする
- 空間にはWhittele-Matern covariance functionがあるとして、２つの変数の分布があるとする。どちらも分散が１とする
  - 一つのランダムフィールドを生成して、それをそのまま使えば、もちろん、２変数の値の分布は同じ

nu <- 2
M1 <- c(1,1)
model <- RMmatrix(M=M1,RMwhittle(nu=nu))
x <- y <- seq(from=-3,to=3,by=0.3)
simu <- RFsimulate(model,x,y)
plot(simu)
plot(simu@data)

- - ２つの変数のそれを違えようと思えば、変数１と変数２とにそれぞれ別のランダムフィールドを作ればよい

nu <- 2
M1 <- c(1,0)
M2 <- c(0,1)
model <- RMmatrix(M=M1,RMwhittle(nu=nu)) + RMmatrix(M=M2,RMwhittle(nu=nu))
x <- y <- seq(from=-3,to=3,by=0.3)
simu <- RFsimulate(model,x,y)
plot(simu)
plot(simu@data)

- - これでは、２変数のばらつき加減は同じにしても、２変数の値に関連が出てこない。なぜなら、同じ $\nu$ を使って、別々に発生させたランダムフィード２つ、そのものだから。混ぜればよい。混ぜるにあたって、一つ目の変数は、一つ目のランダムフィールドと二つ目のランダムフィールドを２：１で混ぜ、二つ目の変数は、１：３で混ぜるとしよう。分散の関係から、 $\sqrt{\frac{2^2}{2^2+1^2}}:\sqrt{\frac{1^2}{2^2+1^2}}=\sqrt{\frac{4}{5}}:\sqrt{\frac{1}{5}}$ と $\sqrt{\frac{1^2}{3^2+1^2}}:\sqrt{\frac{3^2}{3^2+1^2}}=\sqrt{\frac{1}{10}}:\sqrt{\frac{9}{10}}$ に分けるとする

nu <- 2
M1 <- c(sqrt(4/5),sqrt(1/10))
M2 <- c(sqrt(1/5),sqrt(9/10))
model <- RMmatrix(M=M1,RMwhittle(nu=nu)) + RMmatrix(M=M2,RMwhittle(nu=nu))
x <- y <- seq(from=-3,to=3,by=0.3)
simu <- RFsimulate(model,x,y)
plot(simu)
plot(simu@data)

- 分布がシフトする
  - 距離に応じて相関が減ずるモデルがRMstable()モデル。距離の何乗か、それを拡縮するスケールはいくつか、で決まるモデル
  - このRMstable()モデルのランダムフィールドを１つ使って、それをある方向にずらすことができる(雲のパターンが１時間後にどうなるか、みたいな)
  - ずれのベクトルが(4,3)になっているのを確認し

- - そのうえで、１枚の絵の中での、地点間の値のバラツキ具合と、２枚の絵を比較したときの値の相関加減についてのモデル等高線とを示す。ベクトル(4,3)に応じたところに、相関のピークが出るのがわかる

x <- y <- seq(from=-10,to=10,by=0.3)
model <- RMstable(alpha=1.8,scale=2)
model2 <- RMdelay(model, s = c(4,3))
simu <- RFsimulate(model2,x,y)
plot(simu)

plot(model2,dim=2,xlim=c(-5,5))
abline(v=c(0,4))
abline(h=c(0,3))
|