2022年11月8日火曜日

検出エリア内の滞在時間から移動速度を推定する

ゼミ資料

はじめに

　ちょっとしんどいことがあって，仕事にもあまり手がつかない感じなのだが，研究モチベーションの維持のためにも最近取り組んでいたことをまとめてみた．もう少し落ち着いたら後で振り返ろう．

　昔から「RESTで使う正三角形内での滞在時間から移動速度を出せないのか？」と聞かれることが多い．「やればできるだろうなあ」と思いながら，類似の問題（目的と言うより手法が類似）を，元筑波大のAさんが取り組んでいたはずなので手を出さないでいた（あの研究はどうなったんだろう？）．考えてみると，まあ，我々が使っている正三角形バージョンに限定してしまえば，とくに重なるところもないはず．ちょっとやってみることにした．

問いの確認

　取り組みたい問題は，以下の通り．まず，自動撮影カメラの前に設定した正三角形の検出エリアがある．調査の後に，動画を見ながら，動物がその検出エリア内にとどまった時間（あるいは，通過するのに要する時間．以下，滞在時間）を測定したとする（REST法を使えば，滞在時間と撮影頻度から動物の密度が推定できる）．この滞在時間には，動物の移動速度の情報が含まれていることは自明だろう．もし動物が，1. ランダムな場所から侵入しランダムな場所から脱出する，2. 検出エリア内では真っすぐ動くならば，滞在時間から動物の移動速度の推定も可能な気がする．では，具体的にはどうしたらいいのだろうか？

推定の手順

滞在時間の分解

　最初に考えなければならないのは，測定される滞在時間と速度の関係である．これは，難しい話でもなんでもなくて，\(滞在時間t = 距離y / 移動速度x\)になる．

　さて，距離yや移動速度xは，動物が通過するたびに確率的に変わる値，すなわち，確率変数である．今，移動速度xが従う確率分布をf(x)，距離yが従う確率分布をg(y)としよう．すると，滞在時間tは，2つの確率変数の商が生成する確率分布d(t)に従う確率変数と見なせる．一般に，確率変数yとxの「商」は，以下のような確率分布に従う（導出は省略．途中，ガンマ関数とか出てきてややこしいが，以下のような形に落ち着くはず1）．

\[ d(y/x) = \int y f(x) g(y) dy \]

すなわち，f(x)とg(y)が決まれば，yを積分消去することで，滞在時間が従う確率分布をxの関数として定義できることになる．そこまでいけば，滞在時間データから速度が従う確率分布のパラメータを決めることも可能だろう（今回は，最尤法で推定する）．

移動速度の確率密度関数f(x)

　では，移動速度は，どのような確率分布に従うと想定できるだろうか？これについては，イギリスのRowcliffe大先生がパナマの熱帯雨林で詳細に調べた研究から，対数正規分布に従うことが多いということになっている．おそらく，動物は，時々立ち止まるので，裾野長い分布になるということだろう．対数正規分布は，その逆数も対数正規分布に従う（v ~ dlognorm(mealog, sdlog)の場合，1/v ~ dlognorm(-meanlong, sd)になる）という便利な性質があるので，これはうれしい（今回は使わない）．

移動速度の確率密度関数g(x)

　一方，距離の方も（既存の確率分布にはもちろん従わないが）「ある距離になる確率密度」を与える関数を定義することは不可能ではなさそう．最初，これが結構面倒くさそうに感じていたが，超優秀な東大の院生M君と話していたら，割とエレガントな導出方がありそうだったので，やってみることにする．

　正三角形を考える．動物が底辺と左斜辺の任意の点を通過したとする（下図の赤線がその軌跡）．

この時，動物の軌跡bの長さは，正弦定理から以下のようになる．

\[b^2 = a^2 + b^2 - 2 \times a \times c \times cos(pi/3)\]

このbの長さが，どのような分布になるのかをシミュレーションで確かめてみる．辺の長さは，重要ではないので，以下では１辺の長さ1の正三角形の場合を考えてみる．

n <- 100000
z <- rep(0, n)

for(i in 1:n){
  x <- runif(1, 0, 1) # 辺の任意の点
  y <- runif(1, 0, 1)
  z[i] <- sqrt(x^2 + y^2 - x*y) # bの長さ
}

h <- hist(z, freq = F, breaks = 100)
lines(density(z), col = "blue")

　結構面白い形になる．0.9のちょっと前くらいまでは，直線的に増加していき，そのあと，滑らかに減衰していく．

　このようになる理由は，次のように考えるとわかりやすい．長さの決まった棒があり，その棒の両端を正三角形の辺の上に重ねるという遊びを考える．棒の長さが短いと以下の図のように底辺にぎりぎり沿う線から垂線まで動かせる．棒は，左図の黒丸Dから黒丸Bまで動かせると考えるとわかりやすい．棒の置きやすさは，DBの長さに比例するので，棒が長いほど，うまい具合に重ねやすくなる．

　棒がさらに長くなると（三角形の高さ = sqrt(3)/2より長くなると），垂直の棒が置けなくなってしまう（右図の垂線）．この場合，辺ABとの接点である黒丸が動ける範囲は，AからFの範囲になる．

　以上のように考えると，相対的な棒の置きやすさ（=確率密度に比例する値）は，棒の長さをyとすると，

  • 0 < b < sqrt(3)/2の範囲では，左図のDBの長さ．すなわち，\(g(y) = 2y/\sqrt{3}\times 定数\)

  • sqrt(3)/2 <= b < 1の範囲では，右図のAF

ということになる．

　正弦定理を使えば，AFの長さも容易に求まる．最終的にAFは以下のような形になる（綺麗な式になったらうれしいなと思ったけど，そうでもなかった）．

