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
}测试：
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")
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")
生成泊松过程在时刻\(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")
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程序直接使用递推，比较简单。
## 输入:
## - 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) \]
因为R的循环很慢，
所以改用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 <- filter(eps, a, method="recursive", side=1, init=x0)
x <- x2[(n0+1):n2]
x <- ts(x)
attr(x, "model") <- "AR"
attr(x, "a") <- a
attr(x, "sigma") <- sigma
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 <- filter(eps, c(1,a), method="convolution", side=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) \]
对于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,
by.roots=by.roots.ma,
plot.it=FALSE)
x2 <- filter(eps, a, method="recursive", side=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\))模型的模拟程序。