7 随机向量和随机过程的随机数

7.1 条件分布法

产生随机向量的一种方法是条件分布法。 设\(\boldsymbol X = (X_1, X_2, \dots, X_r)\)的分布密度或分布概率 \(p(x_1, x_2, \dots, x_r)\)可以分解为 \[\begin{aligned} p(x_1, x_2, \dots, x_r) = p(x_1) p(x_2 | x_1) p(x_3 | x_1, x_2) \cdots p(x_r | x_1, x_2, \dots, x_{r-1}) \end{aligned}\] 则可以先生成\(X_1\), 由已知的\(X_1\)的值从条件分布\(p(x_2 | x_1)\)产生\(X_2\), 再从已知的\(X_1, X_2\)的值从条件分布\(p(x_3 | x_1, x_2)\)产生\(X_3\), 如此重复直到产生\(X_r\)

7.1.1 多项分布随机数

进行了\(n\)次独立重复试验\(Y_1, Y_2, \dots, Y_n\), 每次的试验结果\(Y_i\)\(1,2,\dots,r\)中取值, \(P(Y_i = j) = p_j, j=1,2,\dots,r; i=1,2,\dots,n\)。 令\(X_j\)为这\(n\)次试验中结果\(j\)的个数(\(j=1,2,\dots,r\)), 称\(\boldsymbol X = (X_1, X_2, \dots, X_r)\)服从多项分布, 其联合概率函数为 \[\begin{aligned} P(X_j=x_j, j=1,2,\dots,r) = \frac{n!}{x_1! x_2! \dots x_r!} p_1^{x_1} p_2^{x_2} \dots p_r^{x_r}. \end{aligned}\] (其中\(\sum_{j=1}^r x_j = n\))。

要生成\(\boldsymbol X\)的随机数,考虑不同结果数\(r\)与试验次数\(n\)的比较。 当\(r\)\(n\)相比很大时, 每个结果的出现次数都不多,而且许多结果可能根本不出现, 这样,可以模拟产生\(Y_i, i=1,2,\dots,n\)然后从\(\{ Y_i \}\)中计数得到 \((X_1, X_2, \dots, X_r)\)。 当\(r\)\(n\)相比较小时, 每个结果出现次数都比较多,可以使用条件分布逐个地产生\(X_1, X_2, \dots, X_r\)

\(X_1\)表示\(n\)次重复试验中结果1的出现次数, 以结果1作为成功,其它结果作为失败, 显然\(X_1 \sim \text{B}(n, p_1)\)。 产生\(X_1\)的值\(x_1\)后, \(n\)次试验剩余的\(n - x_1\)次试验结果只有\(2,3,\dots,r\)可取, 于是在\(X_1 = x_1\)条件下结果2的出现概率是 \[\begin{aligned} P(Y_i = 2 | Y_i \neq 1) = \frac{P(Y_i=2)}{\sum_{k=2}^r P(Y_i=k)} = \frac{p_2}{1 - p_1} \end{aligned}\] 于是剩下的\(n - x_1\)次试验中结果2的出现次数服从\(\text{B}(n-x_1, \frac{p_2}{1-p_1})\)分布, 从这个条件分布产生\(X_2 = x_2\)。 类似地, 剩下的\(n - x_1 - x_2\)次试验中结果3的出现次数服从 \(\textbf{B}(n-x_1-x_2, \frac{p_3}{1-p_1-p_2})\)分布, 从这个条件分布中抽取\(X_3=x_3\), 如此重复可以产生多项分布\(\boldsymbol X\)的随机数\((x_1, x_2, \dots, x_r)\)

测试:

7.2 多元正态分布模拟

设随机向量\(\boldsymbol X = (X_1, X_2, \dots, X_p)^T\)服从多元正态分布 \(\text{N}(\boldsymbol \mu, \Sigma)\), 联合密度函数为 \[\begin{aligned} p(\boldsymbol x) = (2\pi)^{-\frac{p}{2}} |\Sigma|^{-\frac12} \exp\left\{ -\frac12 (\boldsymbol x - \boldsymbol \mu)^T \Sigma^{-1} (\boldsymbol x - \boldsymbol \mu) \right\}, \ \boldsymbol x \in R^p \end{aligned}\] 正定矩阵\(\Sigma\)有Cholesky分解\(\Sigma = C C^T\), 其中\(C\)为下三角矩阵。 设\(\boldsymbol Z = (Z_1, Z_2, \dots, Z_p)^T\)服从\(p\)元标准正态分布 \(\text{N}(\boldsymbol 0, I)\)(\(I\)表示单位阵), 则\(\boldsymbol X = \boldsymbol\mu + C \boldsymbol Z\) 服从\(\text{N}(\boldsymbol \mu, \Sigma)\)分布。

随机数发生器的程序:

测试:

7.3 泊松过程模拟

参数为\(\lambda\)的泊松过程\(N(t), t\geq 0\)是取整数值的随机过程, \(N(t)\)表示到时刻\(t\)为止到来的事件个数, 两次事件到来的时间间隔服从指数分布Exp(\(\lambda\))。