\[ g(y) = (1 - \sqrt{4y^2 - 3})\times 定数 \]

なお，定数と書いているのは，式で算出されるのは，相対的な生じやすさで，確率密度に比例する値だから．まあ，推定上は無視していい値だが，以下では一応値を出している．

念のためRで確認

　Rで確認する．確率密度にするためには，面積を１にする必要があるので，最初に積分して面積を出して調整する．

# 規格化するための定数の算出（総面積を1にするための準備）
# yが0 ~ sqrt(3)/2の範囲
scale.s <- function(x){
  return(2*sqrt(3)/3*x)
}
integrate(scale.s, 0, sqrt(3)/2) # 念のため確認 sqrt(3)/4に一致
## 0.4330127 with absolute error < 4.8e-15
ss <- sqrt(3)/4

# yがsqrt(3)/2 ~ 1の範囲
scale.l <- function(x){
  return(1 - sqrt(4*x^2 - 3))
}
sl <- integrate(scale.l, sqrt(3)/2, 1)$value # これも頑張れば手計算可能

scale.value <- 1/(ss + sl) # これを掛ければ面積1に


# 距離分布の作成
p.dis.s <- function(x){ # y < sqrt(3)/2
  return(2*sqrt(3)/3*x*scale.value)
}
p.dis.l <- function(x){ # sqrt(3)/2 < y < 1
  #　4*((sqrt(3)/2)^2) - 3 = 0にならないのでroundしておく
  return((1 - sqrt(round(4*(x^2) - 3, 10)))*scale.value) 
}

# 一応，シミュレーション結果と一致することを確認
n <- 100000
z <- rep(0, n)

for(i in 1:n){
  x <- runif(1, 0, 1)
  y <- runif(1, 0, 1)
  z[i] <- sqrt(x^2 + y^2 - x*y)
}

h <- hist(z, freq = F, breaks = 100)
x0 <- seq(0, sqrt(3)/2, by = 0.001)
y0 <- p.dis.s(x0)
lines(x0, y0, col = "red", lwd = 2)

x1 <- seq(sqrt(3)/2, 1, by = 0.001)
y1 <- p.dis.l(x1)
lines(x1, y1, col = "red", lwd = 2)

赤の線が式によって得られたモノ．シミュレーションとぴったりと一致したので，大丈夫そう．

尤度関数の定義

　では，相対的な生じやすさが分かったので，最初の式に戻って，滞在時間の確率密度関数を定義することにする．積分が出てくるが，多分陽には解けないので，Rで近似計算する．関数の中でfor関数を回すのは結構気持ち悪いけど，この辺りは慣れの問題．あと，最尤推定の際に必要なので，対数尤度にマイナスをつけた形で定義している．

# マイナス対数尤度関数の作成
minus.loglf <- function(stay){
  lik <- numeric(0)
  return(function(para){
    for(i in 1:length(stay)){
      int.fs <- function(d){
        d * dlnorm(d/stay[i], para[1], para[2]) * p.dis.s(d)
      }
      int.fl <- function(d){
        d * dlnorm(d/stay[i], para[1], para[2]) * p.dis.l(d)
      }
      bon.d <- sqrt(3)/2
      lik[i] <- integrate(int.fs, lower = 0, upper = bon.d)$value + integrate(int.fl, lower = bon.d, upper = 1)$value
    }
    loglik <- -sum(log(lik))
    return(loglik)})
}

最尤推定

　optim関数を使って，最尤推定してみる．

# データ準備
nsim <- 300 # 多めに300個の滞在時間を作成
meanlog <- 1 # これが正解
sdlog <- 0.2 # これが正解
v <- rlnorm(nsim, meanlog = meanlog, sdlog = sdlog) # 速度が対数正規分布に従う

# 移動距離の作成
dist <- rep(0, nsim)
for(i in 1:nsim){
  x <- runif(1, 0, 1)
  y <- runif(1, 0, 1)
  dist[i] <- sqrt(x^2 + y^2 - x*y) # 距離の乱数
}
hist(dist)

stay <- dist/v # 速度と距離から滞在時間を作成
hist(stay,  breaks = 10)

# 距離分布の作成
p.dis.s <- function(x){ # d < sqrt(3)/2
  return(2*sqrt(3)/3*x)
}
p.dis.l <- function(x){ # sqrt(3)/2 < d < 1
   #　4*((sqrt(3)/2)^2) - 3 = 0にならないのでroundしておく
  return((1 - sqrt(round(4*(x^2) - 3, 10))))
}
res <- optim(par = c(0.8, 0.15), fn = minus.loglf(stay)) # parは初期値
res
## $par
## [1] 1.032135 0.201001
## 
## $value
## [1] 899.061
## 
## $counts
## function gradient 
##       55       NA 
## 
## $convergence
## [1] 0
## 
## $message
## NULL
hist(v, freq = F, ylim = c(0,5), breaks = 30)
curve(dlnorm(x, meanlog, sdlog), add=TRUE, col = 'blue') # 正解
curve(dlnorm(x, res$par[1], res$par[2]), add = TRUE, col = 'green') # 推定線

　青が真の確率分布，青が推定されたパラメータによるもの．無事に一致した．

　ということで，無事，滞在時間から移動速度のパラメータを決めることが出来た．実際の動物は，直線的に移動するとは限らないので，実用的な価値は分からないが，とりあえずできてうれしい．

  1. 「はず」というのは，色んなサイト見ながら最後は自分でやったから．割と間違った（少なくとも，推定が上手くいかない）式が書かれていることもあるので要注意．↩︎

投稿者 時刻:

0 件のコメント:

コメントを投稿

登録: コメントの投稿 (Atom)