Simulation Theory

  1. 1. 確率分布
    1. 1.1. (0-1)分布
    2. 1.2. 伯努利试验,二项分布
      1. 1.2.1. 例題
    3. 1.3. 泊松分布
      1. 1.3.1. 例題
    4. 1.4. 几何分布
      1. 1.4.1. 例題
    5. 1.5. 均匀分布
      1. 1.5.1. 例題
    6. 1.6. 指数分布
    7. 1.7. 正态分布
    8. 1.8. 标准正态分布
  2. 2. マルコフ連鎖(Markov chain)
    1. 2.1. 確率過程
    2. 2.2. マルコフ連鎖の例
    3. 2.3. マルコフ連鎖の定義
    4. 2.4. 2状態マルコフ連鎖
    5. 2.5. マルコフ連鎖の時間ダイナミクス
    6. 2.6. 例題2.1
    7. 2.7. 状態の分類
    8. 2.8. 例題2.2
    9. 2.9. 例題2.3
    10. 2.10. 例題2.4
    11. 2.11. 例題2.5
    12. 2.12. マルコフ連鎖の定常分布
    13. 2.13. 例題2.6
    14. 2.14. 吸収マルコフ連鎖
    15. 2.15. 例題2.7
    16. 2.16. 1次元ランダム・ウォーク
    17. 2.17. 例題2.8
    18. 2.18. 例題2.9
  3. 3. ポアソン過程
    1. 3.1. ポアソン過程の定義
    2. 3.2. ポアソン過程の性質
    3. 3.3. ポアソン分布
    4. 3.4. 例題3.1
    5. 3.5. 事象の発生時間間隔の分布
    6. 3.6. 例題3.2
    7. 3.7. n番目の事象が起こるまで時間 SnS_nSn​ が従う分布
    8. 3.8. 例題3.3
    9. 3.9. 例題3.4
    10. 3.10. 純出生過程
    11. 3.11. ユール・ファーリ過程
    12. 3.12. 純死滅過程
    13. 3.13. 例題3.5
    14. 3.14. 出生死滅過程
    15. 3.15. 例題3.6
    16. 3.16. 例題3.7
  4. 4. 待ち行列
    1. 4.1. 待ち行列理論の発展
    2. 4.2. 待ち行列モデル
    3. 4.3. ケンドールの記号
      1. 4.3.1. 客の到着(A)/サービス時間(B)
    4. 4.4. 単一待ち行列モデル
      1. 4.4.1. M/M/1(INF)
      2. 4.4.2. リトルの公式
      3. 4.4.3. 待ち行列で用いられる特性値
      4. 4.4.4. M/M/1(N)
    5. 4.5. 例題4.1
    6. 4.6. 複数待ち行列モデル
      1. 4.6.1. M/M/s(INF)
      2. 4.6.2. M/M/s(N)
    7. 4.7. 例題4.2
    8. 4.8. 例題4.3
    9. 4.9. 例題4.4
  5. 5. 乱数
  6. 6. エージェントモデルシミュレーション
  7. 7. 参考

CS专业课学习笔记

確率分布

(0-1)分布

定义:若X的概率分布律为

XX 00 11
PP 1p1-p pp

其中0<p<10 < p < 1,就称X服从参数为p的0-1分布(或两点分布),记为X01(p)\color{blue}{X \sim 0-1(p)}XB(1,p)\color{blue}{X \sim B(1, p)}

其分布律还可以写成:P(X=k)=pk(1p)1k,k=0,1P(X = k) = p^k(1 - p)^{1-k}, \quad k = 0,1

伯努利试验,二项分布

设试验E只有两个可能的结果:AAA\overline{A},且P(A)=p,0<p<1P(A) = p, \quad 0 < p < 1。将E独立地重复进行n次,则称这一串重复地独立试验为n重伯努利试验

二项分布定义:若XX地的概率分布律为:

P(X=k)=Cnkpk(1p)nk,k=0,1,,n\begin{aligned} \color{blue}{P(X = k) = C_n^k p^k (1 - p)^{n-k}, \quad k = 0, 1, \cdots, n} \end{aligned}

其中n1,0<p<1n \geq 1,\quad 0 < p < 1,就称XX服从参数为n,pn,p二项分布(Binomial),记为XB(n,p)\color{blue}{X \sim B(n, p)}

  • E(X)=npE(X) = np
  • V(X)=np(1p)V(X) = np(1 - p)

例題

〇×で解答する試験において、気まぐれに解答する。10問ある時、正答の個数を XX として、XX の平均と分散を求めよ。また、X8X \geq 8 となる確率を求めよ

  • E(X)=np=100.5=5E(X) = np = 10 * 0.5 = 5
  • V(X)=np(1p)=100.50.5=2.5V(X) = np(1 - p) = 10 * 0.5 * 0.5 = 2.5
  • P(X8)=P(X=8)+P(X=9)+P(X=10)=C108(12)8(12)2+C109(12)9(12)+C1010(12)10(12)00.0547P(X \geq 8) = P(X = 8) + P(X = 9) + P(X = 10) = C_{10}^8 (\frac{1}{2})^8 (\frac{1}{2})^2 + C_{10}^9 (\frac{1}{2})^9 (\frac{1}{2}) + C_{10}^{10} (\frac{1}{2})^{10} (\frac{1}{2})^0 \simeq 0.0547

泊松分布

定义:若XX的概率分布律为

P(X=k)=λkeλk!,k=0,1,2,\begin{aligned} \color{blue}{P(X = k) = \frac{\lambda^k e^{-\lambda}}{k!}, \quad k = 0, 1, 2, \cdots} \end{aligned}

  • E(X)=λE(X) = \lambda
  • V(X)=λV(X) = \lambda

其中λ>0\lambda > 0,就称 XX 服从参数为 λ\lambda泊松分布(Poisson),记为Xπ(λ)\color{blue}{X \sim \pi(\lambda)}XP(λ)\color{blue}{X \sim P(\lambda)}

泊松分布的用途:

  • 某人一天内收到的微信的数量
  • 来到某公交汽车站的乘客
  • 某放射性物质发射出的粒子
  • 显微镜下某区域中的白血球

如果某事件以固定强度 λ\lambda,随机且独立地出现,该事件在单位时间内出现地次数(个数)可以看成是服从泊松分布

二项分布与泊松分布有以下近似关系:

Cnkpk(1p)nkeλλkk!,λ=np,n>10,p<0.1\begin{aligned} C_n^k p^k (1 - p)^{n - k} \simeq \frac{e^{-\lambda}\lambda^k}{k!}, \quad \lambda = np, n > 10, p < 0.1 \end{aligned}

例題

ある製品の製造工程では、微小領域に1万個の点にレーザーを当てていく。位置をずれるのは2万個に一つであるという。この工程ではずれが一つもないという確率を求めよ

  • λ=np=10000120000=0.5\lambda = np = 10000 * \frac{1}{20000} = 0.5
  • E(X)=λ=0.5E(X) = \lambda = 0.5
  • V(X)=λ=0.5V(X) = \lambda = 0.5
  • P(X=0)=P(0,λ)=e0.50.500!0.606531P(X = 0) = P(0, \lambda) = e^{-0.5} \cdot \frac{0.5^0}{0!} \simeq 0.606531
  • B(n,p)=C100000(120000)0(1999920000)100000.606523B(n, p) = C_{10000}^0 (\frac{1}{20000})^0 (\frac{19999}{20000})^{10000} \simeq 0.606523

几何分布

定义:若 XX 的概率分布律为

P(X=k)=p(1p)k1,k=1,2,3,\begin{aligned} \color{blue}{P(X = k) = p(1 - p)^{k - 1}, \quad k = 1, 2, 3, \cdots} \end{aligned}

其中 0<p<10 < p < 1,称 XX 服从参数为 pp几何分布(Geometric),记为 XGeom(p)\color{blue}{X \sim Geom(p)}

几何分布的用途:在重复多次的伯努利实验中,实验进行到某种结果出现第一次为止,此时的试验总次数服从几何分布。如:射击,首次击中目标时射击的次数

  • E(X)=1pE(X) = \frac{1}{p}
  • V(X)=1pp2V(X) = \frac{1 - p}{p^2}

例題

合格率10%の試験を挑み、XX 回目で初めて合格する。XX の平均と分散を求めよ。また、X=3X = 3 となる確率を求めよ

  • E(X)=1p=10.1=10E(X) = \frac{1}{p} = \frac{1}{0.1} = 10
  • V(X)=1pp2=10.10.12=90V(X) = \frac{1 - p}{p^2} = \frac{1 - 0.1}{0.1^2} = 90
  • P(X=3)=0.920.1=0.081P(X = 3) = 0.9^2 * 0.1 = 0.081

均匀分布

XX 的概率密度函数为

