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)\)。

## 输入
##   - n: 需要输出的随机数个数
##   - m: 模型中的独立重复试验次数
##   - prob: 每次试验，各个结果出现的概率
## 输出：n行r列矩阵，r是不同结果个数
rng.multinom <- function(n, m, prob){
  r <- length(prob)
  pcum <- c(0, cumsum(prob))
  xmat <- matrix(0, n, r)
  for(i in seq(n)){
    m1 <- m
    for(k in 1:(r-1)){
      xmat[i, k] <- rbinom(1, m1, prob[k]/(1 - pcum[k]))
      m1 <- m1 - xmat[i, k]
    }
    xmat[i,r] <- m1
  }
  
  xmat
}

测试：

prob <- c(0.1, 0.3, 0.6)
x <- rng.multinom(100000, 5, prob)
# (1,2,2)的概率：
p122 <- 30*prod(prob^c(1,2,2)); p122
## [1] 0.0972
mean(x[,1]==1 & x[,2]==2 & x[,3]==2)
## [1] 0.0992
# (1,1,3)的概率：
p113 <- 20*prod(prob^c(1,1,3)); p113
## [1] 0.1296
mean(x[,1]==1 & x[,2]==1 & x[,3]==3)
## [1] 0.12805

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)\)分布。

随机数发生器的程序：

## 生成多元正态分布随机数，
## n是需要的个数，
## mu是均值向量，
## Sigma是协方差阵，需要是对称正定阵。
rng.mnorm <- function(n, mu, Sigma){
  m <- length(mu)
  M <- chol(Sigma) # Sigma = M' M
  y <- matrix(rnorm(n*m), n, m) %*% M
  for(j in seq(along=mu)){
    y[,j] <- y[,j] + mu[j]
  }
  y
}

测试：

x <- rng.mnorm(1000, c(3,2), rbind(c(4, 1), c(1, 1)))
plot(x[,1], x[,2], type="p", cex=0.1)
abline(v=3, h=2, col="green")

var(x)
##          [,1]     [,2]
## [1,] 4.151643 1.125089
## [2,] 1.125089 1.070800

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程序：

rng.poisproc <- function(T=100, lambda=1){
  nmax <- max(c(100, round(2*T*lambda)))
  xvec <- numeric(nmax) # 用来保存事件到来时刻
  S <- 0   # 当前时刻
  nev <- 0 # 已发生事件个数
  repeat{
    X <- -log(runif(1))/lambda
    S <- S + X
    if(S <= T){
      nev <- nev + 1
      xvec[nev] <- S
      if(nev == nmax){ # 需要扩大存储空间
        xvec <- c(xvec, numeric(nmax))
        nmax <- nmax*2
      }
    } else {
      break
    }
  }
  xvec[1:nev]
}

测试：

T <- 100; lambda <- 0.5
x <- rng.poisproc(T, lambda)
plot(c(0, T), c(0, 1), type="n", xlab="time", ylab="")
abline(v=x, col="red")

1/mean(diff(x)) # 估计的速率
## [1] 0.4375458

