【R】エーレンフェストモデル

1.エーレンフェストモデル

エーレンフェストモデル(Ehrenfest model)は, 外界から閉ざされた容器内での気体の拡散を説明する場面に使われる確率モデルの一つで、マルコフ連鎖の例としてしばしば取り上げられます.

以下、Rでのシミュレーションと合わせてこのモデルの基本性質について取り上げます.

1.1 定式化

2つの箱(それぞれA, Bとします)内に合計 $r$ 個のボールが入っており, ボールには $1$ ~ $r$ までの番号が書かれているとします. 次のルールに従ってボールの遷移を考えます.

ルール: $r$ 個のボールから1つをランダムに選び, もう他方の箱へ移す.

この遷移が $t$ 回行われたときの箱A内のボールの個数 $X(t)$ は、 $S = \{ 0, 1, \ldots, r \}$ 上のマルコフ連鎖になります.

この遷移確率は次で与えられることがわかります,

$p_{ij}:=p(X(t+1)=j \ | \ X(t)= i) = \begin{cases} \frac{i}{r} \quad (j=i-1) \\ \frac{r-i}{r} \quad (j=i+1) \\ 0 \quad (\text{otherwise}). \end{cases}$

1.2 定常分布

十分に長い時間が経過したとき, $X(t)$ の分布は定常分布 $\pi$ に近づきます. $\pi_i = \pi(X=i)$ としたとき, 定常分布(equilibrium distribution)は次の詳細釣り合いの式から求めることが可能です,

$\pi_i p_{ij} = \pi_j p_{ji} \quad (0 \leq i,j \leq r)$

特に $j=i-1$ とすることで次が得られます,

$\pi_i = \frac{p_{i-1i}}{p_{ii-1}}\pi_{i-1} = \frac{r-i+1}{i}\pi_{i-1}=\cdots = \binom{r}{i}\pi_0.$

分布の条件 $1 = \pi_0+\cdots+\pi_r$ から $\pi_0 = 1 / 2^{r}$ ともとまる.

エーレンフェストモデルの定常分布 エーレンフェストモデルの定常分布は二項分布 $\begin{equation}　{\sf Binom}\bigl( r, \frac{1}{2} \bigr)　\end{equation}$ である.

1.3 Rを使った定常分布のシミュレーション

require(tidyverse)
ehrenfest.path <- function(r, len, init) {
    path <- c(init)
    for (t in 1:(len-1)) {
        x_next <- ifelse(rbinom(1, 1, path[t]/r) == 1, path[t]-1, path[t]+1)
        path <- c(path, x_next)
    }
    return(path)
}

# 長さ4000のパス
sample.path <- ehrenfest.path(r=100, len=4000, init=0)

data.frame(value=sample.path) %>%
    ggplot(aes(x=value)) + 
    geom_density(fill="blue", alpha=.4) +
    geom_col(aes(x=x, y=y), data.frame(x=c(0:100), y=dbinom(c(0:100), size=100, prob=.5)), alpha=.7) +
    xlim(c(0,100)) +
    ggtitle('histgram of path v.s. equilibrium distribution')

f:id:kimigayoseishou:20200427200246p:plain — サンプルの密度分布(薄青)と, 定常分布(試行回数=100, 成功確率=1/2の二項分布)の比較

2.1 ギブスエントロピー

有限個の状態をとりうる系について, ギブスエントロピー(Gibbs entropy)をつぎの式で定義されます.

$H = - \sum_i p_i \ln p_i$

熱力学第二法則によれば, 不可逆な変化のもとでギブスエントロピーは増大します. エーレンフェストモデルでは箱A内のボール個数の分布は定常分布へ近づくことになり, 不可逆な変化が発生します. 以下では, エーレンフェストモデルのギブスエントロピーが時間と共に増加することをシミュレーションにより観察します.

2.1 エーレンフェストモデルのシミュレーション

エーレンフェストモデルはマルコフ連鎖に従うため, $X(t)$ の確率分布から $X(t+1)$ の確率分布を次のように明示的に表すことができます,

$p(X(t+1) = i) \\ = \sum_{j=0}^r p(X(t+1)=i,\ X(t) = j) \\ = \sum_{j=0}^r p(X(t+1)=i | \ X(t) = j) p(X(t) = j) \\ = p(X(t)=i-1)p_{i-1i} + p(X(t)=i+1)p_{i+1i} \\ = \frac{r-(i-1)}{r}p(X(t)=i-1) + \frac{i+1}{r}p(X(t)=i+1)$

この式による, エーレンフェストモデルの時間変化のシミュレーションのRコード以下のようになります.

#' @param T_max 
#' @param N 時刻を刻む個数
#' @param r ボールの個数
#' @param init r-vetor 初期値
Ehrenfest.sim <- function(T_max, N, r, init) {
    # 横ベクトルの要素を左にずらす
    shift.p <- function(x) {
        return(c(x, 0)[c(-1)])
    }
    # 横ベクトルの要素を右にずらす
    shift.m <- function(x) {
        return(c(0, x)[c(-length(x)-1)])
    }

    sim <- matrix(0, nrow = T_max, ncol = r+1)
    sim[1,] <- init
    for (i in 2:T_max) {
        sim[i,] <- (1 + 1 / r - c(0:r)/r) * shift.m(sim[i-1,]) + (1 / r + c(0:r) / r) * shift.p(sim[i-1,])
    }
    return(sim)
}

ehrenfest.sim <- Ehrenfest.sim(T_max=200, N=100, r=100, init=c(c(1, 1)/2, rep(0, 99)))

# シミュレーション結果をプロット
plt.sim <- function(sim, time) {
    x.range <- c(0:(ncol(sim)-1))
    plot(x=x.range, y=sim[1,], xlim = c(0, max(x.range)), ylim = c(0, 0.6), type="l", 
         ylab = 'probability')
    
    for (t in time) {
        lines(x=x.range, y=sim[t,], lty="solid")
    }
}

# 10回の遷移づつ確率分布をプロット
# (横軸 - ボールの個数, 縦軸 - 確率)
plt.sim(ehrenfest.sim, c(1:10) * 20)

f:id:kimigayoseishou:20200427200154p:plain — 確率質量分布の時間変化

プロットの結果より, 左端に局在していた確率分布が時間の経過と共に定常分布に近づく様子が確認できます.

次にギブスエントロピーのシミュレーションを行います.

gibbs_entropy <- function(x) - sum(x * log(x), na.rm = TRUE)
sim.entropy <- apply(ehrenfest.sim, MARGIN = 1, gibbs_entropy)
plot(sim.entropy, type="l", xlab = 'time', ylab='gibbs entropy')