RとStanで微分方程式入りのモデル

Tokyo.R #105

伊東宏樹

2023-04-22

自己紹介

名前: 伊東宏樹

勤務先: 森林総合研究所 北海道支所

共訳書:

BUGSで学ぶ階層モデリング入門 生態学のための階層モデリング

ロジスティック方程式

生物の個体数の時間変化をあらわす微分方程式

\[ \frac{dx}{dt} = rx\left(1-\frac{x}{K}\right) \]

  • x : 個体数
  • t : 時間
  • r : 内的自然増加率
  • K : 環境収容力

模擬データの生成

  • 用意したロジスティック曲線にそった値に、対数正規分布にしたがうノイズが加わって点のデータが生成されたとする。
  • 微分方程式を組み込んだモデルで、生成された点のデータから、曲線のパラメータ (r, K) を推定する。

Stanで微分方程式

Stanには常微分方程式の数値解を求める関数が組み込まれている。アルゴリズムの異なる関数がいくつか用意されているが、今回は一般的なode_rk45を使用。

array[] vector ode_rk45(function ode, 
                        vector initial_state,
                        real initial_time,
                        array[] real times, ...)

引数として与える関数の形式

vector ode(real time, vector state, ...)

Stanモデル: functionsブロック

ロジスティック方程式を定義

functions {
  vector logistic(real t, vector x, array[] real theta) {
    vector[1] dx_dt;
    real r = theta[1];
    real K = theta[2];

    dx_dt[1] = r * x[1] * (1 - x[1] / K);
    return dx_dt;
  }
}

data & parametersブロック

data {
  int<lower = 0> N;            // number of measurements
  array[N] real ts;            // measurement times
  real<lower = 0> y0;          // initial measured value
  array[N] real<lower = 0> y;  // measured values
}

parameters {
  real<lower = 0> r;        // intrinsic growth rate
  real<lower = 0> K;        // carrying capacity
  vector<lower = 0>[1] z0;  // initial value
  real<lower = 0> sigma;    // noise scale
}

modelブロック

model {
  array[2] real theta = {r, K};
  array[N] vector[1] z = ode_rk45(logistic, z0, 0, ts,
                                  theta);

  y0 ~ lognormal(log(z0), sigma);
  for (n in 1:N)
    y[n] ~ lognormal(log(z[n]), sigma);
  // priors
  r ~ normal(0, 5);
  K ~ normal(0, 100);
  z0 ~ normal(0, 100);
  sigma ~ normal(0, 5);
}

結果

パラメータの事後平均値をつかってロジスティック曲線を描画(赤線)。

本日の資料

🥳