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是不同结果个数
<- function(n, m, prob){
rng.multinom <- length(prob)
r <- c(0, cumsum(prob))
pcum <- matrix(0, n, r)
xmat for(i in seq(n)){
<- m
m1 for(k in 1:(r-1)){
<- rbinom(1, m1, prob[k]/(1 - pcum[k]))
xmat[i, k] <- m1 - xmat[i, k]
m1
}<- m1
xmat[i,r]
}
xmat }
测试:
<- c(0.1, 0.3, 0.6)
prob <- rng.multinom(100000, 5, prob)
x # (1,2,2)的概率:
<- 30*prod(prob^c(1,2,2)); p122
p122 ## [1] 0.0972
mean(x[,1]==1 & x[,2]==2 & x[,3]==2)
## [1] 0.0992
# (1,1,3)的概率:
<- 20*prod(prob^c(1,1,3)); p113
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是协方差阵,需要是对称正定阵。
<- function(n, mu, Sigma){
rng.mnorm <- length(mu)
m <- chol(Sigma) # Sigma = M' M
M <- matrix(rnorm(n*m), n, m) %*% M
y for(j in seq(along=mu)){
<- y[,j] + mu[j]
y[,j]
}
y }
测试:
<- rng.mnorm(1000, c(3,2), rbind(c(4, 1), c(1, 1)))
x 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程序:
<- function(T=100, lambda=1){
rng.poisproc <- max(c(100, round(2*T*lambda)))
nmax <- numeric(nmax) # 用来保存事件到来时刻
xvec <- 0 # 当前时刻
S <- 0 # 已发生事件个数
nev repeat{
<- -log(runif(1))/lambda
X <- S + X
S if(S <= T){
<- nev + 1
nev <- S
xvec[nev] if(nev == nmax){ # 需要扩大存储空间
<- c(xvec, numeric(nmax))
xvec <- nmax*2
nmax
}else {
} break
}
}1:nev]
xvec[ }
测试:
<- 100; lambda <- 0.5
T <- rng.poisproc(T, lambda)
x 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程序如下:
<- function(T=100, lambda=function(t) 1, M=1){
rng.poisproc.nh <- max(c(100, round(2*T*M)))
nmax <- numeric(nmax) # 用来保存事件到来时刻
xvec <- 0 # 当前时刻
S <- 0 # 已发生事件个数
nev repeat{
<- -log(runif(1))/M
X <- S + X
S if(S <= T){
if(runif(1) <= lambda(S)/M){
<- nev + 1
nev <- S
xvec[nev] # 否则不记录此事件
} if(nev == nmax){ # 需要扩大存储空间
<- c(xvec, numeric(nmax))
xvec <- nmax*2
nmax
}else {
} break
}
}1:nev]
xvec[ }
测试。假设前50分钟速率为0.2, 后50分钟速率为1。
<- 100; lambda <- function(t) ifelse(t<=50, 0.4, 1)
T <- rng.poisproc.nh(T, lambda, M=1)
x 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: 需要舍弃的初始部分
<- function(n, a, sigma=1.0, n0=1000){
ar.gen.v1 <- n0 + n
n2 <- length(a)
p <- rev(a)
arev <- numeric(n2)
x2 <- rnorm(n2, 0, sigma)
eps for(tt in seq(p+1, n2)){
<- eps[tt] + sum(x2[(tt-p):(tt-1)] * arev)
x2[tt]
}<- x2[(n0+1):n2]
x <- ts(x)
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) \]
<- ar.gen.v1(100, c(0.5), sigma=2)
x plot(x)
pacf(x)
因为R的循环很慢,
所以改用stats::filter()
函数写成如下版本:
<- function(n, a, sigma=1.0,
ar.gen n0=1000,
x0=numeric(length(a))){
<- n0 + n
n2 <- length(a)
p <- rnorm(n2, 0, sigma)
eps <- stats::filter(
x2 method="recursive", sides=1, init=x0)
eps, a, <- x2[(n0+1):n2]
x <- ts(x)
x attr(x, "model") <- "AR"
attr(x, "a") <- a
attr(x, "sigma") <- sigma
x }
测试:
<- ar.gen(100, c(0.5), sigma=2)
x 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程序:
<- function(n, a, sigma=1.0){
ma.gen <- length(a)
q <- n+q
n2 <- rnorm(n2, 0, sigma)
eps <- stats::filter(
x2 c(1,a), method="convolution", sides=1)
eps, <- x2[(q+1):n2]
x <- ts(x)
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) \]
<- ma.gen(1000, c(0.5), sigma=2)
x 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\))模型的一次实现。
模拟程序如下:
<- function(
arma.gen sigma=1.0,
n, a, b, n0=1000, x0=numeric(length(a))){
<- n0 + n
n2 <- length(a)
p ## first generate n0+n MA(q) series
<- ma.gen(n2, b, sigma=sigma)
eps <- stats::filter(
x2 method="recursive", sides=1, init=x0)
eps, a, <- x2[(n0+1):n2]
x <- ts(x)
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\))模型的模拟程序。