純粋跳躍過程の生成作用素と区分的確定的 Markov 過程

ジャンプと確定的な動きによる新たな MCMC 手法

Process
Sampling
R
Author

司馬博文

Published

1/31/2024

Modified

7/02/2024

概要
PDMP は，A 型の Lévy 過程を含む，複合 Poisson 点過程が定めるジャンプと決定論的なドリフトのみからなる確率過程のクラスをいう．この性質をよく理解するために，まずは，有界なレートを持つ純粋に跳躍のみで動く過程の生成作用素を調べる．確率核 \(\mu\) とレート \(\lambda\) という２つのパラメータは，それぞれ各地点からのジャンプ先を定める確率核と，ジャンプの起こりやすさを表す．最後に，現状もっとも活発に研究されている２つの PDMP である Zig-Zag Sampler と Bouncy Particle Sampler とを紹介する．

1 純粋跳躍過程

1.1 導入

次の生成作用素はどのような Markov 過程に対応するか？

\(E\) を距離空間，\(\mu:E\times\mathcal{B}(E)\to[0,1]\) を確率核，\(\lambda\in\mathcal{L}_b(E)_+\) を有界可測関数とする． \[ Af(x):=\lambda(x)\int_E\biggr(f(y)-f(x)\biggl)\mu(x,dy),\quad f\in\mathcal{L}_b(E), \tag{1}\] は有界作用素 \(A\in B(\mathcal{L}_b(E))\) を定め，一様連続半群 \[ \{T_t:=e^{tA}\}_{t\in\mathbb{R}_+}\subset B(\mathcal{L}_b(E)) \] を生成する．1 これに対応する Markov 過程はどのようなものだろうか？

第 1.4 節をご覧ください．

初期分布を \(\nu\in\mathcal{P}(E)\) とする．

1.2 構成１

累積するジャンプを足し合わせた値として，畳み込み半群 \(\{\mu^{\otimes k}\}_{k\in\mathbb{N}}\) に対応する初期分布 \(\nu\) の Markov 連鎖 \(\{Y_k\}_{k=0}^\infty\) を用意する．

これと独立な指数確率変数列 \(\Delta_k\overset{\text{iid}}{\sim}\operatorname{Exp}(1)\) を用いて， \[ X_t:=\begin{cases} Y_0&0\le t<\frac{\Delta_0}{\lambda(Y_0)}\\ Y_k&\sum_{j=0}^{k-1}\frac{\Delta_j}{\lambda(Y_j)}\le t<\sum_{j=0}^{k}\frac{\Delta_j}{\lambda(Y_j)} \end{cases} \] と構成した過程 \(\{X_t\}_{t\in\mathbb{R}_+}\) が，\(\{e^{tA}\}\) に対応する Markov 過程になる．2

なお，\(\lambda(x)=0\) の場合は，ジャンプは起きないもの \(\frac{\Delta}{\lambda(x)}=\infty\) と解する．この場合は零過程である．一般に関数 \(\lambda\in\mathcal{L}_b(E)\) は位置 \(x\in E\) からのジャンプの起こりやすさを表していると思える．

1.3 構成２

