SDE のベイズ推定入門

R の YUIMA パッケージに関する詳細は，次の記事も参照：

1 ベイズ推定の実行

1.1 はじめに

確率過程の統計推測を行うための R パッケージ yuima では多次元の SDE \[ dX_t=b_t^\theta(X_t)\,dt+\sigma_t^\phi(X_t)\,dW_t,\qquad X_0\in\mathbb{R}^d, \] に対してパラメータ \(\theta\in\mathbb{R}^{d_\theta},\phi\in\mathbb{R}^{d_\phi}\) の推定を実行することができる（YUIMA の記事も参照）．

YUIMA においてこの \(\theta,\varphi\) の推定に，qmle() による擬似最尤推定量を用いることも，adaBayes() によるベイズ事後平均による推定も実行可能である．

1.2 確率的プログラミング言語 `Stan` との連携

Stan は Hamiltonian Monte Carlo 法を用いた事後分布サンプリングを，確率モデルを定義するだけで実行することができる汎用的な言語である．

加えて，バックグラウンドで C++ を用いているため，非常に高速な MCMC 計算が可能である．

確率的プログラミング言語 Stan については次の記事も参照：

Stan 入門

rstan による Stan の利用

2024-05-17

R 上の Stan インターフェイス

RStan と CmdStanR

2024-09-19

RStan は確率的プログラミング言語 Stan とのインターフェイスを提供するパッケージであり，これを用いることで Stan を通じた HMC を用いた事後分布からのサンプリングが実行できる．

ここでは，yuima から Stan を使うことを可能にする関数 adaStan() の実装を試みる．

1.3 `RStan` とのインターフェイス

具体的には，R において，次のような関数を定義する：

library(yuima)
library(rstan)