生成泊松过程在时刻\(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程序如下:

rng.poisproc.nh <- function(T=100, lambda=function(t) 1, M=1){
  nmax <- max(c(100, round(2*T*M)))
  xvec <- numeric(nmax) # 用来保存事件到来时刻
  S <- 0   # 当前时刻
  nev <- 0 # 已发生事件个数
  repeat{
    X <- -log(runif(1))/M
    S <- S + X
    if(S <= T){
      if(runif(1) <= lambda(S)/M){
        nev <- nev + 1
        xvec[nev] <- S
      } # 否则不记录此事件
      if(nev == nmax){ # 需要扩大存储空间
        xvec <- c(xvec, numeric(nmax))
        nmax <- nmax*2
      }
    } else {
      break
    }
  }
  xvec[1:nev]
}

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

T <- 100; lambda <- function(t) ifelse(t<=50, 0.4, 1)
x <- rng.poisproc.nh(T, lambda, M=1)
plot(c(0, T), c(0, 1), type="n", xlab="time", ylab="")
abline(v=50, lwd=2)
abline(v=x, col="red")

1/mean(diff(x[x<=50])) # 估计的第一段速率
## [1] 0.3227109
1/mean(diff(x[x>50])) # 估计的第二段速率
## [1] 0.9067809

7.4 布朗运动模拟

时间区间\([0, T]\)上的一维的标准布朗运动为随机过程\(\{ W(t), 0 \leq t \leq T \}\)，满足：

\(W(0)=0\);
以概率1地，\(W(t)\)作为\(t\)的函数，（每条轨道）在\([0,T]\)连续；
独立增量；
\(W(t) - W(s) \sim \text{N}(0, t-s)\), 对任意\(0 \leq s < t \leq T\)。

\(W(t)\)分布为 \[ W(t) \sim \text{N}(0, t), \ 0 < t \leq T. \]

对标准布朗运动\(W(t)\)，令 \[ X(t) = \mu t + \sigma W(t), 0 \leq t \leq T, \] 称\(X(t)\)为有漂移\(\mu\)和扩散系数\(\sigma^2\)的布朗运动，简记为BM(\(\mu, \sigma^2\))。其中\(\mu\)，\(\sigma\)为常数，\(\sigma>0\)。 \(X(t)\)仍是高斯过程，从0出发，有独立增量，轨道连续，边缘分布为 \[ X(t) \sim \text{N}(\mu t, \sigma^2 t) . \]

也可以构造\(X(0)=x\)的解，只要给每个\(X(t)\)加\(x\)即可。

可以定义漂移和扩散系数时变的布朗运动，存在时变的函数\(\mu(t)\)和\(\sigma(t)\)，使得增量\(X(t) - X(s)\)(\(0 \leq s < t \leq T\))服从正态分布，期望为 \[ E[X(t) - X(s)] = \int_s^t \mu(u) \,du, \] 方差为 \[ \text{Var}[X(t) - X(s)] = \int_s^t \sigma^2(u) \,du . \]

设\(0 < t_1 < t_2 < \dots < t_n \leq T\)，要生成\((W(t_1), \dots, W(t_n))\)或者\((X(t_1), \dots, X(t_n))\)的抽样，最简单的方法是利用布朗运动的独立增量性和正态性。设\(Z_1, \dots, Z_n\)是独立同标准正态分布的随机变量（随机数）， \(t_0=0\), \(W(0)=0\), \(X(0)=0\)。

令 \[\begin{align} W(t_{i+1}) = W(t_i) + \sqrt{t_{i+1} - t_i} Z_{i+1}, \ i=0,1,\dots, n-1 . \tag{7.1} \end{align}\]

对\(X(t) \sim \text{BM}(\mu, \sigma^2)\)，令 \[\begin{align} X(t_{i+1}) = X(t_i) + \mu \cdot (t_{i+1} - t_i) + \sigma \sqrt{t_{i+1} - t_i} Z_{i+1}, \ i=0,1,\dots, n-1 . \tag{7.2} \end{align}\]

对于漂移和扩散系数时变的\(X(t)\)，令 \[\begin{align} X(t_{i+1}) = X(t_i) + \int_{t_i}^{t_{i+1}}\mu(u) \,du + \left( \int_{t_i}^{t_{i+1}}\sigma^2(u) \,du \right)^{1/2} Z_{i+1}, \ i=0,1,\dots, n-1 . \tag{7.3} \end{align}\]

只要输入的\(Z_1, \dots, Z_n\)的联合分布正确，则输出的\((W(t_1), \dots, W(t_n))\)和\((X(t_1), \dots, X(t_n))\)的联合分布也是正确的，没有离散化误差。

这只输出了\(t_1, t_2, \dots, t_n\)这有限个点上的轨道的值，中间的值未知，如果用通常的线性插值等方法近似，就会产生离散化误差，分布也不准确。

在漂移和扩散系数时变时，如果将(7.3)中的积分替换成如下的欧拉(Euler)近似： \[ X(t_{i+1}) = X(t_i) + \mu(t_i) \cdot (t_{i+1} - t_i) + \sigma(t_i) \sqrt{t_{i+1} - t_i} Z_{i+1}, \ i=0,1,\dots, n-1 , \] 则输出的\((X(t_1), \dots, X(t_n))\)的联合分布也是近似的，带有离散化误差(discretization error)，而且这种误差是累积的，时间越往后误差越大。

7.5 平稳时间序列模拟

平稳时间序列中的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程序直接使用递推，比较简单。

## 输入: 
##   - n: 需要输出的序列长度
##   - a: 自回归系数$(a_1, a_2, \dots, a_p)，是递推方程右侧的系数
##   - sigma: 新息项的标准差
##   - n0: 需要舍弃的初始部分
ar.gen.v1 <- function(n, a, sigma=1.0, n0=1000){
  n2 <- n0 + n
  p <- length(a)
  arev <- rev(a)
  x2 <- numeric(n2)
  eps <- rnorm(n2, 0, sigma)
  for(tt in seq(p+1, n2)){
    x2[tt] <- eps[tt] + sum(x2[(tt-p):(tt-1)] * arev)
  }
  x <- x2[(n0+1):n2]
  x <- ts(x)
  attr(x, "model") <- "AR"
  attr(x, "a") <- a
  attr(x, "sigma") <- sigma
  
  x
}

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

x <- ar.gen.v1(100, c(0.5), sigma=2)
plot(x)

pacf(x)

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

ar.gen <- function(n, a, sigma=1.0, 
                   n0=1000,
                   x0=numeric(length(a))){
  n2 <- n0 + n
  p <- length(a)
  eps <- rnorm(n2, 0, sigma)
  x2 <- stats::filter(
    eps, a, method="recursive", sides=1, init=x0)
  x <- x2[(n0+1):n2]
  x <- ts(x)
  attr(x, "model") <- "AR"
  attr(x, "a") <- a
  attr(x, "sigma") <- sigma

  x
}

测试：

x <- ar.gen(100, c(0.5), sigma=2)
plot(x)

pacf(x)

对于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程序：

ma.gen <- function(n, a, sigma=1.0){
  q <- length(a)
  n2 <- n+q
  eps <- rnorm(n2, 0, sigma)
  x2 <- stats::filter(
    eps, c(1,a), method="convolution", sides=1)
  x <- x2[(q+1):n2]
  x <- ts(x)
  attr(x, "model") <- "MA"
  attr(x, "b") <- a
  attr(x, "sigma") <- sigma
  x
}

测试，模型为

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

x <- ma.gen(1000, c(0.5), sigma=2)
plot(x)

acf(x)

对于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\))模型的一次实现。

模拟程序如下:

arma.gen <- function(
  n, a, b, sigma=1.0,
  n0=1000, x0=numeric(length(a))){
  n2 <- n0 + n
  p <- length(a)
  ## first generate n0+n MA(q) series
  eps <- ma.gen(n2, b, sigma=sigma)
  x2 <- stats::filter(
    eps, a, method="recursive", sides=1, init=x0)
  x <- x2[(n0+1):n2]
  x <- ts(x)
  attr(x, "model") <- "ARMA"
  attr(x, "a") <- a
  attr(x, "b") <- b
  attr(x, "sigma") <- sigma
  x
}

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\))模型的模拟程序。