\(\lambda=0\) の場合は零過程であるから， \[ \lambda:=\sup_{x\in E}\lambda(x)>0 \] とし，新たな確率核 \(\mu':E\to E\)を \[ \mu'(x,\Gamma):=\left(1-\frac{\lambda(x)}{\lambda}\right)\delta_x(\Gamma)+\frac{\lambda(x)}{\lambda}\mu(x,\Gamma) \] と定めると，生成作用素 \(A\) は \[ Af(x)=\lambda\int_E\biggr(f(y)-f(x)\biggl)\mu'(x,dy) \] とも表せる．

このとき，畳み込み半群 \(\{\mu'^{\otimes k}\}_{k\in\mathbb{N}}\) に対応する初期分布 \(\nu\) の Markov 連鎖 \(\{Y_k'\}_{k=0}^\infty\) は \(\{Y_k\}\) とは分布同等でない．

だが，この \(\{Y_k'\}\) に対しては，独立な Poisson 過程 \(\{V_t\}_{t\in\mathbb{R}_+}\) に対して \[ X'_t:=Y'_{V_t}\quad t\in\mathbb{R}_+ \] と構成される過程 \(\{X_t'\}_{t\in\mathbb{R}_+}\) はやはり \(\{e^{tA}\}\) に対応する Markov 過程である．

命題3
  1. \(\{X'_t\}\) は \(\{X_t\}\) に分布同等である．
  2. \(\{Y'_k\}\) は Markov 性 \[ \operatorname{E}[f(Y'_{k+V_t})|\mathcal{F}_t]=P^kf(X'_t) \] を満たす．ただし，\(P\) は \(\mu'\) が定める作用，\(\mathcal{F}_t:=\mathcal{F}_t^V\lor\mathcal{F}^{X'}_t\) とした．
  3. \(X'\) は \(\{T_t\}\) に対応する Markov 過程である： \[ \operatorname{E}[f(X'_{t+s})|\mathcal{F}_t]=T_sf(X'_t). \]

1.4 まとめ

すなわち \(X\) は非一様な複合 Poisson 過程になる．

各点 \(x\in E\) に於て，レート \(\lambda(x)\) で，分布 \(\mu(x,dy)\) に従う点へとジャンプをする．4

具体的には，次のように整理できる：

強度測度 \(\lambda(x)\mu(x,dy)dx\) を持つ \(E^2\) 上の Poisson 点過程 \(\eta\) に対して， \[ \xi(A)=\int_{A\times E}y\,\eta(dxdy) \] で定まる非一様な複合 Poisson 点過程 \(\xi\) に対応する．

例えば \(E=\mathbb{R}\) のとき， \[ M_t:=\xi([0,t])=\int_0^t\xi(ds) \] が，所定の生成作用素 \(\{T_t\}\) を持つ加法過程になる．

\(\lambda\) を有界としているから，有界区間上でのジャンプ数は有限である．Gamma 過程 のように，\(\mathbb{R}\) の稠密部分集合上でジャンプを繰り返す，というようなことは起こり得ない．

Poisson 過程については次の稿も参照：

Article Image

Poisson 過程を見てみよう

YUIMA パッケージを用いたシミュレーションを通じて

2 区分的確定的 Markov 過程

2.1 導入

前節で調べた純粋跳躍過程に，決定論的な動きを加えた Markov 過程のクラスを，(Davis, 1984) 以来 区分的確定的 Markov 過程 (PDMP: Piecewise Deterministic Markov Process) と呼ぶ．

PDMP は拡散項を持たない Lévy 過程とも理解できる．5

実際，(Davis, 1984) は拡散過程に相補的なクラスとして導入し，(Davis, 1993) は PDMP が待ち行列，ポートフォリオ最適化，確率的スケジューリング，標的追跡，保険数理，資源最適化など，広く確率的モデリングと最適化において重要な役割を演じることを見事に描き出した．6

the class of “piecewise-deterministic” Markov processes, newly introduced here, provides a general family of models covering virtually all non-diffusion applications. (Davis, 1984)

実際，PDMP を用いた MCMC である Piecewise Deterministic Monte Carlo または連続時間 MCMC は，高次元データと大規模データセットに対する効率的なサンプリング法開発の鍵と目されている．

次の記事も参照：

Article Image

新時代の MCMC を迎えるために

連続時間アルゴリズムへの進化

2.2 設定

2.2.1 空間

PDMP は動きのモードが移り変わっても良い．これを表すパラメータ \(v\in K\) を導入する．7

このパラメータ空間上に，PDMP が動く範囲の次元を表す関数 \(d:K\to\mathbb{N}^+\) を導入し，パラメータが \(v\in K\) である際は，\(M_v\subset\mathbb{R}^{d(v)}\) として，積空間 \(M_v\times\{v\}\) 上を運動するものとする．

総合して，PDMP の状態空間を \(M_v\) の直和 \[ E:=\bigcup_{v\in K}M_v\times\{v\} \] と定める．8

2.2.2 フロー

\(\mathfrak{X}_v\) を \(M_v\) 上のベクトル場とし，積分曲線 \[ \frac{d \phi_v(x,t)}{d t}=\mathfrak{X}_v(\phi_v(x,t)), \] \[ \phi_v(x,0)=x \] が存在するとする．9

\(\lambda:E\to\mathbb{R}_+\) は \(s\mapsto\lambda(\phi_v(x,s),v)\) が局所可積分であるような関数とし， \[ \Lambda(t):=\int^t_0\lambda(\phi(s,x))ds \] で表す．

2.2.3 境界

\(M_v\subset\mathbb{R}^{d(v)}\) の境界のうち，積分曲線が到達し得る部分を \[ \partial^*M_v:=\left\{y\in\partial M_v\,\middle|\,\exists_{t>0}\;\exists_{x\in M_v}\;\phi_v(x,t)=y\right\} \] で表し，その直和を \(\Gamma^*:=\bigcup_{v\in K}\partial^*M_v\times\{v\}\) とする．

さらにそのうち，正な確率を持って到達可能な部分を \[ \Gamma:=\left\{(y,v)\in\Gamma^*\mid\lim_{x\to y}\operatorname{P}_{(x,v)}[T_1=t^*(x,v)]=1\right\} \]

点 \((x,v)\in E\) からフローに従って運動した場合の境界 \(\Gamma^*\) への到達時刻を \[ t^*(x,v):=\inf\left\{t\ge 0\mid\phi_v(x,t)\in\partial^*M_v\right\} \] とする．

2.3 定義

定義 (PDMP: Piecewise Deterministic Markov Process)10

\((\phi(t,-))\) をフロー，\(\lambda:E\to\mathbb{R}_+\) をレート関数，\(Q:(E\cup\Gamma^*)\times\mathcal{E}\to[0,1]\) を確率核とする．

３組 \((\phi,\lambda,Q)\) が定める PDMP \(Z\) とは，次のように帰納的に定まる確率過程をいう：

  1. あるスタート地点 \(z\in E\) から，一様乱数 \(U_1\sim\mathrm{U}([0,1])\) とフロー \(\phi\) が定めるジャンプ時刻 \[ T_1:=\inf\left\{t\ge 0\,\middle|\,e^{-\Lambda(t)}\le U_1\right\}\lor t^*(x,v) \] まで，フローの通りに移動するとする： \[ Z_t=\phi(t,x),\quad t< T_1, \]
  2. 時刻 \(T_1\) での位置は確率核 \(Q\) に従って決まるとする： \[ Z_{T_1}\sim Q\biggr((\phi(T_1,x),v),-\biggl). \]
  3. この手続きを，\(Z_{T_1}\) を次のスタート地点として繰り返すのが特性量 \((\phi,\lambda,Q)\) が定める PDMP \(Z\) である．

\(\lambda\) は局所有限であるから，有限時区間上でのジャンプ数は有限になる．境界への到達によるジャンプ数も同様の性質を満たすと仮定する．11

こうして得た PDMP は，Feller-Dynkin 過程であるとは限らないにも拘らず，時間的に一様な強 Markov 過程である．12

2.4 PDMP の拡張生成作用素

PDMP の拡張生成作用素は，ちょうど（有限な Lévy 測度を持つ）純粋跳躍過程の生成作用素 1 にドリフト項を加えたものになる．

\(X\) を PDMP とする．このとき，\(X\) の拡張生成作用素 \(\widehat{A}\) について，可測関数 \(f\in\mathcal{L}(E\cup\Gamma)\) が \(f\in\mathcal{D}(\widehat{A})\) であるとは，次と同値：

  1. 絶対連続性：任意の \(x\in E\) について， \[ t\mapsto f(\phi(t,x)) \] は \[ s_*(x):=\inf\left\{t\ge 0\mid\operatorname{P}_x[T_1>t]=0\right\} \] に関して，\(t\in[0,s_*(x))\) 上で絶対連続である．
  2. 境界条件： \[ \Gamma:=\left\{x\in E\mid\exists_{t>0}\;\exists_{z\in E}\;x=\phi(t,z)\right\} \] 上では \[ f(x)=\int_Ef(y)Q(x,dy),\quad x\in\Gamma \] を満たす．
  3. 跳躍の局所可積分性：\(E\times\mathbb{R}_+\) 上のランダム関数 \[ (x,s)\mapsto Bf(x,s):=f(x_s)-f(X_{s-}) \] は，\(X\) の跳躍が定める \(E\times\mathbb{R}_+\) 上の Poisson 点過程 \(\eta\) について，ある発散する停止時の増加列 \(\tau_n\) が存在して， \[ \operatorname{E}\left[\int_0^{\tau_n}\int_E\lvert Bf(x,s)\rvert\eta(dxds)\right]<\infty \] \[ n=1,2,\cdots \] が成り立つ．
命題（PDMP の拡張生成作用素の特徴付け）14

\(X\) を PDMP とする．このとき，\(X\) の拡張生成作用素 \(\widehat{A}\) について，\(\mathcal{D}(\widehat{A})\) の全域で次が成り立つ： \[ \widehat{A}f(x)=\mathfrak{X}f(x)+\lambda(x)\int_E\biggr(f(y)-f(x)\biggl)Q(x,dy). \]

\(\mathfrak{X}\) の標準座標 \(\frac{\partial }{\partial x_1},\cdots,\frac{\partial }{\partial x_d}\) に関する成分表示を \(g^i\) とすると，\(\mathfrak{X}f(x)\) は \((g(x)|\nabla f(x))\) とも表せる．

この \(\nabla f\) は，\(f\in C^1(M_v)\) である場合は問題がないが，一般に単に絶対連続である場合は問題になる可能性がある．

しかしそれでも，次を満たす関数 \(\nabla f\) として用意することができる： \[ f(\phi_v(x,t),v)-f(x,v)=\int^t_0\nabla f(\phi_v(x,s),v)ds. \]

2.5 PDMP が Feller-Dynkin 過程であるとき

命題（PDMP の Feller 性の特徴付け）16

\(X\) を PDMP とする．次の全てが成り立つならば，\(X\) は Feller-Dynkin 過程である：

  1. 任意の \(x\in E\) について， \[ t_*(x):=\inf\left\{t>0\mid\phi(t,x)\in\partial E\right\}=\infty. \]
  2. \(\lambda\in C_b(E)\)．
  3. \(Q\) も Feller 性を持つ：\(Q(C_b(E))\subset C(E)\)．

しかし，PDMP は常に (Meyer, 1966) にいう right process ではある．すなわち，次の４条件を満たす：

  1. \(E\) は Lusin 空間である．
  2. \(P_t(\mathcal{L}_b(E))\subset\mathcal{L}_b(E)\)．
  3. \(X\) は殆ど確実に右連続な見本道を持つ．
  4. \(f\) が \(\alpha\)-excessive function ならば，\(t\mapsto f\circ X_t\) も殆ど確実に右連続な見本道を持つ．

2.6 PDMP のシミュレーション

2.6.1 到着時刻のシミュレーション

レート関数が \(\lambda>0\) であるとき，Poisson イベントのシミュレーションは次のように，指数分布に従う確率変数のシミュレーションによって行える．

累積関数を \[ \Lambda(t):=\int^t_0\lambda(s)ds \] とすると，\(N_t\sim\mathrm{Pois}(\Lambda(t))\) であるから，最初のイベントの到着時刻 \(T_1\) は，次の生存関数によって特徴付けられる： \[ \operatorname{P}[T_1\ge t]=\operatorname{P}[N_t=0]=\exp\left(-\Lambda(t)\right). \]

ここで，\(E\sim\operatorname{Exp}(1)\) とすると， \[ \operatorname{P}[\Lambda^{-1}(E)\ge t]=\operatorname{P}\left[E\ge\Lambda(t)\right]=e^{-\Lambda(t)}. \] 従って，\(T_1\overset{\text{d}}{=}f^{-1}(E)\) である．

すなわち，\(E\) を通じて \[ T_1':=\inf\left\{t\ge 0\,\middle|\,\int^t_0\lambda(s)ds=E\right\} \] を計算すれば，最初の到着時刻が計算できる．17

2.6.2 PDMP のためのパッケージ

Joris Bierkens ら開発の R パッケージ RZigZag (GitHub / CRAN) を通じて実行してみる．

install.packages("Rcpp")
install.packages("RcppEigen")
install.packages("RZigZag")

3 Zig-Zag Sampler (Bierkens et al., 2019)

3.1 導入

１次元の Zig-Zag 過程は元々，Curie-Weiss 模型 における Glauber 動力学を lifting により非可逆化して得る Markov 連鎖の，スケーリング極限として特定された Feller-Dynkin 過程である (Bierkens and Roberts, 2017)．

ただし，(Goldstein, 1951) で電信方程式と関連して，同様の過程が扱われた歴史もある．

3.2 設定

Zig-Zag 過程 \(Z=(X,\Theta)\) の状態空間は \(E=\mathbb{R}^d\times\{\pm1\}^d\) と見ることが多い．

\(\theta\in\{\pm1\}^d\) は速度を表す．すなわち，全座標系と \(45\) 度をなす方向に，常に一定の単位速度で \(\mathbb{R}^d\) 上を運動するとする．

換言すれば，決定論的なフローは次のように定める： \[ \frac{d \phi_{(x,\theta)}(t)}{d t}=\theta, \] \[ \frac{d \Theta}{d t}=0,\qquad\Theta_0=0. \]

3.3 アルゴリズム

後述の関数 \(\lambda\) に対して，各座標 \(i\in[d]\) におけるレートを \[ m_i(t):=\lambda(x+\theta t,\theta) \] で定めた \(\mathbb{R}^d\) 上の Poisson 過程を考える．

多次元の Poisson 過程の各成分の跳躍は独立だから，18 それぞれの成分ごとに Poisson 到着時刻 \(T_i\;(i\in[d])\) をシミュレーションし，最初に到着したものを \(T_j\) とすると，関数 \[ F_j(\theta)_i=\begin{cases}-\theta_i&i=j\\\theta_i&i\ne j\end{cases} \] に従ってジャンプすると考えて良い．

状態空間は \[ E=\bigcup_{\theta\in\{\pm1\}^d}\mathbb{R}^d\times\{\theta\} \] とみなすところから始める必要がある．

加えて，レート関数は \[ \lambda(x,\theta):=\sum_{i=1}^d\lambda_i(x,\theta) \] であり，跳躍測度は点過程 \[ Q((x,\theta),-):=\sum_{i=1}^d\frac{\lambda_i(x,\theta)}{\lambda(x,\theta)}\delta_{(x,F_i(\theta))}(-) \] とみなした場合に当たる．

3.4 レート関数

PDMP の妙は全てレート関数に宿っている．

レート \(\lambda:E\to\mathbb{R}_+\) は，負の対数密度 \(U\in C^1(\mathbb{R}^d)\) が定める目標分布 \(\pi(dx)\,\propto\,e^{-U(x)}dx\) に対して， \[ \lambda_i(x,\theta):=(\theta_i\partial_iU(x))_++\gamma_i(x,\theta_{-i})\quad(i\in[d]) \] と定める．ただし，\(\gamma_i\) は，\(\theta_i\) のみには依らない任意の連続関数 \(\gamma_i:E\to\mathbb{R}_+\) とした．19

(Bierkens et al., 2019, pp. 1294 定理2.2)

\(U\in C^1(\mathbb{R}^d)\) は \(e^{-U}\in\mathcal{L}^1(\mathbb{R}^d)\) を満たすとする．このとき，\(Z\) は次の分布 \(\mu=\pi\otimes\mathrm{U}(\{\pm1\}^{d})\in\mathcal{P}(E)\) を不変にする： \[ \mu(dxd\theta)=\frac{1}{2^d\mathcal{Z}}e^{-U(x)}\,dxd\theta \]

3.5 部分サンプリング

MCMC の計算複雑性のボトルネックは，尤度の評価にある．各ステップで全てのデータを用いて尤度を計算する必要がある点が，MCMC を深層学習などの大規模データの設定への応用を難しくしている (Murphy, 2023, p. 647)．

\(p(x)\) を事前分布，\(p(y|x)\) を観測のモデル（または尤度）とし，データ \(y_1,\cdots,y_n\) は互いに独立であるとする．このとき，事後分布 \(\pi(x):=p(x|y)\) は \[ \pi(x)\,\propto\,\left(\prod_{k=1}^n p(y_k|x)\right)p(x) \] より，Hamiltonian \(U\) は \[\begin{align*} U(x)&=-\sum_{k=1}^n\log p(y_k|x)-\log p(x)\\ &=\frac{1}{n}\sum_{k=1}^n\biggr(-n\log p(y_k|x)-\log p(x)\biggl). \end{align*}\] と表せる．

すなわち，各微分係数 \(\partial_i U(x)\) は，独立な観測 \(y_1,\cdots,y_n\) が定める統計量 \[ E^i_k(x):=\frac{\partial }{\partial x_i}\biggr(-n\log p(y_k|x)-\log p(x)\biggl) \] の平均により推定される値となっている．

よって，精度は劣るかもしれないが，一様に選んだ \(K\sim\mathrm{U}([n])\) から定まる \(E^i_K\) の値は \(\partial_i U(x)\) の不偏推定量となっている．20

これは，元々のレート関数 \[ \lambda_i(x,\theta)=\frac{1}{n}\sum_{k=1}^n(\theta E^i_k(x))_+ \] に対して， \[ \gamma_i(x,\theta):=\frac{1}{n}\sum_{k=1}^n(\theta_iE^i_k(x))_+-\left(\theta_i\frac{1}{n}\sum_{k=1}^nE^i_k(x)\right)_+ \] という項を加えて得る Zig-Zag サンプラーとみなせるためである．

こうして定義された \(\gamma_i\) が非負関数である限り，平衡分布は \(\mu\) のままであるが，これは関数 \((-)_+\) の凸性から従う．

こうして，サブサンプリングの実行による精度の劣化が，(Andrieu and Livingstone, 2021) の枠組みで捉えられる，ということでもある．

加えて，レート関数が大きくなっており，したがってそもそも尤度関数の評価の回数が増えていることにも注意．

3.6 制御変数による分散低減

上述の \(\partial_iU(x)\) の不偏推定量の分散は，制御変数の方法を用いて低減できる．これにより，事前処理の部分を除けば，データのサイズに依存しない計算複雑性で事後分布からの正確なサンプリングが可能になる．

3.7 シミュレーション

Code 
library(RZigZag)
library(ggplot2)
V <- matrix(c(3,1,1,3),nrow=2)
mu <- c(2,2)
result <- ZigZagGaussian(V, mu, 100)
ggplot() +
   geom_path(aes(x=result$Positions[1,], y=result$Positions[2,]), color="#78C2AD") +
   geom_point(aes(x=result$Positions[1,], y=result$Positions[2,]), color="#78C2AD") +
   labs(x="", y="", title="Zig-Zag Sampler") +
   theme_void() +
   theme(text=element_text(size=12), axis.title=element_text(color="#78C2AD"), plot.title=element_text(color="#78C2AD"))

\(\mathrm{N}_2\left(\begin{pmatrix}2\\2\end{pmatrix},\begin{pmatrix}3&1\\1&3\end{pmatrix}\right)\) に対する Zig-Zag 過程
Code 
set.seed(123)
dim <- 2
dof <- 1
result <- ZigZagStudentT(dof, dim, n_iter=1000, sphericallySymmetric = TRUE)
ggplot() +
   geom_path(aes(x=result$Positions[1,], y=result$Positions[2,]), color="#78C2AD") +
   geom_point(aes(x=result$Positions[1,], y=result$Positions[2,]), color="#78C2AD") +
   labs(x="", y="", title="Zig-Zag Sampler") +
   theme_void() +
   theme(text=element_text(size=12), axis.title=element_text(color="#78C2AD"), plot.title=element_text(color="#78C2AD"))

Cauchy 分布 \(\mathrm{C}(0,1)=\mathrm{t}(1)\) に対する Zig-Zag 過程

続きは次の稿も参照：

Article Image

Zig-Zag サンプラー

ジャンプと確定的な動きによる新たな MCMC 手法

4 Bouncy Particle Sampler (Bouchard-Côté et al., 2018)

4.0.1 部分サンプリング

Zig-Zag サンプラーに対応するサブサンプリングの技術を (Pakman et al., 2017) が提案している．

4.0.2 シミュレーション

Code 
V <- matrix(c(3,1,1,3),nrow=2)
mu <- c(2,2)
x0 <- c(0,0)
result <- BPSGaussian(V, mu, n_iter = 100, x0 = x0)
ggplot() +
   geom_path(aes(x=result$Positions[1,], y=result$Positions[2,]), color="#78C2AD") +
   geom_point(aes(x=result$Positions[1,], y=result$Positions[2,]), color="#78C2AD") +
   labs(x="", y="", title="Bouncy Particle Sampler") +
   theme_void() +
   theme(text=element_text(size=12), axis.title=element_text(color="#78C2AD"), plot.title=element_text(color="#78C2AD"))

\(\mathrm{N}_2\left(\begin{pmatrix}2\\2\end{pmatrix},\begin{pmatrix}3&1\\1&3\end{pmatrix}\right)\) に対する Zig-Zag 過程

References

Andrieu, C., and Livingstone, S. (2021). Peskun–Tierney ordering for Markovian Monte Carlo: Beyond the reversible scenario. The Annals of Statistics, 49(4), 1958–1981.
Bierkens, J., Fearnhead, P., and Roberts, G. (2019). The Zig-Zag Process and Super-Efficient Sampling for Bayesian Analysis of Big Data. The Annals of Statistics, 47(3), 1288–1320.
Bierkens, J., and Roberts, G. (2017). A Piecewise Deterministic Scaling Limit of Lifted Metropolis-Hastings in the Curie-Weiss Model. The Annals of Applied Probability, 27(2), 846–882.
Bouchard-Côté, A., Vollmer, S. J., and Doucet, A. (2018). The bouncy particle sampler: A nonreversible rejection-free markov chain monte carlo method. Journal of the American Statistical Association, 113(522), 855–867.
Davis, M. H. A. (1984). Piecewise-deterministic markov processes: A general class of non-diffusion stochastic models. Journal of the Royal Statistical Society. Series B (Methodological), 46(3), 353–388.
Davis, M. H. A. (1993). Markov models and optimization,Vol. 49. Chapman & Hall.
Ethier, S. N., and Kurtz, T. G. (1986). Markov processes: Characterization and convergence. John Wiley & Sons, Inc.
Goldstein, S. (1951). On Diffusion by Discontinuous Movements, and on the Telegraph Equation. The Quarterly Journal of Mechanics and Applied Mathematics, 4(2), 129–156.
Meyer, P. A. (1966). Probability and potentials. Blaisdell Publishing Company.
Murphy, K. P. (2023). Probabilistic machine learning: Advanced topics. MIT Press.
Pakman, A., Gilboa, D., Carlson, D., and Paninski, L. (2017). Stochastic bouncy particle sampler. In D. Precup and Y. W. Teh, editors, Proceedings of the 34th international conference on machine learning,Vol. 70, pages 2741–2750. PMLR.
Pazy, A. (1983). Semigroups of linear operators and applications to partial differential equations,Vol. 44. Springer New York.
Revuz, D., and Yor, M. (1999). Continuous martingales and brownian motion. Springer Berlin, Heidelberg.
Sen, D., Sachs, M., Lu, J., and Dunson, D. B. (2020). Efficient posterior sampling for high-dimensional imbalanced logistic regression. Biometrika, 107(4), 1005–1012.
Vasdekis, Georgios. (2021). On zig-zag extensions and related ergodicity properties (PhD thesis). University of Warwick. Retrieved from http://webcat.warwick.ac.uk/record=b3714913
Vasdekis, G., and Roberts, G. O. (2023). Speed up Zig-Zag. The Annals of Applied Probability, 33(6A), 4693–4746.

Footnotes

  1. (Pazy, 1983, p. 2) 参照．↩︎

  2. この事実は第 1.3 節で２つ目の構成と同時に証明される．↩︎

  3. (Ethier and Kurtz, 1986, pp. 163–164) 参照．↩︎

  4. \(f(y+x)-f(x)\) ではなくて \(f(y)-f(x)\) であるので，\(\mu(x,dy)\) は必ずしも現在地点 \(x\) からみた変位の分布ではないことに注意．↩︎

  5. ただし，ジャンプが頻繁すぎないことも必要である．(Davis, 1993, p. 60) 仮定24.4 では，有界区間上のジャンプは有限回であると過程している：\(\operatorname{E}[N_t]<\infty\;(t\in\mathbb{R}_+)\)．すなわち，A 型までの加法過程を PDMP というのであって，B 型と C 型は，シミュレーションが容易な連続時間過程であるという美点を逃してしまう．加法過程の分類については，この 稿 も参照．↩︎

  6. 当時は連続時間確率過程といえば拡散過程であり，SDE によるモデリングが興隆した時代であった．(Davis, 1993) では，必ずしも SDE を使うことが自然なモデリング方法でないにも拘らず，無理やり SDE の枠組みに落とし込もうとする当時の慣行を批判し，PDMP はこのギャップを埋めるために開発した，としている．なお，(Davis, 1993) では PDMP ではなく PDP と呼んでいる．↩︎

  7. (Davis, 1993, p. 57) では \(K\) は可算という仮定をおいている．↩︎

  8. 積空間としての \(\sigma\)-代数を導入する．↩︎

  9. (Davis, 1993) では，\(\mathfrak{X}_v\) は局所 Lipcthiz 連続である上に，（暗黙のうちに）完備で，\(\phi_v\) は任意の \(t\in\mathbb{R}\) について定義されるとしていることに注意．(G. Vasdekis and Roberts, 2023) では違う．↩︎

  10. (Davis, 1993, p. 58) と (Georgios Vasdekis, 2021, pp. 15–16) に倣った．↩︎

  11. このことを (Davis, 1993, p. 60) 仮定24.4としている．↩︎

  12. (Davis, 1984)，(Davis, 1993, p. 64) 定理25.5も参照．↩︎

  13. (Davis, 1993, p. 69) 定理26.14．↩︎

  14. (Davis, 1993, p. 69) 定理26.14 の後半．↩︎

  15. (Georgios Vasdekis, 2021, p. 18) 系2.3.4 も参照．↩︎

  16. (Davis, 1993, p. 77) 定理27.6．↩︎

  17. 簡単な対象分布では，\(\Lambda(t)=E\) の解が解析的に求まる事が多い．これが intractable である場合は，剪定 を用い，\(\lambda^*\) としては affine 関数を用いる事が多いが，一般に３次の多項式までならば解の公式があるために，\(\Lambda(t)=E\) を解析的に解く方法から効率的にシミュレーションできる．どれくらい \(\lambda^*\) として \(\lambda\) に近いものを選べば良いかは不明．(Georgios Vasdekis, 2021, p. 13) 命題2.2.2も参照．↩︎

  18. (Revuz and Yor, 1999, p. 473) 命題XII.1.7．↩︎

  19. 従って，レート関数 \(\lambda\) は連続とする．この関数 \(\gamma_i\) は，\(U\) の情報には依らない追加のリフレッシュ動作を仮定に加える．実際，\(\lambda_i(x,\theta)-\lambda_i(x,F_i(\theta))=\theta_i\partial_iU(x)\) である限り，\(\theta\) と \(F_i(\theta)\) の往来には影響を与えず釣り合っているため，どのような \(\gamma_i\) をとっても，平衡分布には影響を与えない．しかし，高くするごとにアルゴリズムの対称性が上がるため，\(\gamma\equiv0\) とすることが Monte Carlo 推定量の漸近分散を最小にするという (Andrieu and Livingstone, 2021)．↩︎

  20. このとき，必ずしも \(K\sim\mathrm{U}([n])\) とする必要はなく，特定の観測に重みをおいても良い (Sen et al., 2020)．↩︎