所以,如果要生成泊松过程前\(n\)个事件到来的时间, 只要生成\(n\)个独立的Exp(\(\lambda\))随机数\(X_1, X_2, \dots, X_n\), 则\(S_k = \sum_{i=1}^k X_i, k=1,2,\dots,n\)为各个事件到来的时间。

如果要生成泊松过程在时刻\(T\)之前的状态, 只要知道发生在\(T\)之前的所有事件到来时间就可以了。

R程序:

测试:

生成泊松过程在时刻\(T\)之前的状态的另外一种方法是 先生成\(N(T) \sim \text{Poisson}(\lambda)\), 设\(N(T)\)的值为\(n\), 再生成\(n\)个独立的U(0,1)随机变量\(U_1, U_2, \dots, U_n\), 从小到大排序为\(U_{(1)} \leq U_{(2)} \leq \dots U_{(n)}\), 则\((T U_{(1)} , T U_{(2)} , \dots, T U_{(n)})\)为时刻\(T\)之前的所有事件到来时间。

为了生成强度函数为\(\lambda(t), t\geq 0\)的非齐次泊松过程到时刻\(T\)为止的状态, 如果\(\lambda(t)\)满足 \[\begin{aligned} \lambda(t) \leq M, \ \forall t \geq 0\end{aligned}\] 则可以按照生成参数为\(M\)的齐次泊松过程的方法去生成各个事件到来时刻, 但是以\(\lambda(t)/M\)的概率实际记录各个时刻。 R程序如下:

测试。假设前50分钟速率为0.2, 后50分钟速率为1。

7.4 平稳时间序列模拟

平稳时间序列中的ARMA模型可以递推生成。

对AR(\(p\))模型 \[\begin{aligned} X_t = a_1 X_{t-1} + a_2 X_{t-2} + \dots + a_p X_{t-p} + \varepsilon_t, \ t \in \mathbb{Z} \end{aligned}\] (\(\mathbb{Z}\)表示所有整数的集合), 其中\(\{ \varepsilon_t \}\)为白噪声WN(0,\(\sigma^2\)), 即方差都为\(\sigma^2\)、彼此不相关的随机变量序列, \(\{ \varepsilon_t \}\)可以用N(0,\(\sigma^2\))分布的独立序列来模拟, 也可以用其它有二阶矩的分布。 因为AR(\(p\))模型具有所谓“稳定性”, 所以我们从任意初值出发按照递推\(N_0\)步(\(N_0\)是一个较大正整数)后, 再继续递推\(N\)步后得到的\(N\)\(X_t\)就可以作为上述AR(\(p\))模型的一次实现。 根据模型稳定性的好坏, \(N_0\)可取为\(50 \sim 1000\)之间。

下面的R程序直接使用递推,比较简单。

测试,用如下AR(1)模型: \[ X_t = 0.5*X_{t-1} + \varepsilon_t, \ \varepsilon_t \sim \text{N}(0, 2^2) \]

因为R的循环很慢, 所以改用filter()函数写成如下版本:

测试:

对于MA(\(q\))模型 \[\begin{aligned} X_t = \varepsilon_t + b_1 \varepsilon_{t-1} + \dots + b_q \varepsilon_{t-q}, \ t \in \mathbb{Z} \end{aligned}\] 生成\(\varepsilon_{1-q}, \varepsilon_{2-q}, \dots, \varepsilon_0, \varepsilon_1, \dots, \varepsilon_N\)后可以直接对\(t=1,2,\dots,N\)按公式 计算得到\(\{X_t, t=1,2,\dots,N \}\)

模拟MA(\(q\))的R程序:

测试,模型为

\[ X_t = \varepsilon_t + 0.5 \varepsilon_{t-1}, \ \varepsilon_t \sim \text{N}(0, 2^2) \]

对于ARMA(\(p,q\))模型 \[\begin{aligned} X_t = a_1 X_{t-1} + a_2 X_{t-2} + \dots + a_p X_{t-p} + \varepsilon_t + b_1 \varepsilon_{t-1} + \dots + b_q \varepsilon_{t-q}, \ t \in \mathbb{Z} \end{aligned}\] 也可以从初值零出发递推生成\(N_0 + N\)个然后只取最后的\(N\)个作为 ARMA(\(p,q\))模型的一次实现。

模拟程序如下:

R函数filter可以进行递推和卷积计算。 R函数arima.sim可以进行ARMA模型和ARIMA模型模拟。

习题

习题1

设坛子中有\(n\)个不同颜色的球, 第\(i\)中颜色的球有\(n_i\)个(\(i=1,2,\dots,r\))。 从坛子中随机无放回地抽取\(m\)个球, 设随机变量\(X_i\)表示取出的第\(i\)种颜色球的个数, 设计高效的算法模拟\((X_1, X_2, \dots, X_r)\)的值。

习题2

利用稀松法编写模拟到时刻\(T=10\)为止的强度为 \[\begin{aligned} \lambda(t) = 3 + \frac{4}{t+1} \end{aligned}\] 的非齐次泊松过程的R程序; 设法改进这个算法的效率。

习题3

用R的filter()函数编写AR(\(p\))模型的模拟程序。

习题4

用R的filter()函数编写MA(\(q\))模型的模拟程序。