はじめに

動物の移動はGPS発信機を装着し、位置取得することは可能です。しかし、その時に何をしていたのかということは分かりません。
また、GPSで取得した位置は計算過程やノイズによって誤差が生じ、本当にいた位置とは異なります。特にGPSで取得した動物の移動データは誤差が大きいものが多く、どうすれば本当の位置データを得ることができ、何をしていたのか知ることができるのだろうと疑問に思ったことがありました。その時出会った手法が、MCMCを用いてSCRWモデルのパラメータを求め、真の位置と移動モードとの推定する手法でした*1。論文を読むにあたって、MCMCやSCRWモデルや状態空間モデルについて、初歩的なところから理解する必要があると感じました。
そこで、RとjagsでMCMCを用いてSCRWモデルのパラメータを推定し、真の位置と移動モードを推定する方法を勉強したので、自分が本当に理解できているのかの確認も含めてアウトプットしてみようと思います。*2

SCRWモデル

ちなみに、SCRWとはSwitching Correlated Random Walkの略です。 Correlated Random WalkモデルとSwitchingモデルを組み合わせたモデルのことをSCRWモデルと呼びます。 Correlated Random Walkモデルとは、ある時刻の点が1つ前の時刻の点の速度の影響を受け、ある一定の角度で曲がるモデルを指します。
また、Switchingモデルとは、ある状態から別の状態へ遷移するモデルを指します。
例えば動物の移動の場合、ある場所に向かって移動する場合と索餌を行い大きな移動がない2パターンが考えられます。ある場所に向かって移動する場合は、移動角度が小さく、前の移動速度の影響を多く受けます。索餌の場合は、移動角度が大きく、前の移動速度の影響をあまり受けません。
今回は移動と索餌の2つの状態があると仮定し、MCMC（マルコフ連鎖モンテカルロ法）して状態の遷移確率・移動と索餌の際の前の速度の影響および移動角度を推定します。

jagsでは、.bugという拡張子でモデルを定義し、そのモデルを読み込んでMCMCします。
SCRWモデルは以下のように定義しました*3。

model{
  # def pi
  pi <- 3.141592653589;
  npi <- -1 * pi;
  # 移動の際の誤差
  # 緯度経度方向に正規分布の観測誤差が発生すると仮定する
  sigma_lon ~ dnorm(1, 0.1);
  sigma_lat ~ dnorm(1, 0.1);
  # rhは相関係数
  rh ~ dunif(-1, 1);
  # 多変量正規分布
  Sigma[1, 1] <- sigma_lon * sigma_lon;
  Sigma[1, 2] <- rh * sigma_lon * sigma_lat;
  Sigma[2, 1] <- rh * sigma_lon * sigma_lat;
  Sigma[2, 2] <- sigma_lat * sigma_lat;
  iSigma[1:2, 1:2] <- inverse(Sigma[,])
  # xは真の位置、yは観測値とする
  # 1st position
  for(k in 1:2){
    x[1, k] ~ dnorm(y[1, k], 1);
  }
  # 2nd position
  x[2, 1:2] ~ dmnorm(x[1,], iSigma[,]);
  # theta[1]は移動の時の移動角度、theta[2]は索餌の時の移動角度
  # 移動の時は移動角度が小さい（まっすぐ移動する）はず
  theta[1] ~ dunif(-0.1, 0.1);
  # 索餌の時は移動角度が大きい（ぐるぐる移動する）はず
  theta[2] ~ dunif(npi, pi);
  # gamma[1]は移動の時の移動速度の影響
  # gamma[2]は索餌の時の移動速度の影響
  # 移動の時は前の移動速度の影響を大きく受けるはず
  gamma[1] ~ dunif(0.5, 1);
  # 索餌の時は前の移動速度の影響をほとんど受けないはず
  gamma[2] ~ dunif(0, 0.5);
  # alpha[1]は移動の時の次の状態が移動となる確率
  # alpha[2]が索餌の時の次の状態が索餌となる確率
  alpha[1] ~ dbeta(1, 1);
  alpha[2] ~ dbeta(1, 1);
  # 初期値
  lambda[1] ~ dunif(0, 1);
  lambda[2] <- 1 - lambda[1];
  # bmode[]は移動モードを表す
  bmode[1] ~ dcat(lambda[]);
  #
  for(t in 2:(N-1)){
    phi[t, 1] <- alpha[bmode[t-1]];
    phi[t, 2] <- 1 - alpha[bmode[t-1]];
    # カテゴリカル変数
    bmode[t] ~ dcat(phi[t,]);
    # 回転行列
    T[t, 1, 1] <- cos(theta[bmode[t]]);
    T[t, 1, 2] <- (-sin(theta[bmode[t]]));
    T[t, 2, 1] <- sin(theta[bmode[t]]);
    T[t, 2, 2] <- cos(theta[bmode[t]]);
    for(k in 1:2){
      Tdx[t, k] <- T[t, k, 1] * (x[t, 1] - x[t-1, 1]) + T[t, k, 2] * (x[t, 2] - x[t-1, 2]);
      x.mn[t, k] <- x[t, k] + gamma[bmode[t]] * Tdx[t, k];
    }
    x[t+1, 1:2] ~ dmnorm(x.mn[t,], iSigma[,]);
  }
  # Measurement equation
  # 観測値は真の位置に誤差を加えたものになる
  itau ~ dgamma(0.001, 0.001);
  for(t in 2:N){
    for(k in 1:2){
      y[t, k] ~ dnorm(x[t, k], itau);
    }
  }
}