f(x)={1ba,x(a,b)0,其他\begin{aligned} \color{blue}{f(x) = \begin{cases} \frac{1}{b-a}, \quad &x \in (a, b) \\ 0, \quad &\text{其他} \end{cases}} \end{aligned}

其中 a<ba < b,就称 XX 服从 (a,b)(a, b) 上的均匀分布(Uniform),记为 XU(a,b)\color{blue}{X \sim U(a, b)}XUnif(a,b)\color{blue}{X \sim \text{Unif}(a, b)}

均匀分布的性质:均匀分布具有等可能性,即,服从 U(a,b)U(a, b) 上的均匀分布的随机变量 XX 落入 (a,b)(a, b) 中的任意子区间上的概率只与其区间长度有关与区间所处的位置无关。即 XX 落入 (a,b)(a, b) 中的等长度的任意子区间上是等可能的。

  • E(X)=b+a2E(X) = \frac{b + a}{2}
  • V(X)=(ba)212V(X) = \frac{(b - a)^2}{12}

例題

ある自動車のブレーキテストで時速50kmでブレーキをかけた時、止まるまでに要する距離 XX は区間 (a,b)(a, b) 上の一様分布に従うという。bb よりも aa に近いところで止まる確率を求めよ

P(a<X<a+b2)=aa+b21badx=12P(a < X < \frac{a+b}{2}) = \int_a^{\frac{a+b}{2}} \frac{1}{b - a} dx = \frac{1}{2}

指数分布

XX 的概率密度函数为

f(x)={λeλx,x>00,x0\begin{aligned} \color{blue}{f(x) = \begin{cases} \lambda e^{-\lambda x}, \quad &x > 0 \\ 0, \quad &x \leq 0 \end{cases}} \end{aligned}

其中 λ>0\lambda > 0,就称 XX 服从参数为 λ\lambda指数分布(Exponential),记为 XE(λ)\color{blue}{X \sim E(\lambda)}XExp(λ)\color{blue}{X \sim \text{Exp}(\lambda)}

  • E(X)=1λE(X) = \frac{1}{\lambda}
  • V(X)=1λ2V(X) = \frac{1}{\lambda^2}

随机变量 XX分布函数为:

F(x)={1eλx,x>00,x0\begin{aligned} \color{blue}{F(x) = \begin{cases} 1 - e^{-\lambda x}, \quad &x > 0 \\ 0, \quad &x \leq 0 \end{cases}} \end{aligned}

指数分布的性质:指数分布具有无记忆性

指数分布的用途:

  • 指数分布可以用来表示独立随机事件发生的时间间隔,比如旅客进机场的时间间隔,中文维基百科新条目出现的时间间隔等等
  • 在排队论中,一个旅客接受服务的时间长短也可以用指数分布来近似
  • 无记忆性的现象(连续时)

正态分布

XX 的概率密度函数为

f(x)=12πσe(xμ)22σ2,<x<+\begin{aligned} \color{blue}{f(x) = \frac{1}{\sqrt{2\pi}\sigma} e^{-\frac{(x - \mu)^2}{2\sigma^2}}, \quad -\infty < x < +\infty} \end{aligned}

其中 <μ<+,σ>0-\infty < \mu < +\infty, \sigma > 0 ,就称 XX 服从参数为 μ,σ\mu, \sigma正态分布(或高斯分布),记为 XN(μ,σ2)\color{blue}{X \sim N(\mu, \sigma^2)}

  • E(X)=μE(X) = \mu
  • V(X)=σ2V(X) = \sigma^2

正态分布的特征:

  • f(x)f(x) 关于 x=μx = \mu 对称
  • xμx \leq \mu 时,f(x)f(x) 是严格单调递增函数
  • fmax=f(μ)=12πσf_{max} = f(\mu) = \frac{1}{\sqrt{2\pi}\sigma}
  • limxμf(x)=0\lim_{|x - \mu| \rightarrow \infty} f(x) = 0

正态分布的参数:

  • μ\mu位置参数,决定对称轴位置
  • σ\sigma尺度参数,决定曲线分散程度

正态分布的用途:

  • 自然界和人类社会中很多现象可以看做正态分布
    • 例如:人的生理尺寸(身高,体重),医学检验指标(红细胞数,血小板),测量误差等等
  • 多个随机变量的和可以用正态分布来近似
    • 例如:注册MOOC的某位同学完成所有作业的时间,二项分布等等

标准正态分布

ZN(0,1)Z \sim N(0, 1),称 ZZ 服从标准正态分布

  • ZZ 的概率密度函数:ϕ(z)=12πez22\phi(z) = \frac{1}{\sqrt{2\pi}}e^{-\frac{z^2}{2}}
  • ZZ 的分布函数:Φ(z)=Z12πet22dt\Phi(z) = \int_{-\infty}^Z \frac{1}{\sqrt{2\pi}} e^{-\frac{t^2}{2}} \, dt

标准正态分布的分布函数有一个重要的性质:

Φ(z0)=1Φ(z0)\Phi(-z_0) = 1 - \Phi(z_0)

对于任意的实数 z0z_0 都成立

此外,当 XN(μ,σ2)X \sim N(\mu, \sigma^2) 时,XμσN(0,1)\frac{X - \mu}{\sigma} \sim N(0, 1),由此可知,当 XN(μ,σ2)X \sim N(\mu, \sigma^2) 时,对于任意实数 aa,有

FX(a)=P(Xa)=P(Xμσaμσ)=Φ(aμσ)\color{blue}{F_X(a) = P(X \leq a) = P(\frac{X - \mu}{\sigma} \leq \frac{a - \mu}{\sigma}) = \Phi(\frac{a - \mu}{\sigma})}

マルコフ連鎖(Markov chain)

马尔可夫过程是研究信号多级传输,分子的布朗运动,顾客服务,计算机网络流量等诸多问题时使用的经典模型

確率過程

時間 tt に応じて偶然に変化する結果を X(t)X(t) とする

  • X(t)X(t):確率変数
  • {X(t);t0}\lbrace X(t); t \geq 0 \rbrace:確率過程

例:

  • tt 時に届くemailの数
  • tt 日の株価の終値
  • コイン投げ tt 回における表の出現回数

マルコフ過程:ある時点における現象の確率が直前の結果だけで決まる確率過程

  • マルコフ連鎖:離散時間 {t=0,1,2,}\lbrace t = 0, 1, 2, \cdots \rbrace、離散値 {i;i=0,1,2,}\lbrace i;i = 0, 1, 2, \cdots \rbrace をとるマルコフ過程
  • ポアソン過程:連続時間 t0t \geq 0、離散値をとるマルコフ過程

マルコフ連鎖の例

0, 1のラベルがついた球が入っている箱A, 箱Bを考える

初めにどちらかの箱から球を取り出し、玉のラベルを書き出す。次にラベルに応じて次に取り出す箱を決める。球はその都度元の箱へ戻す

  • ラベルが0なら箱A
  • ラベルが1なら箱B

マルコフ連鎖の定義

Xn;n=0,1,2,X_n;n = 0, 1, 2, \cdots が任意の n0n \geq 0 と、任意の SS の要素 x0,x1,,i,jx_0, x_1, \cdots, i, j について次のような性質を持つとき、マルコフ連鎖という

P(Xn+1=j  X0=x0,X1=x1,,Xn=i)=P(Xn+1=j  Xn=i)P(X_{n+1} = j \ |\ X_0 = x_0, X_1 = x_1, \cdots, X_n = i) = P(X_{n+1} = j \ |\ X_n = i)

2状態マルコフ連鎖

ある学習モデルを考える。一連の実験において正答した次は確率 aa で誤答し、誤答した次は確率 bb で正答する。このマルコフ連鎖の推移図と推移行列を求めよ。(ただし、正答の状態を0、誤答の状態を1とおく)

マルコフ連鎖の時間ダイナミクス

  • pijp_{ij} は状態 ii から状態 jj へ一段階で推移する確率
  • 状態 ii から状態 jj へちょうど nn 段かかって推移する確率を pij(n)p^{(n)}_{ij} とおく (nn 段階の推移確率)
  • 時間に一様な場合なので mm に無関係

pij(n)=P(Xm+n=j  Xm=i)p_{ij}^{(n)} = P(X_{m+n} = j \ |\ X_m = i)

查普曼-柯尔莫格洛夫方程(Chapman-Kolmogorov equation,C-K equation)

チャップマン・コルモゴロフの等式(Chapman–Kolmogorov equation):一般に (n+m)(n + m) 段階では

pij(n+m)=pi0(n)p0j(m)+pi1(n)p1j(m)++pik(n)pkj(m)+=kpik(n)pkj(m)p_{ij}^{(n + m)} = p_{i0}^{(n)} p_{0j}^{(m)} + p_{i1}^{(n)} p_{1j}^{(m)} + \cdots + p_{ik}^{(n)} p_{kj}^{(m)} + \cdots = \sum_k p_{ik}^{(n)} p_{kj}^{(m)}

nn 段階の推移確率 pij(n)p_{ij}^{(n)}(i,j)(i, j) 成分とする行列を P(n)P^{(n)} (nn 段階の推移行列)と書く

P(n+m)=P(n)P(m)P^{(n+m)} = P^{(n)} P^{(m)}

  • P(1)=PP^{(1)} = P
  • P(2)=P(1+1)=PP=P2P^{(2)} = P^{(1+1)} = PP = P^2
  • P(3)=P(2+1)=P2P=P3P^{(3)} = P^{(2+1)} = P^2P = P^3
  • P(n)=PnP^{(n)} = P^n

初期確率 αi\alpha_i (状態 ii から始まる確率)が次のように与えられているとする

αi=P(X0=i),i=0,1,2,(iαi=1)\alpha_i = P(X_0 = i), \quad i = 0, 1, 2, \cdots \quad (\sum_i \alpha_i = 1)

マルコフ連鎖が nn 段階後に状態 jj にある確率 P(Xn=j)P(X_n = j) は、

P(Xn=j)=iαipij(n)P(X_n = j) = \sum_i \alpha_i p_{ij}^{(n)}

したがって、初期確率ベクトル α=(α0,,αi,)\mathbf{\alpha} = (\alpha_0, \cdots, \alpha_i, \cdots) とおくと P(Xn=j)P(X_n = j)αPn\mathbf{\alpha}P^n(i,j)(i, j) 成分である

例題2.1

α=(0,0,1)\alpha = (0, 0, 1)

と置く

αP2=(0.50.40.40.30.40.30.20.30.5)(0.50.40.40.30.40.30.20.30.5)=(0.29,0.35,0.36)\begin{aligned} \mathbf{\alpha}P^2 = \begin{pmatrix} 0.5 & 0.4 & 0.4 \\ 0.3 & 0.4 & 0.3 \\ 0.2 & 0.3 & 0.5 \end{pmatrix} \begin{pmatrix} 0.5 & 0.4 & 0.4 \\ 0.3 & 0.4 & 0.3 \\ 0.2 & 0.3 & 0.5 \end{pmatrix} = (0.29, 0.35, 0.36) \end{aligned}

状態の分類

ある整数 n0n \geq 0 に対し、pi,j(n)>0p_{i,j}^{(n)} > 0 の時、状態 jj は状態 ii から到達可能であるといい、iji \rightarrow j と表す。i,ji, j が互いに到達可能であるとき、iji \leftrightarrow j と表す

関係 \leftrightarrow は次の同値関係を満たす

  • iji \leftrightarrow j:反射法則
  • iji \leftrightarrow j ならば jij \leftrightarrow i:対称法則
  • iji \leftrightarrow j および jkj \leftrightarrow k ならば、iki \leftrightarrow k:推移法則

状態の全体は、同値な組に分割することができる(同値な組に含まれる状態は互いに到達可能)

同値な組がただ一つなとき、そのマルコフ連鎖は既約であるという

  • 少なくとも1個のループができている = 既約

例題2.2

pii=1p_{ii} = 1 の時、状態 ii吸収的である

第一个状态转移图有错(14,34\frac{1}{4}, \frac{3}{4} 颠倒),图片待修改!!!
第二个状态转移图也有错,p00=1p_{00} = 1 ,图片待修改!!!

状態 ii から出発するマルコフ連鎖が状態 ii初めて戻るまでに要する時間(段階数)を τi\tau_i とする。すなわち

τi=min{n1;Xn=i}\tau_i = \min \lbrace n \geq 1;X_n = i \rbrace

この時、τi\tau_i が有限となる確率が1の時、状態 ii再帰的であるという。一方、1よりも小さい時、状態 ii一時的であるという

  • 状態 ii再帰的、状態 ii と状態 jj は互いに到達可能 \rightarrow 状態 jj も再帰的
  • 状態 ii一時的、状態 ii と状態 jj は互いに到達可能 \rightarrow 状態 jj も一時的

再帰的、一時的の性質:

  • 規約なマルコフ連鎖においては全ての状態は再帰的か一時的かのどちらか
    • どれか一つの状態をチェックすればよい
  • 有限マルコフ連鎖(有限個の状態からなる連鎖)においては、状態 ii が一時的であるための必要十分条件は、iji \rightarrow j であるが iji \nLeftarrow j となる jj があることである
    • 有限マルコフ連鎖の一時的な組は"閉じてない"
    • 少なくとも一つの再帰的な状態が存在する(規約なら再帰的)
  • 状態 ii が再帰的であれば、ii から出発すると何回も ii に戻ってくる。したがって、状態 ii が再帰的であるための必要十分条件は、ii から出発して ii へ戻る回数の平均が無限大になることである

式で書くと

ln={1,Xn=10,Xn1T=n=1ln\begin{aligned} ln = \begin{cases} 1, \quad X_n = 1 \\ 0, \quad X_n \neq 1 \end{cases} \qquad T = \sum_{n=1}^{\infty} l_n \end{aligned}

とおき、マルコフ連鎖が状態 ii から出発するときの条件的確率で測った TT の期待値を求めると、

Ei[T]=n=1Ei[ln]=n=11×P(Xn=i  X0=i)=n=1pii(n)E_i[T] = \sum_{n=1}^{\infty} E_i[l_n] = \sum_{n=1}^{\infty} 1 \times P(X_n = i \ |\ X_0 = i) = \sum_{n=1}^{\infty}p_{ii}^{(n)}

となる。したがって、状態 ii

  • n=1pii(n)=\sum_{n=1}^{\infty}p_{ii}^{(n)} = \infty ならば再帰的
  • n=1pii(n)<\sum_{n=1}^{\infty}p_{ii}^{(n)} < \infty ならば一時的

別の表現として

状態 ii から出発して、nn 段階目ではじめて jj に到達する確率を fij(n)f^{(n)}_{ij}初到達確率)とおく

fij(1)=pijf_{ij}^{(1)} = p_{ij}

fij=n=1fij(n):最終到達確率f_{ij} = \sum_{n=1}^{\infty}f_{ij}^{(n)}: \text{最終到達確率}

  • n=1fii(n)=\sum_{n=1}^{\infty}f_{ii}^{(n)} = \infty ならば再帰的
  • n=1fii(n)<\sum_{n=1}^{\infty}f_{ii}^{(n)} < \infty ならば一時的

例題2.3

マルコフ連鎖の再帰性と一時性

思路:求解 n=1pii(n)\sum_{n=1}^{\infty}p_{ii}^{(n)} 即可

  • p00(n)=1n=1p00(n)=p_{00}^{(n)} = 1 \rightarrow \sum_{n=1}^{\infty} p_{00}^{(n)} = \infty:状態 0 は再帰的
  • n=1p11(n)=n=1(12)n=12112=1\sum_{n=1}^{\infty} p_{11}^{(n)} = \sum_{n=1}^{\infty} (\frac{1}{2})^n = \frac{\frac{1}{2}}{1 - 1 - 2} = 1:状態 1 は一時的

进入0的话一直循环,如果是1的话,只会进入1次0,一旦进入0的话就出不来了

状態 ii再帰的なとき、τi\tau_i平均値

μi=n=1nfii(n)\mu_i = \sum_{n=1}^{\infty} nf_{ii}^{(n)}

平均再帰時間といい、μi<,μi=\mu_i < \infty, \mu_i = \infty に応じて ii をそれぞれ正状態、零状態という

状態 ii に対して、pii(n)>0p_{ii}^{(n)} > 0 を満たす整数 n1n \geq 1全体の最大公約数を状態 ii周期といい、did_i で表す。di=1d_i = 1 の時、ii非周期的であるという。また、全ての nn に対して pii(n)=0p_{ii}^{(n)} = 0 の時は di=1d_i = 1 とする。非周期的、再帰的、正状態をエルゴード状態という

例題2.4

マルコフ連鎖の周期

对于状态 ii,满足 pii(n)>0p_{ii}^{(n)} > 0 的整数 n1n \geq 1 的所有最大公约数为状态 ii周期di=1d_i = 1 的时候,ii 是非周期。

第二个周期为2的例子,是以01循环,故周期为2

例題2.5

今日の天気が晴れなので

α=(1,0,0)\alpha = (1, 0, 0)

  • αP=(0.4,0.6,0)\alpha \mathbf{P} = (0.4, 0.6, 0)
  • αP2=(0.28,0.54,0.18)\alpha \mathbf{P}^2 = (0.28, 0.54, 0.18)

マルコフ連鎖の定常分布

マルコフ連鎖が既約でエルゴード状態である場合,推移行列の極限値が存在する

πj=limnpij(n)\pi_j = \lim_{n \rightarrow \infty}p_{ij}^{(n)}

これは状態 ii無関係であり、

πj=1μj,j=0,1,2,\pi_j = \frac{1}{\mu_j}, \quad j = 0, 1, 2, \cdots

μj\mu_j は状態 jj平均再帰時間、この πj\pi_j は次の方程式を解いて得られ,解は正の値としてただ一通りに決まる

πP=π,π=(π0,π1,),j=0πj=1\color{blue}{\mathbf{\pi} P = \mathbf{\pi}, \qquad \mathbf{\pi} = (\pi_0, \pi_1, \cdots), \qquad \sum_{j=0}^{\infty} \pi_j = 1}

時間が経過しても確率分布が変化しない場合,その分布を定常分布という

πP=πi=0πipij=πj,j=0,1,2,\mathbf{\pi} P = \mathbf{\pi} \quad \leftrightarrow \quad \sum_{i=0}^{\infty} \pi_i p_{ij} = \pi_j, \quad j = 0, 1, 2, \cdots

この方程式の解 π\piマルコフ連鎖の定常分布という

確率分布 π\pi が存在するための必要十分条件はマルコフ連鎖が再帰的で正状態であること

既約なマルコフ連鎖のすべての状態は再帰的な正状態、再帰的な零状態、一時的のいずれかである

既約な有限マルコフ連鎖は再帰的な正状態

例題2.6

πP=π\pi \mathbf{P} = \pi

を解くと

π=(1985,4885,1885)\pi = (\frac{19}{85}, \frac{48}{85}, \frac{18}{85})

π=(a,b,c)\pi = (a, b, c) 计算即可,注意 a+b+c=1a + b + c = 1

吸収マルコフ連鎖

吸収マルコフ連鎖:状態空間が一時的もしくは吸収的な状態からなるマルコフ連鎖

kk 個の吸収状態nkn - k 個の一時状態がある吸収マルコフ連鎖

P=(I0RQ)\begin{aligned} \color{blue}{\mathbf{P} = \begin{pmatrix} \mathbf{I} & \mathbf{0} \\ \mathbf{R} & \mathbf{Q} \end{pmatrix}} \end{aligned}

  • I\mathbf{I}k×kk \times k単位行列
  • R\mathbf{R}一時状態から吸収状態への推移に関する (nk)×k(n - k) \times k 行列
  • Q\mathbf{Q}一時状態から一時状態への推移に関する (nk)×(nk)(n - k) \times (n - k) 行列
  • 0\mathbf{0}:ゼロ行列

Pn=(I0RnQn)\begin{aligned} \mathbf{P}^n = \begin{pmatrix} \mathbf{I} & \mathbf{0} \\ \mathbf{R}_n & \mathbf{Q}^n \end{pmatrix} \end{aligned}

ただし、

Rn=(I+Q++Qn1)R=(IQn)(IQ)1R\begin{aligned} \mathbf{R}_n &= (\mathbf{I} + \mathbf{Q} + \cdots + \mathbf{Q}^{n-1}) \mathbf{R} \\ &= (\mathbf{I} - \mathbf{Q}^n) (\mathbf{I} - \mathbf{Q})^{-1} \mathbf{R} \end{aligned}

Q\mathbf{Q}非負行列で,行の和は1以下であるため、limnQn=0\lim_{n \rightarrow \infty} \mathbf{Q}^n = 0 より

limnRn=(IQ)1R=MR\lim_{n \rightarrow \infty} \mathbf{R}_n = (\mathbf{I} - \mathbf{Q})^{-1} \mathbf{R} = \mathbf{MR}

P=(I0MR0)\begin{aligned} \color{blue}{\mathbf{P}^{\infty} = \begin{pmatrix} \mathbf{I} & \mathbf{0} \\ \mathbf{MR} & \mathbf{0} \end{pmatrix}} \end{aligned}

  • M\mathbf{M}基本行列
  • MR\mathbf{MR}吸収確率

状態 i,ji, j が一時状態の時、ii から出発して吸収状態に達するまでの間に jj を訪問する平均回数mijm_{ij} とする

mijm_{ij} は,一時状態 ii から出発して吸収状態に達するまでに平均何時間(何段階)を一時状態 jj で過ごすかを表している

mij=δij+tpitmtjm_{ij} = \delta_{ij} + \sum_t p_{it} m_{tj}

状態 ii から一段階で吸収状態に達するか、途中の一時状態 tt を経て最終的に吸収状態に達するかのどちらか

M=(mij)\mathbf{M} = (m_{ij}) と置いて上式を行列で表せば、平均滞在時間(平均訪問回数)の行列 M\mathbf{M} は次式で表される

M=I+QMM=(IQ)1\mathbf{M} = \mathbf{I} + \mathbf{QM} \quad \leftrightarrow \quad \color{blue}{\mathbf{M} = (\mathbf{I} - \mathbf{Q})^{-1}}

例題2.7

吸収マルコフ連鎖

先按照吸収マルコフ連鎖定义将推移矩阵分为4部分,并通过 M=(IQ)1\mathbf{M} = (\mathbf{I} - \mathbf{Q})^{-1} 求解 M\mathbf{M}MR\mathbf{MR}

与单位矩阵同一列的矩阵是表示从一时状态到吸收状态

从一时状态到吸收状态所需平均时间看矩阵 M\mathbf{M}。例如状态2到吸收状态的平均时间记为矩阵 M\mathbf{M} 中表示状态2的那一行加起来即可

到吸收状态的概率看吸收概率矩阵 MR\mathbf{MR}

1次元ランダム・ウォーク

数直線上を確率 pp で右へ1だけ移動するか,確率 q=1pq = 1 - p で左へ1だけ移動する粒子の運動

状態 ii から状態 ii に戻る推移確率を考える

  • 奇数ステップ (2n1)(2n - 1) で戻る確率:pii(2n1)=0,n=1,2,p_{ii}^{(2n - 1)} = 0, \quad n = 1, 2, \cdots
  • 偶数ステップ (2n)(2n) で戻る確率:pii(2n)=(2n)!n!n!(p(1p))np_{ii}^{(2n)} = \frac{(2n)!}{n!n!} (p(1 - p))^n

スターリングの公式

n!2πn(ne)nn! \simeq \sqrt{2\pi n}(\frac{n}{e})^n

より、p=12p = \frac{1}{2} の時

n=1pii(2n)=n=1(4p(1p))nπn=n=11πn>1πn=11n=\sum_{n=1}^{\infty}p_{ii}^{(2n)} = \sum_{n=1}^{\infty} \frac{(4p (1 - p))^n}{\sqrt{\pi n}} = \sum_{n=1}^{\infty} \frac{1}{\sqrt{\pi n}} > \frac{1}{\sqrt{\pi}} \sum_{n=1}^{\infty}\frac{1}{n} = \infty

よって,再帰的

p12p \neq \frac{1}{2} の時、r=4p(1p),(r<1)r = 4p(-1p), (|r| < 1) と置くと

n=1pii(2n)=n=1rnπn<1πn=1rn=1πr1r<\sum_{n=1}^{\infty}p_{ii}^{(2n)} = \sum_{n=1}^{\infty} \frac{r^n}{\sqrt{\pi n}} < \frac{1}{\sqrt{\pi}} \sum_{n=1}^{\infty} r^n = \frac{1}{\sqrt{\pi}} \frac{r}{1 - r} < \infty

よって,一時的

左右等確率で移動する場合を、1次元の対称なランダム・ウォークという

1次元と2次元の対称な (p=1/2)(p = 1/2) ランダム・ウォークは再帰的、3次元以上の対称なランダム・ウォークは一時的である

  • 原点回りをうろつくよりも、一方の側に偏っていることが多い。原因:原点右側にいる割合を xx,左側にいる割合を 1x1 - x とすると、その確率密度関数 f(x)f(x)f(x)=1πx(1x)f(x) = \frac{1}{\pi \sqrt{x(1 - x)}} となることが知られている(逆正弦法則)

例題2.8

P=(3107102313)\begin{aligned} \mathbf{P} = \begin{pmatrix} \frac{3}{10} & \frac{7}{10} \\\\ \frac{2}{3} & \frac{1}{3} \end{pmatrix} \end{aligned}

定常分布:πP=π\pi \mathbf{P} = \pi

例題2.9

ポアソン過程

泊松过程是一种典型的独立增量过程,它具有连续参数 tt 与分离散状态取值,因而也是连续参数马尔科夫链。在通信,交通,日常零售业务等各个领域的研究中,泊松过程是各类问题建模时最常用的一种输入模型

e-mailの着信数,交通事故の発生数,新築の住宅数,外国人の流入数,機械の修理数,サービス窓口の行列人数などの事象を考える

t=0t = 0 で観測を始め,時間 tt までに起こる発生数を N(t)N(t) とすると、N(t)N(t) はどのような性質を持つ?

時間10までに起こる交通事故の発生数 N(10)N(10) は、時間10と15の間に起こる交通事故の発生数 N(15)N(10)N(15) - N(10) とは無関係

独立増分:互いに排反な時間区分に起こる事象の発生数が独立

N(15)N(10)N(15) - N(10) は時間の幅 1510=515 - 10 = 5 のみにより定まる確率変数

定常増分:時間 sst(s<t)t(s < t) の間に起こる事象の発生数の分布が時間の幅 (ts)(t - s) のみに依存して決まる

ポアソン過程の定義

連続時間 t0t \geq 0 で離散値 {i;i=0,1,2,}\lbrace i; i = 0, 1, 2, \cdots \rbrace を取るマルコフ過程のうち、次のような性質を持つものをポアソン過程という

  1. N(0)=0N(0) = 0
  2. N(t)N(t)定常独立増分を持つ
  3. 十分小さい時間の幅 hh に対して,P(N(h)=1)=λh+o(h)P(N(h) = 1) = \lambda h + o(h).ただし λ>0\lambda > 0
  4. 十分小さい時間の幅 hh に対して,P(N(h)2)=o(h)P(N(h) \geq 2) = o(h)

ここで,limh0=f(h)h=0\lim_{h \rightarrow 0} = \frac{f(h)}{h} = 0 が成り立つとき,f(h)=o(h)f(h) = o(h) と表す hh が十分小さいとき,o(h)o(h) は無視できる

P(N(h)=1)λhP(N(h)2)0\begin{aligned} &P(N(h) = 1) \simeq \lambda h \\ &P(N(h) \geq 2) \simeq 0 \end{aligned}

ポアソン過程の性質

  1. 観測開始時点 t=0t = 0 では事象の発生はない
  2. t0<t1<<tnt_0 < t_1 < \cdots < t_n の時,確率変数 N(ti)N(ti1)N(t_i) - N(t_{i - 1})N(tj)N(tj1)N(t_j) - N(t_{j - 1})独立 (ij)(i \neq j) (独立増分)
  3. 確率変数 N(t2+h)N(t1+h)N(t_2 + h) - N(t_1 + h)N(t2)N(t1)N(t_2) - N(t_1) は,任意の h>0h > 0 に対し同じ分布を持つ(定常増分)
  4. 十分小さい時間幅 hh の間に発生する事象の数が1回である確率は hh に比例
  5. hh の間に発生する事象の数が2回以上発生する確率は0
  6. hh の間に事象が発生しない確率は 1λh1 - \lambda h

y=N(t)y = N(t) は,高さ1で飛躍する階段状の関数となる.τi\tau_i は確率変数である

ポアソン分布

時間 tt 内に事象が kk 個発生する確率 Pk(t)P_k(t)

N(t)N(t) は平均 λt\lambda tポアソン分布をしている

Pk(t)=eλt(λt)kk!,k=0,1,2,\color{blue}{P_k(t) = e^{-\lambda t} \frac{(\lambda t)^k}{k!}, k = 0, 1, 2, \cdots}

  • E[N(t)]=V[N(t)]=λtE[N(t)] = V[N(t)] = \lambda t

例題3.1

事象の発生確率–ポアソン分布

アクセス数は λ=0.5\lambda = 0.5 のポアソン過程をなすと考える

Pk(t)=eλt(λt)kk!,k=0,1,2,\color{blue}{P_k(t) = e^{-\lambda t} \frac{(\lambda t)^k}{k!}, k = 0, 1, 2, \cdots}

(1) 丁度3

P3(2)=e113!=0.0613P_3(2) = e^{-1} \frac{1}{3!} = 0.0613

(2) 3以下

P0(2)+P1(2)+P2(2)+P3(2)=e1(1+11!+12!+13!)=0.981P_0(2) + P_1(2) + P_2(2) + P_3(2) = e^{-1} (1 + \frac{1}{1!} + \frac{1}{2!} + \frac{1}{3!}) = 0.981

(3) 3以上

1(P0(2)+P1(2)+P2(2))=10.9197=0.08031 - (P_0(2) + P_1(2) + P_2(2)) = 1 - 0.9197 = 0.0803

事象の発生時間間隔の分布

発生時間間隔 TiT_i は独立,かつ,パラメータ λ\lambda指数分布に従う

E[Ti]=1λ,V[Ti]=1λ2E[T_i] = \frac{1}{\lambda}, \quad V[T_i] = \frac{1}{\lambda^2}

P(Ti>t)=eλtP(T_i > t) = e^{-\lambda t}

指数分布の無記憶性

P(X>t+s  X>s)=P(X>t)P(X > t + s \ |\ X > s) = P(X > t)

を満たす確率分布は指数分布のみ

事象の発生確率がポアソン分布、事象の発生時間間隔は指数分布

例題3.2

事象の発生時間間隔–指数分布

P(Ti>t)=eλt\color{blue}{P(T_i > t) = e^{-\lambda t}}

(1) 2分以上

P(T2)=1P(T<2)=1(1e1)=0.3679P(T \geq 2) = 1 - P(T < 2) = 1 - (1 - e^{-1}) = 0.3679

(2) 4分以内

P(T4)=1P(T>4)=1e2=0.8647P(T \leq 4) = 1 - P(T > 4) = 1 - e^{-2} = 0.8647

n番目の事象が起こるまで時間 SnS_n が従う分布

SnS_n はパラメータ n,λn, \lambdaガンマ分布 gn,λg_{n, \lambda} に従う

gn,λ(t)=λeλt(λt)n1(n1)!,t0g_{n, \lambda}(t) = \lambda e^{-\lambda t} \frac{(\lambda t)^{n-1}}{(n - 1)!}, \quad t \geq 0

E[Sn]=nλ,V[Sn]=nλ2E[S_n] = \frac{n}{\lambda}, \quad V[S_n] = \frac{n}{\lambda^2}

例題3.3

(1) E[Sn]=nλE[S_n] = \frac{n}{\lambda} より

E[S100]=1001=100E[S_{100}] = \frac{100}{1} = 100

(2) P(Ti>t)=eλtP(T_i > t) = e^{-\lambda t} より

P(T101>2)=eλt=e2=0.1353P(T_{101} > 2) = e^{-\lambda t} = e^{-2} = 0.1353

(3) Pk(t)=eλt(λt)kk!,k=0,1,2,P_k(t) = e^{-\lambda t} \frac{(\lambda t)^k}{k!}, k = 0, 1, 2, \cdots より

P(N(6)=0)=P0(6)=e6=0.0025P(N(6) = 0) = P_0(6) = e^{-6} = 0.0025

例題3.4

[1] ポアソン分布

(1)

Pn(t)=eλt(λt)nn!P_n(t) = e^{-\lambda t} \cdot \frac{(\lambda t)^n}{n!}

(2)

P0(T)=eλtP_0(T) = e^{-\lambda t}

[2] 一時間あたり受け取るメール数は λ=4824=2\lambda = \frac{48}{24} = 2

(1) 事象の発生確率がポアソン分布に従う。

P3(6)=e26(26)33!=0.00177P_3(6) = e^{-2 * 6} \cdot \frac{(2 * 6)^3}{3!} = 0.00177

(2) 事象の発生時間間隔は指数分布に従う。P(Ti>t)=eλtP(T_i > t) = e^{-\lambda t} より

P(T2)=e22=0.01832P(T \geq 2) = e^{-2 * 2} = 0.01832

(3) 事象の発生確率がポアソン分布に従う。E[N(t)]=V[N(t)]=λtE[N(t)] = V[N(t)] = \lambda t より

E[N(t)]=λt=E[N(3)]=23=6E[N(t)] = \lambda t = E[N(3)] = 2 * 3 = 6

純出生過程

ポアソン過程の最も簡単な拡張を考える

  • 独立性の仮定の除外
  • 時間 tt 内に起こる事象の発生数が jj (状態 jj )のとき、引き続く微小区間で新た1つ発生する確率が jj に依存する場合を考える

状態遷移が隣り合う状態だけで起きる

  • 純出生過程
  • 純死滅過程

に分けて考える

純出生過程

時刻 tt で状態 jj の時

  1. 十分小さい時間幅 hhjj+1j \rightarrow j + 1 に変わる確率 λjh+o(h)\lambda_j h + o(h)
  2. 十分小さい時間幅 hhjj+1j \rightarrow j + 1 以外に変わる確率 o(h)o(h)
  3. 変化しない確率 1λjh+o(h)1 - \lambda_j h + o(h)

ただし,各 λj\lambda_j は正数(出生率に相当)

Pj(t)=λjPj(t)+λλj1Pj1(t),j=1,2,P0(t)=λ0P0(t)\begin{aligned} &P'_j(t) = -\lambda_j P_j(t) + \lambda_{\lambda_{j-1}}P_{j-1}(t), \quad j = 1, 2, \cdots \\ &P'_0(t) = -\lambda_0 P_0(t) \end{aligned}

初期条件 P0(0)=1,Pj(0)=0,j=1,2,P_0(0) = 1, P_j(0) = 0, j = 1, 2, \cdots を用いて,

P0(t)=eλ0tP_0(t) = e^{-\lambda_0 t}

となり,P1(t),P2(t)P_1(t), P_2(t) を順に求めることができる

ユール・ファーリ過程

λj=jλ,(λ>0)\lambda_j = j \lambda, (\lambda > 0) とすると,初期条件

Pj(0)=1,j=1Pj(0)=0,j1\begin{aligned} &P_j(0) = 1, \quad j = 1 \\ &P_j(0) = 0, \quad j \neq 1 \end{aligned}

のもとで,Pj(t)P_j(t) は次式で表される.

Pj(t)=eλt(1eλt)j1,j1P_j(t) = e^{-\lambda t}(1 - e^{\lambda t})^{j - 1}, \quad j \geq 1

P0(t)=0P_0(t) = 0. 線形出生過程ともいわれる

j=1Pj(t)=j=1eλt(1eλt)j1=eλt{11(1eλt)}=1\sum_{j=1}^{\infty} P_j(t) = \sum_{j=1}^{\infty} e^{-\lambda t}(1 - e^{\lambda t})^{j-1} = e^{-\lambda t} \lbrace \frac{1}{1 - (1 - e^{-\lambda t})} \rbrace = 1

純死滅過程

時刻 tt で状態 jj の時

  1. 十分小さい時間幅 hhjj1j \rightarrow j - 1 に変わる確率 μjh+o(h)\mu_j h + o(h)
  2. 十分小さい時間幅 hhjj1j \rightarrow j - 1 以外に変わる確率 o(h)o(h)
  3. 変化しない確率 1μjh+o(h)1 - \mu_j h + o(h)

ただし,各 μj\mu_j は正数(死亡率に相当)

Pj(t)=μjPj(t)+μμj+1Pj+1(t),j=1,2,Pn(t)=μnPn(t)P0(t)=μ1P1(t)\begin{aligned} &P'_j(t) = -\mu_j P_j(t) + \mu_{\mu_{j+1}}P_{j+1}(t), \quad j = 1, 2, \cdots \\ &P'_n(t) = -\mu_n P_n(t) \\ &P'_0(t) = \mu_1 P_1(t) \end{aligned}

初期条件 Pn(0)=1,Pj(0)=0,j=1,2,P_n(0) = 1, P_j(0) = 0, j = 1, 2, \cdots を用いて,

Pn(t)=eμntP_n(t) = e^{-\mu_n t}

となり,Pn1(t),Pn2(t),,P0(t)P_{n-1}(t), P_{n-2}(t), \cdots, P_0(t) を順に求めることができる

例題3.5

  1. 時刻0は Pn(0)P_n(0) 以外は全部0
  2. 通式:Pj(t)=eμt(μt)nj(nj)!,j=1,2,,nP_j(t) = e^{-\mu t} \cdot \frac{(\mu t)^{n-j}}{(n - j)!}, \quad j = 1, 2, \cdots, n

出生死滅過程

  • 純出生過程では状態数の増加のみが発生
  • 純死滅過程では状態数の減少のみが発生

実際のシステムでは,出生と死滅がランダムに発生する場合が多い

時刻 tt で状態 jj の時

  1. 十分小さい時間幅 hhjj+1j \rightarrow j + 1 に変わる確率 λjh+o(h)\lambda_j h + o(h)
  2. 十分小さい時間幅 hhjj1j \rightarrow j - 1 に変わる確率 μjh+o(h)\mu_j h + o(h)
  3. 隣以外に変わる確率 o(h)o(h)
  4. 変化しない確率 1(λj+μj)h+o(h)1 - (\lambda_j + \mu_j) h + o(h)

時間 tt で状態 jj にある確率を Pj(t)P_j(t) とおくと,次式を得る

ddtPj(t)=(λj+μj)Pj(t)+λj1Pj1(t)+μj+1Pj+1(t)\frac{d}{dt}P_j(t) = -(\lambda_j + \mu_j)P_j(t) + \lambda_{j-1}P_{j-1}(t) + \mu_{j+1} P_{j+1}(t)

t=0t = 0 で状態 ii にあるとする.この時,初期条件は

Pj(0)=δij,δij={1,j=i0,ji\begin{aligned} P_j(0) = \delta_{ij}, \quad \delta_{ij} = \begin{cases} 1, \quad j = i \\ 0, \quad j \neq i \end{cases} \end{aligned}

に従う Pj(t)P_j(t)Pi,j(t)P_{i, j}(t) とする

ddtPi,0(t)=λ0Pi,0(t)+μ1Pi,1(t)ddtPi,j(t)=(λi+μj)Pi,j(t)+λj1Pi,j1(t)+μj+1Pi,j+1(t)Pi,j(0)=1,j=iPi,j(0)=0,ji\begin{aligned} &\frac{d}{dt}P_{i, 0}(t) = -\lambda_0 P_{i, 0}(t) + \mu_1 P_{i, 1}(t) \\\\ &\frac{d}{dt}P_{i, j}(t) = -(\lambda_i + \mu_j)P_{i, j}(t) + \lambda_{j - 1} P_{i, j-1}(t) + \mu_{j+1}P_{i, j+1}(t) \\\\ &P_{i, j}(0) = 1, \quad j = i \\\\ &P_{i, j}(0) = 0, \quad j \neq i \end{aligned}

パラメータ {λj>0,μj+1>0.j=0,1,2,}\lbrace \lambda_j > 0, \mu_{j+1} > 0. j = 0, 1, 2, \cdots \rbrace の出生死滅過程において,極限の推移確率

pj=limtPi,j(t),i,j=0,1,2,p_j = \lim_{t \rightarrow \infty} P_{i, j}(t), \quad i, j = 0, 1, 2, \cdots

が初期状態 ii と独立に存在するための必要十分条件

j=0aj<,aj=λj1λj=2λ0μjμj1μ1,a0=1\sum_{j=0}^{\infty} a_j < \infty, a_j = \frac{\lambda_{j-1} \lambda_{j=2} \cdots \lambda_0}{\mu_j \mu_{j-1} \cdots \mu_1}, \quad a_0 = 1

この時,極限の推移確率は次式となる

p0=(j=0)1,pj=ajp0\color{blue}{p_0 = (\sum_{j=0}^{\infty})^{-1}, p_j = a_j p_0}

例題3.6

(1) nn 番目の事象が起こるまで時間 SnS_n が従う分布はガンマ分布なので

E[Sn]=nλ=E[S180]=1803=60[min]E[S_n] = \frac{n}{\lambda} = E[S_{180}] = \frac{180}{3} = 60[min]

(2) P(Snt)=P(N(t)n)=j=1eλt(λt)jj!\color{blue}{P(S_n \leq t) = P(N(t) \geq n) = \sum_{j=1}^{\infty} e^{-\lambda t} \frac{(\lambda t)^j}{j!}}

P(S5>2)=1P(S52)=e6(1+6+622!+633!+644!)=0.2851P(S_5 > 2) = 1 - P(S_5 \leq 2) = e^{-6}(1 + 6 + \frac{6^2}{2!} + \frac{6^3}{3!} + \frac{6^4}{4!}) = 0.2851

(3) 事象の発生時間間隔–指数分布

P(T40.5)=1P(T40.5)=1e30.5=0.7769P(T_4 \leq 0.5) = 1 - P(T_4 \geq 0.5) = 1 - e^{-3 * 0.5} = 0.7769

例題3.7

待ち行列

待ち行列理論の発展

  • 1838: ポアソン分布 ランダムに発生するイベントの統計分布
  • 1900前半: アーランの計算 ポアソン分布を電話の呼び出しが発生する確率に応用し,電話がつながらない確率などを計算(待ち行列理論の始まり)
  • 1953: ケンドールの記号 ケンドールの記号を導入し,待ち行列理論を分類・整理
  • 1957: ジャクソンによるネットワーク問題の解析 待ち行列理論を組み合わせたネットワークの問題を解析(空港の手続きのようにチェックインで並び,セキュリティ検査で並び,というような引き続き起こる待ち行列モデルの解析)
  • 1961: リトルの公式 リトルが待ち行列,待ち時間に関する一般的な公式を導出

待ち行列モデル

待ち行列理論(Queueing Theory):不確定な需要と供給という両面から待ちができる行列モデルを数学モデルとして扱う

  • スーパーのレジ台
    • 到着する客が支払いサービスを受ける
    • サービスが終わると出ていく
    • 順番を待っている次の客がサービスを受ける
  • 情報処理センターのスパコン
    • 到着するジョブが処理のサービスを受ける
    • 終わって空きができると次のジョブが処理を受ける

客、ジョブの到着、支払い/処理のサービス時間は偶然に変化する.どのようにシステムを設計する?

  • 客:サービスを受けようとする個々の要求
  • 窓口:サービスを提供するもの
  • 待ち行列:サービス中の客を除いた実際の待ち行列
  • システム:サービス中の人も含めた全体の待ち行列

行列が一つでサービス窓口が一つの場合を単一待ち行列という(行列が一つで窓口が複数の場合は複数待ち行列)

ケンドールの記号

  • AA:客はどれくらいの間隔で到着する?
  • BB:サービスの時間はどれくらいばらつく?
  • CC:窓口の数?
  • KK:行列の上限は?

ケンドールの記号: A/B/C/KA/B/C/K もしくは A/B/C(K)A/B/C(K)

この表記の場合,

  • 母集団(システムに来る客の数)は (=m)\infty (= m)
  • 待ち行列のルールは先着順(First Come, First Served) (=Z)(= Z)

全て書き出す場合は,

A/B/C/K/m/ZA/B/C/K/m/Z

と表される

客の到着(A)/サービス時間(B)

ポアソン到着(記号:MM

  • 到着がパラメータ λ\lambdaポアソン過程に従う
  • 到着時間間隔は指数分布に従う

F(t)=1eλt,t>0F(t) = 1 - e^{-\lambda t}, \quad t > 0

一定分布(記号:DD

  • 到着間隔が t=λt = \lambda で一定

F(t)={0,t1λ1,t1λ\begin{aligned} F(t) = \begin{cases} 0, \quad t \leq \frac{1}{\lambda} \\ 1, \quad t \geq \frac{1}{\lambda} \end{cases} \end{aligned}

一般分布(記号:GG)

  • 到着間隔は独立で同一の分布 F(t)F(t) に従う

単一待ち行列モデル

M/M/1(INF)

  • 客の到着:パラメータ λ\lambdaポアソン分布
  • サービス時間:パラメータ μ\mu指数分布
  • 窓口の数:1つ
  • 行列の上限:制限なし

客の到着を個体の出生,サービスを受けて退出することを死亡とすると,M/MM/M 型待ち行列モデルは

λi=λ,μi+1=μ(i=0,1,2,)\lambda_i = \lambda, \quad \mu_{i+1} = \mu \quad (i = 0, 1, 2, \cdots)

出生死滅過程で定式化でき,極限の推移確率 pjp_j を求めることができる

極限の推移確率が存在するための必要十分条件は

j=0aj<,aj=λj1λj=2λ0μjμj1μ1=(λμ)j,a0=1\sum_{j=0}^{\infty} a_j < \infty, a_j = \frac{\lambda_{j-1} \lambda_{j=2} \cdots \lambda_0}{\mu_j \mu_{j-1} \cdots \mu_1} = (\frac{\lambda}{\mu})^j, \quad a_0 = 1

ρ=λμ\rho = \frac{\lambda}{\mu} とおくと,j=0ρj<\sum_{j=0}^{\infty} \rho^j < \infty であるから,ρ<1\rho < 1 でなければならない.この時,

j=0ρj=11ρ<\sum_{j=0}^{\infty}\rho^j = \frac{1}{1 - \rho} < \infty

これは,到着率(=出生率) λ\lambdaサービス率(=死亡率) μ\mu より小さいときである

ρ<1\rho < 1 の時,極限の推移確率は,

pj=ajp0=ρjp0,p0=(j=0aj)1=(j=0ρj)1=1ρp_j = a_j p_0 = \rho^j p_0, \quad p_0 = (\sum_{j=0}^{\infty} a_j)^{-1} = (\sum_{j=0}^{\infty} \rho^j)^{-1} = 1 - \rho

これはシステム内に jj 人の客(1人のサービス中の客 +(j1)+ (j - 1) 人の待ち行列の客)がいる確率を表している.ここで,

ρ=λμ=到着率サービス率=1/μ1/λ=平均サービス時間平均到着時間間隔\color{blue}{\rho = \frac{\lambda}{\mu} = \frac{\text{到着率}}{\text{サービス率}} = \frac{1/\mu}{1/\lambda} = \frac{\text{平均サービス時間}}{\text{平均到着時間間隔}}}

トラフィック密度もしくは利用率と呼ばれる

システム内のいる客の平均値および分散は,

L=E(X)=j=0j×pj=ρ1ρ\color{blue}{L = E(X) = \sum_{j=0}^{\infty} j \times p_j = \frac{\rho}{1 - \rho}}

Var(X)=E(X2)E(X)2=ρ(1ρ)2\color{blue}{\text{Var}(X) = E(X^2) - E(X)^2 = \frac{\rho}{(1 - \rho)^2}}

待ち行列の平均の長さは,サービス中の客を除いた平均の待ち行列の長さであるから,

Lq=j=1(j1)×pj=L(1p0)=Lρ=ρ21ρ\color{blue}{L_q = \sum_{j=1}^{\infty} (j - 1) \times p_j = L - (1 - p_0) = L - \rho = \frac{\rho^2}{1 - \rho}}

リトルの公式

待ち人数 = 待ち時間 ×\times 到着率

システム内の客の平均待ち時間を WW待ち行列での客の平均待ち時間を WqW_q とすると,

L=λWLq=λWq\begin{aligned} &L = \lambda W \\ &L_q = \lambda W_q \end{aligned}

が成り立つ.したがって

W=1λL=1μ(1ρ)(平均系内時間)Wq=1λLq=ρμ(1ρ)(平均待ち時間)\begin{aligned} &\color{blue}{W = \frac{1}{\lambda}L = \frac{1}{\mu(1 - \rho)} \quad \text{(平均系内時間)}} \\ &\color{blue}{W_q = \frac{1}{\lambda}L_q = \frac{\rho}{\mu(1 - \rho)} \quad \text{(平均待ち時間)}} \end{aligned}

待ち行列で用いられる特性値

注:7的待ち行列での客の平均待ち時間为 WqW_q

M/M/1(N)

パラメータ λi=λ,μi+1=μ\lambda_i = \lambda, \mu_{i+1} = \mu,有限状態 i=0,1,,Ni = 0, 1, \cdots, N からなる出生死滅過程を考える.

j=1,2,,N1j = 1, 2, \cdots, N - 1 においては,

ddtPi,0(t)=λPi,0(t)+μPi,1(t)ddtPi,j(t)=(λ+μ)Pi,j(t)+λPi,j1(t)+μPi,j+1(t)Pi,j(0)=1,j=iPi,j(0)=0,ji\begin{aligned} &\frac{d}{dt}P_{i, 0}(t) = -\lambda P_{i, 0}(t) + \mu P_{i, 1}(t) \\\\ &\frac{d}{dt}P_{i, j}(t) = -(\lambda + \mu)P_{i, j}(t) + \lambda P_{i, j-1}(t) + \mu P_{i, j+1}(t) \\\\ &P_{i, j}(0) = 1, \quad j = i \\\\ &P_{i, j}(0) = 0, \quad j \neq i \end{aligned}

が成り立つ.j=Nj = N における境界条件を考えると

ddtPi,N(t)=λPi,N1(t)μPi,N(t)\color{blue}{\frac{d}{dt}P_{i, N}(t) = \lambda P_{i, N-1}(t) - \mu P_{i, N}(t)}

が成り立つ.

出生死滅過程の極限の推移行列が存在するための条件は,

j=0aj<,aj=λj1λj=2λ0μjμj1μ1=(λμ)j,a0=1\sum_{j=0}^{\infty} a_j < \infty, a_j = \frac{\lambda_{j-1} \lambda_{j=2} \cdots \lambda_0}{\mu_j \mu_{j-1} \cdots \mu_1} = (\frac{\lambda}{\mu})^j, \quad a_0 = 1

であるが,有限状態の出生死滅過程では常に満足することが明らか j=0Naj<\sum_{j=0}^N a_j < \infty である.従って,

pj=ρjp0\color{blue}{p_j = \rho^j p_0}

ただし,

p0+p1++pN=j=0Npj=1p_0 + p_1 + \cdots + p_N = \sum_{j=0}^N p_j = 1

需要注意的是,与 M/M/1()M/M/1(\infty) 不同的是,M/M/1(N)M/M/1(N) 并没有 ρ<1\rho < 1 的条件

ρ=1\rho = 1 の時,

p0(1+N)=1p0=11+Np_0 (1 + N) = 1 \quad \leftrightarrow \quad p_0 = \frac{1}{1 + N}

ρ1\rho \neq 1 の時,

p0(1+ρ++ρN)=ρ(1ρN+11ρ)=1p0=1ρ1ρN+1p_0(1 + \rho + \cdots + \rho^N) = \rho(\frac{1 - \rho^{N+1}}{1 - \rho}) = 1 \quad \leftrightarrow \quad p_0 = \frac{1 - \rho}{1 - \rho^{N+1}}

よって

pj={11+N,ρ=11ρ1ρN+1ρj,ρ1\begin{aligned} p_j = \begin{cases} \frac{1}{1 + N}, \quad &\rho = 1 \\\\ \frac{1 - \rho}{1 - \rho^{N+1}} \rho^j, \quad &\rho \neq 1 \end{cases} \end{aligned}

pN=1ρ1ρN+1ρN呼損率(棄損率)\color{blue}{p_N = \frac{1 - \rho}{1 - \rho^{N+1}} \rho^N} \quad \text{呼損率(棄損率)}

システム内にいる客の平均値ρ1\rho \neq 1 に対して,

L=j=0Nj×pj=ρ1ρ1(N+1)ρN+NρN+11ρN+1\color{blue}{L = \sum_{j=0}^N j \times p_j = \frac{\rho}{1 - \rho} \frac{1 - (N + 1)\rho^N + N\rho^{N+1}}{1 - \rho^{N+1}}}

ρ\rho が十分小さく ρN0,ρN+10\rho^N \simeq 0, \rho^{N+1} \simeq 0 とできる場合は,

L=ρ1ρ1(N+1)0+N010L = \frac{\rho}{1 - \rho} \frac{1 - (N + 1)0 + N0}{1 - 0}

となり近似的に M/M/1()M/M/1(\infty) モデルの LL と同じになる.

待ち行列の平均長さは,

Lq=j=1N(j1)×pj=ρ21ρ1NρN1+(N1)ρN1ρN+1\color{blue}{L_q = \sum_{j=1}^N (j - 1) \times p_j = \frac{\rho^2}{1 - \rho} \frac{1 - N\rho^{N-1} + (N - 1) \rho^N}{1 - \rho^{N+1}}}

ρ=1\rho = 1 に対しては,

L=j=0Nj×pj=11+NN(N+1)2=N2\color{blue}{L = \sum_{j=0}^N j \times p_j = \frac{1}{1 + N} \frac{N(N+1)}{2} = \frac{N}{2}}

Lq=j=1N(j1)×pj=L(1p0)=N(N1)2(N+1)\color{blue}{L_q = \sum_{j=1}^N (j - 1) \times p_j = L - (1 - p_0) = \frac{N(N - 1)}{2(N + 1)}}

M/M/1(N)M/M/1(N) 待ち行列モデルでもリトルの公式は成り立つ.ただし,システムに NN 人の客がいるとき,新たに待ち行列に加わることができないため,実質到着率 λa\lambda_a

λa=j=0N1pjλ=(1PN)λ=1ρN1ρN+1λ,ρ1\lambda_a = \sum_{j=0}^{N-1} p_j \lambda = (1 - P_N)\lambda = \frac{1 - \rho^N}{1 - \rho^{N+1} \lambda}, \quad \rho \neq 1

システム内の客の平均待ち時間WW待ち行列での客の平均待ち時間WqW_q とすると,

W=1λaL=1μ(1ρ)1(N+1)ρN+NρN+11ρN(平均系内時間)Wq=1λaLq=ρμ(1ρ)1NρN1+(N1)ρN1ρN(平均待ち時間)\begin{aligned} &\color{blue}{W = \frac{1}{\lambda_a}L = \frac{1}{\mu(1 - \rho)} \frac{1 - (N + 1)\rho^N + N\rho^{N+1}}{1 - \rho^N} \quad \text{(平均系内時間)}} \\\\ &\color{blue}{W_q = \frac{1}{\lambda_a}L_q = \frac{\rho}{\mu(1 - \rho)} \frac{1 - N\rho^{N-1} + (N - 1)\rho^N}{1 - \rho^N} \quad \text{(平均待ち時間)}} \end{aligned}

例題4.1

複数待ち行列モデル

M/M/s(INF)

  • 客の到着:パラメータ λ\lambdaポアソン分布
  • サービス時間:パラメータ μ\mu指数分布
  • 窓口の数:s(s=2,3,)s(s = 2, 3, \cdots)
  • 行列の上限:制限なし

サービス率 μk\mu_k は,窓口の数によって異なる.システム内に kk 人の客がいた場合,

μk={kμ,kssμ,k>s\begin{aligned} \mu_k = \begin{cases} &k\mu, \quad k \leq s \\ &s\mu, \quad k > s \end{cases} \end{aligned}

となり,客の数が窓口数よりも多いか少ないかでサービス率が異なる

極限の推移確率が存在するための必要十分条件は

j=0aj=j=0k=1jλk1μk=j=0s1ujj!+uss!n=0ρn<,u=λμ,ρ=λsμ\sum_{j=0}^{\infty} a_j = \sum_{j=0}^{\infty} \prod_{k=1}^j \frac{\lambda_{k-1}}{\mu_k} = \sum_{j=0}^{s - 1} \frac{u^j}{j!} + \frac{u^s}{s!} \sum_{n=0}^{\infty} \rho^n < \infty, \quad \color{blue}{u = \frac{\lambda}{\mu}, \rho = \frac{\lambda}{s\mu}}

であるから,ρ<1\rho < 1 でなければならない.この時,極限の推移確率は,

pj={ujj!p0,jsuss!(us)jsp0,j>s\begin{aligned} p_j = \begin{cases} &\frac{u^j}{j!}p_0, \quad j \leq s \\\\ &\frac{u^s}{s!} (\frac{u}{s})^{j-s}p_0, \quad j > s \end{cases} \end{aligned}

p0=[j=0s1ujj!+uss!(1ρ)]1\color{blue}{p_0 = [\sum_{j=0}^{s-1} \frac{u^j}{j!} + \frac{u^s}{s!(1 - \rho)}]^{-1}}

待ち行列の平均の長さ LqL_q は,窓口の個数が ss であるから,

Lq=j=s(js)pj=p0uss!(0+1ρ+2ρ2+)=p0usρs!(1ρ)2\color{blue}{L_q = \sum_{j=s}^{\infty} (j - s)p_j = p_0 \frac{u^s}{s!}(0 + 1\rho + 2\rho^2 + \cdots) = \frac{p_0 u^s \rho}{s!(1 - \rho)^2}}

システム内にいる客の平均値 LL

L=j=0jpj=j=0s1ujj!p0+j=suss!(us)jsp0=u+Lq\color{blue}{L = \sum_{j=0}^{\infty} jp_j = \sum_{j=0}^{s-1} \frac{u^j}{j!} p_0 + \sum_{j=s}^{\infty} \frac{u^s}{s!} (\frac{u}{s})^{j-s} p_0 = u + L_q}

実質到着率 k=0λkpk=λ\sum_{k=0}^{\infty} \lambda_k p_k = \lambda より,リトルの公式から

W=1λL(平均系内時間)Wq=1λLq(平均待ち時間)\begin{aligned} &\color{blue}{W = \frac{1}{\lambda}L \quad \text{(平均系内時間)}} \\ &\color{blue}{W_q = \frac{1}{\lambda}L_q \quad \text{(平均待ち時間)}} \end{aligned}

M/M/s(N)

  • 客の到着:パラメータ λ\lambdaポアソン分布
  • サービス時間:パラメータ μ\mu指数分布
  • 窓口の数:s(s=2,3,)s(s = 2, 3, \cdots)
  • 行列の上限NsN - s

1j<s1 \leq j < s のとき

pj=ujj!p0,u=λμ\color{blue}{p_j = \frac{u^j}{j!} p_0, u = \frac{\lambda}{\mu}}

sjNs \leq j \leq N のとき

pj=uss!(us)jsp0\color{blue}{p_j = \frac{u^s}{s!} (\frac{u}{s})^{j-s} p_0}

p0=[j=0s1ujj!+uss!(1ρNs+1(1ρ))]1\color{blue}{p_0 = [\sum_{j=0}^{s-1} \frac{u^j}{j!} + \frac{u^s}{s!}(\frac{1 - \rho^{N-s+1}}{(1 - \rho)})]^{-1}}

待ち行列の平均の長さ LqL_q は,窓口の個数が ss であるから,

Lq=j=s(js)pj=uss!ρ(1ρ)2p0{1[(Ns)(1ρ)+1]ρNs}\color{blue}{L_q = \sum_{j=s}^{\infty} (j - s)p_j = \frac{u^s}{s!} \frac{\rho}{(1 - \rho)^2} p_0 \lbrace 1 - [(N - s) (1 - \rho) + 1] \rho^{N - s} \rbrace}

システム内にいる客の平均値 LL は,

L=j=0jpj=Lq+λPNμ\color{blue}{L = \sum_{j=0}^{\infty} jp_j = L_q + \frac{\lambda P_N}{\mu}}

リトルの公式から

W=1λL(平均系内時間)Wq=1λLq(平均待ち時間)\begin{aligned} &\color{blue}{W = \frac{1}{\lambda}L \quad \text{(平均系内時間)}} \\ &\color{blue}{W_q = \frac{1}{\lambda}L_q \quad \text{(平均待ち時間)}} \end{aligned}

例題4.2

例題4.3

平均サービス率を固定して,平均到着率を変化させる. ρ\rho が大きくなるにつれて待ち行列は長くなる.特に0.9を超えると急激に増加する

例題4.4

乱数

エージェントモデルシミュレーション

参考