ESTRELA連載 第77回

今回使うデータ

データの作成と表示

  • パッケージの読み込み
    library(spatstat)
  • データの読み込み
    X600 <- read.table("X600.txt", sep=",", header=T)
    m <- sample(c("c","m","y","b"), 600, replace=TRUE)
    m <- factor(m, levels=c("c","m","y","b"))
    X <- ppp(X600$x, X600$y, c(0,1), c(0,1), marks=m)
  • データ作成例
    x <- runif(600)
    y <- runif(600)
    m <- sample(c("c","m","y","b"), 600, replace=TRUE)
    m <- factor(m, levels=c("c","m","y","b"))
    X <- ppp(x, y, c(0,1), c(0,1), marks=m)
  • データ表示
    plot(X)
    plot(split(X)$b, main="color dots")
    points(split(X)$m, col=14)
    points(split(X)$y, col=7)
    points(split(X)$c, col=5)

密度の計算

  • 密度
    plot(density(X), 0.1)
    contour(density(X), add=T)
    plot(density(split(X)$c), 0.1)
    contour(density(split(X)$c), add=T)
    D <- density(X)
    persp(D)
  • 最小二乗誤差法によるバンド幅の決定
    library(splancs)
    poly0 <- cbind(c(0,1,1,0), c(0,0,1,1))
    Mse2d <- mse2d(as.points(X), poly0, nsmse=50, range=1)
    plot(Mse2d$h[5:50],Mse2d$mse[5:50], type="l")
    h.mse2d <- min(Mse2d$h)
  • 交差検証対数尤度関数によるバンド幅の決定
    library(spatialkernel)
    x <- runif(600)
    y <- runif(600)
    mks <- sample(c("a","b", "c"), 600, replace=TRUE)
    pts <- cbind(x, y)
    h <- seq(0.01, 1, by=0.01)
    cv <- cvloglk(pts, mks, h=h)$cv
    plot(h, cv, type="l")
    hcv <- h[which.max(cv)]
  • コドラート法
    quadratcount(X, nx=5, ny=5)
    plot(X)
    plot(quadratcount(X, nx=5, ny=5), add=T)
  • ティーセン分割
    plot(dirichlet(X))
  • ドロネー三角網図
    plot(delaunay(X))

ランダム性の検定 (Complete Spatial Randomness: CSR)

  • コドラート法によるランダム性のχ2検定
    # all points
    quadrat.test(X, nx=5, ny=5)
    # plot(X)
    # plot(quadrat.test(X, nx=5, ny=5), add=T)
    plot(quadrat.test(X, nx=5, ny=5),  col=2)
    points(X, col=16, cex=0.2)
    plot(quadrat.test(X, nx=5, ny=5),  col=2, add=T)
    # 'cyan' points
    quadrat.test(split(X)$c, nx=5, ny=5)
    plot(split(X)$c)
    plot(quadrat.test(split(X)$c, nx=5, ny=5), add=T)
    • 訂正:本文中、「p値が非常に小さいとき帰無仮説を棄却し、点過程がランダムでないと判断します」となっていますが、正しくは「p値が大きいとき帰無仮説を棄却し、点過程がランダムでないと判断します」または「p値が非常に小さいとき帰無仮説を採択し、点過程がランダムであると判断します」です。
  • コルモゴロフ・スミルノフ検定
    xcoord <- function(x,y) {x}
    kstest(split(X)$c, xcoord)
    plot(kstest(split(X)$c, xcoord))
  • ポアソンモデルの適用 λ_θ(x,y)=θ_0+θ_1*x+θ_2*y
    modelc <- ppm(split(X)$c, ~x+y)
    plot(modelc, how="image", se=FALSE)
    plot(predict(modelc, type="trend", ngrid=512))
    • 訂正:本文中「均一なポアソン過程」となっていますが、正しくは「不均一なポアソン過程」です。
  • 境界効果
    x <- runif(100)
    y <- runif(100)
    X <- ppp(x, y, c(0,1), c(0,1))
    Y <- ppp(x,y,c(0.1,0.9),c(0.1,0.9))
    plot(X, pch=19)
    rect(0.1,0.1,0.9,0.9, lty=2)
    plot(distmap(Y), add=T)
    points(Y, col="white", pch=19)
    #
    plot(distmap(X))
    points(X) 
    #
    plot(X)
    plot(density(Y), add=T)
    points(Y, pch=19)
    #
    plot(density(X))
    points(X) 
  • F関数:一定半径以内
    plot(distmap(split(X)$c))
    plot(split(X)$c, add=T)
    plot(Fest(split(X)$c))
    # Fest(X)
    # plot(Fest(X))
    plot(Fest(split(X)$c), theo~r)
  • G関数:最近隣距離
    plot(Gest(X))
  • K関数:ペアワイズ距離
    K <- Kest(X)
    plot(K)
    Kc <- Kest(split(X)$c)
    plot(Kc)
  • L関数:一般に用いられるK関数
    plot(Lest(X))
  • ペア相関関数:K'/πr^2
    plot(pcf(X))
  • J関数:{1-G}/{1-F}
    plot(Jest(X)

不均一なK関数

  • inhomogeneous K-function
    predc <- predict(modelc, locations=X)
    Ki <- Kinhom(X, predc)
    plot(Ki)
    Ki2 <- Kinhom(X)
    plot(Ki2)

マーク付き点過程

  • marked data
    plot(split(X))
    plot(density(split(X)))
    plot(smooth.ppp(X))
  • 要約統計量
    M <- marktable(X, R=0.1)
  • multi-type
  • Gij関数
    plot(Gcross(X, "c", "b"))
    plot(alltypes(X, "G"))
  • Kij関数
    plot(Kcross(X, "c", "b"))
    plot(alltypes(X, "K"))
  • Lij関数
    plot(Lcross(X, "c", "b"))
  • Jij関数
    plot(Jcross(X, "c", "b"))
    plot(alltypes(X, "J"))
  • i-to-any
  • Gdot関数
    plot(alltypes(X, "Gdot"))
  • マーク相関
    eqfun <- function(m1, m2) { m1 == m2}
    M <- markcorr(X, eqfun, correction="translate", method="density", kernel="epanechnikov")
    plot(M)

ランダム性の検定

  • 帰無仮説:マーク付き点過程がポアソン均一過程に従う
    E <- envelope(X, Kcross, nsim=99, i="b", j="c")
    plot(E)
  • 要素の独立性
    E <- envelope(X, Kcross, nsim=99, i="b", j="c", simulate=expression(rshift(X, radius=0.1)))
    plot(E)
  • ランダムラベリング

包絡分析と適合度検定

  • ポイントワイズ包絡分析
    E <- envelope(X, Kest, nsim=99)
    plot(E)
  • モンテカルロ検定
    E <- envelope(X, Kest, nsim=99, global=T)
    plot(E)

ギブズ過程


トップ   編集 凍結解除 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2010-09-16 (木) 22:53:04 (3228d)