yuima_to_stan <- function(yuima){
  excode <- 'data {\n  int N;\n  array[N+1] real x;\n  real T;\n  real h;\n}
parameters {\n'

  for(i in 1:length(yuima@model@parameter@all)){
    excode <- paste(excode, " real", yuima@model@parameter@all[i], ";\n")
  }

  excode <- paste(excode,"\n}")

  excode <- paste(excode,'\nmodel {\n  x[1] ~ normal(0,1);\n  for(n in 2:(N+1)){')

  excode <- paste(excode,
    "\n  x[n] ~ normal(x[n-1] + h *", gsub("x", "x[n-1]", yuima@model@drift), ",sqrt(h) *", gsub("x", "x[n-1]", yuima@model@diffusion[[1]]),");\n  }")

  excode <- paste(excode,'\n}\n')
}

adaStan <- function(yuima){
  # To Do: yuima オブジェクトのスロットの存在をチェック
  # 変数名の例外処理

  excode <- yuima_to_stan(yuima)

  sde_dat <- list(N =  yuima@sampling@n,
    x = as.numeric(yuima@data@original.data), 
    T=yuima@sampling@Terminal,
    h=yuima@sampling@Terminal/yuima@sampling@n)

  fit <- stan(model_code=excode,
    data = sde_dat, 
    iter = 1000,
    chains = 4)

  return(fit)
}

1: Stan モデルのコード（パラメータ部分は未定）を文字列として excode 変数に格納する．
2: ここからが adaStan 関数の本体である．Yuima モデルの全てのパラメータについてループを開始して，excode にパラメータの宣言を追加していく．
3: ここでついに Stan モデルのパラメータの定義部分が完成する．
4: 最後はモデルの定義部分を追加して，Stan モデルのコードが完成する．最初の観測値 x[1] は \(\mathrm{N}(0,1)\) に従う．
5: それ以降の観測値 x[n] は，前の観測値 x[n-1] に drift 項と diffusion 項を加えたものに従う．これを実装するために，Yuima モデルの drift 項と diffusion 項の定義文を呼び出し，x を x[n-1] に置換することで Stan モデルのコードに埋め込む．
6: adaStan という関数を定義する．この関数は，Yuima パッケージのオブジェクトを引数として受け取り，Stan での推定を行い，その結果を fit オブジェクトとして返す．
7: Stan での推定を実行するために，Yuima モデルのデータを Stan モデルに渡すためのリスト sde_dat を作成する．
8: 最後に Stan モデルをコンパイルして実行し，結果を fit オブジェクトとして返す．

これが最初に思いつく最も直接的な方法かも知れないが，このままではいくつかの問題がある：

問題点

関数内部で Stan model のコードを文字列として生成していることが苦しい．

Stan を使う以上，どこかで Stan モデルの情報を受け渡すことは必要になるが，できることならばもっと良い方法を考えたい．

1.4 問題点

adaStan() 関数の挙動を詳しく見るために，次の具体例を考える．

YUIMA を通じて１次元 OU 過程

\[ dX_t=\theta(\mu-X_t)\,dt+\sigma\,dW_t \]

をシミュレーションをするためには，次のようにモデル定義をする：

model <- setModel(drift = "theta*(mu-x)", diffusion = "sigma", state.variable = "x", solve.variable = "x")

これだけで，YUIMA は勝手にパラメータを識別してくれる：

model@drift

expression((theta * (mu - x)))

model@diffusion[[1]]

expression((sigma))

これを通じて生成される Stan モデル文は

data {
  int N;
  real x[N+1];
  real T;
  real h;
}

parameters {
  real<lower=0> theta;
  real<lower=0> mu;
  real<lower=0> sigma;
}

model {
  x[1] ~ normal(0,1);
  for(n in 2:(N+1)){
    x[n] ~ normal(x[n-1] + h * theta * (mu - x[n-1]),
                  sqrt(h) * sigma);
  }
}

となるべきであるが，実際その通りになる：

x <- setYuima(model = model)
stancode <- yuima_to_stan(x)
cat(stancode)

data {
  int N;
  array[N+1] real x;
  real T;
  real h;
}
parameters {
  real theta ;
  real mu ;
  real sigma ;
 
} 
model {
  x[1] ~ normal(0,1);
  for(n in 2:(N+1)){ 
  x[n] ~ normal(x[n-1] + h * (theta * (mu - x[n-1])) ,sqrt(h) * (sigma) );
  } 
}

1.5 `CmdStanR` による解決

library(cmdstanr)
stan_file_variables <- write_stan_file(stancode)
mod <- cmdstan_model(stan_file_variables)

mod$print()

data {
  int N;
  array[N+1] real x;
  real T;
  real h;
}
parameters {
  real theta ;
  real mu ;
  real sigma ;
 
} 
model {
  x[1] ~ normal(0,1);
  for(n in 2:(N+1)){ 
  x[n] ~ normal(x[n-1] + h * (theta * (mu - x[n-1])) ,sqrt(h) * (sigma) );
  } 
}

などとすることで，一時ファイル上で stan ファイルとバイナリファイルを作成・操作することができる．

mod$stan_file()

[1] "/var/folders/7c/j9mzb7pn0wn1k_f9j58c8y480000gn/T/RtmpeJ3FD3/model_3268e3fdcc16f1dc945d0e0e9d6f4e39.stan"

2 Stan インターフェイス

2.1 `RStan` パッケージ

詳しくは次稿参照：

R 上の Stan インターフェイス

RStan と CmdStanR

2024-09-19

前述の OU 過程 1.4

\[ dX_t=\theta(\mu-X_t)\,dt+\sigma\,dW_t \]

で stan 関数でベイズ推定を実行してみます．

パラメータは \[ \begin{pmatrix}\theta\\\mu\\\sigma\end{pmatrix} = \begin{pmatrix}1\\0\\0.5\end{pmatrix} \] として YUIMA を用いてシミュレーションをし，そのデータを与えてパラメータが復元できるかをみます．

sampling <- setSampling(Initial = 0, Terminal = 3, n = 1000)
yuima <- setYuima(model = model, sampling = sampling)
simulation <- simulate(yuima, true.parameter = c(theta = 1, mu = 0, sigma = 0.5), xinit = rnorm(1))

# シミュレーション結果
plot(simulation)

rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
fit <- adaStan(simulation)

print(fit)

Inference for Stan model: anon_model.
4 chains, each with iter=1000; warmup=500; thin=1; 
post-warmup draws per chain=500, total post-warmup draws=2000.

         mean se_mean   sd    2.5%     25%     50%     75%   97.5% n_eff Rhat
theta    0.46    0.15 0.55   -0.15    0.04    0.28    0.70    1.62    13 1.24
mu      -1.62    1.60 9.33  -23.94   -3.00   -0.37    0.30   16.67    34 1.10
sigma    0.49    0.00 0.01    0.48    0.49    0.49    0.50    0.52   197 1.04
lp__  3109.08    0.11 1.16 3106.30 3108.59 3109.28 3109.88 3110.69   116 1.05

Samples were drawn using NUTS(diag_e) at Thu Sep 19 20:42:29 2024.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

plot(fit)

ci_level: 0.8 (80% intervals)

outer_level: 0.95 (95% intervals)

library("bayesplot")
library("rstanarm")
library("ggplot2")

posterior <- as.matrix(fit)
plot_title <- ggtitle("Posterior distributions",
                      "with medians and 80% intervals")

mcmc_areas(posterior,
           pars = c("theta", "mu", "sigma"),
           prob = 0.8) + plot_title

\(\sigma\) の推定はよくできているが \(\mu\) の精度はあまりよくなく，\(\theta\) はバイアスがある様で，自信を持って間違えることも多い．

2.2 `brms` パッケージ

brms や rethinking も，背後で Stan を利用している．これらが文字式をどのように取り扱っているかを調査する．

R によるベイズ混合モデリング入門

brms を用いた混合効果モデリング入門

2024-05-12

Stan コードを扱っている関数は .stancode() であった．

最終的に，.compile_model_rstan() と .fit_model_rstan() が呼ばれるようになっている．

最終的にはこれらの関数も，第 1.3 節のサンプルコードと同様の要領で paste0 を使っていた．

paste0 と paste の違い

# paste0() の使用例
result1 <- paste0("Hello", "world")
print(result1)  # "Helloworld"

[1] "Helloworld"

# paste() の使用例
result2 <- paste("Hello", "world")
print(result2)  # "Hello world"

[1] "Hello world"

result3 <- paste("Hello", "world", sep = "-")
print(result3)  # "Hello-world"

[1] "Hello-world"

3 テンプレート操作

R の Expression について

オブジェクト志向言語ではコード自体もオブジェクトであり，これを R では Expression と呼ぶ．

R（４）メタプログラミング

Expression について

2021-05-07

１つのクラスからなるわけではなく，call, symbol, constant, pairlist の４つの型からなる．¹

次のような操作ができる²

rlang::expr がコンストラクタである：

library(rlang)
z <- expr(y <- x*10)
z

y <- x * 10

expression オブジェクトは base::eval() で評価できる：

x <- 4
eval(z)
y

[1] 40

expression には list のようにアクセス可能である：³

f <- expr(f(x = 1, y = 2))

f$z <- 3
f

f(x = 1, y = 2, z = 3)

f[[2]] <- NULL
f

f(y = 2, z = 3)

4 To Do List

adaStan()

References

Wickham, H. (2019). Advanced r. Chapman & Hall.

Footnotes

(Wickham, 2019) 第17章２節．↩︎
(Wickham, 2019) 第18章 ↩︎
(Wickham, 2019) 第17章２節．↩︎

YUIMA 入門

Stan 入門

R による記号微分入門

1 ベイズ推定の実行

1.1 はじめに

1.2 確率的プログラミング言語 Stan との連携

Stan 入門

R 上の Stan インターフェイス

1.3 RStan とのインターフェイス

1.4 問題点

1.5 CmdStanR による解決

2 Stan インターフェイス

2.1 RStan パッケージ

R 上の Stan インターフェイス

2.2 brms パッケージ

R によるベイズ混合モデリング入門

3 テンプレート操作

R（４）メタプログラミング

4 To Do List

References

Footnotes

1.2 確率的プログラミング言語 `Stan` との連携

1.3 `RStan` とのインターフェイス

1.5 `CmdStanR` による解決

2.1 `RStan` パッケージ

2.2 `brms` パッケージ