北京大学Julia语言入门讲义第16章： 随机模拟

平均值法

随机模拟的基本思想可以用随机模拟积分来介绍，
考虑积分问题

∫01h(x)dx，
若取

X∼U(0,1)，
产生

X1,X2,…,XN独立同

U(0,1)分布，
则因为

E[h(X)]=∫01h(x)dx，而

E[h(X)]≈1N∑i=1Nh(Xi),

所以可以用方程的右侧公式来估计积分

∫01h(x)dx。

例如，

∫01exdx=e−1≈1.718282，
用Julia程序模拟估计，取

N=106：

using Random
using Statistics 
Random.seed!(101)
mean(exp.(rand(1_000_000)))
## 1.7178412037826383

为了计算

∫abh(x)dx，
设

X∼U(a,b)，
即

X服从

(a,b)上的均匀分布，
因为

E[h(X)]=∫abh(x)1b−adx=1b−a∫abh(x)dx,

所以可以生成服从

(a,b)上均匀分布的随机变量

Xi，

i=1,2,…,N,
计算

b−aN∑i=1Nh(Xi)作为

∫abh(x)dx的估计值。

# 重复调用随机数函数的错误
function sim_demo(N=1_000_000)
    s = 0.0
    for i=1:N
        if rand() < 0.5
            s += rand()
        else
            s += 1 - rand()
        end
    end

    s / N
end

using Random
Random.seed!(101)
sim_demo()
## 0.5001518755483094

function sim_demo(N=1_000_000)
    s = 0.0
    for i=1:N
        u = rand()
        if u < 0.5
            s += u
        else
            s += 1 - u
        end
    end

    s / N
end
sim_demo()
## 0.24983656295114867

function sim_demo(N=1_000_000)
    sum(x -> x < 0.5 ? x : 1-x, rand(N)) / N
end
sim_demo()
## 0.25019553248724047

重要抽样法

产生均匀分布随机数估计定积分的方法无法计算无穷区间的积分

∫−∞∞f(x)dx；
另外，
在计算有限区间积分

∫01f(x)dx时，
如果

f(x)在

[0,1]中某一段的绝对值很大，
其它位置的绝对值很小，
则当

Xi落入

f(x)绝对值很小的区间时，

Xi对最后估计结果的贡献很小。
所以，
希望不使用均匀分布抽样，
而是抽取随机数

Xi使得

Xi的取值在

f(x)绝对值大的位置出现更密集，
在

f(x)绝对值小的位置出现较稀疏，
这样能够合理利用所有的抽样点

Xi。

设

X∼g(⋅)，

g(x)是区间

(a,b)上的一个分布密度（

(a,b)可以是无穷区间），

i=1,2,…,N,
为了计算

∫abh(x)dx，
注意

∫abh(x)dx=∫abh(x)g(x)g(x)dx=E[h(X)g(X)],

所以可以取

Xi∼g(⋅),

i=1,2,…,N，
估计

∫abh(x)dx为

1N∑i=1Nh(Xi)g(Xi).

这种方法称为重要抽样法(importance sampling)，
其中

g(x)称为试抽样密度(trial density)，
应取

g(x)使其形状与

|h(x)|尽可能接近。

例如，仍考虑

∫01exdx=1.718282的问题。
函数

h(x)=ex是增函数，
所以应取试抽样密度

g(x)在左边较小，
在右边较大。
对

ex用一阶泰勒展开近似有

ex≈1+x，
所以取试抽样密度为

g(x)=c(1+x)=23(1+x),

要产生

g(x)的随机数可以用逆变换法，
密度

g(x)的分布函数

G(x)的反函数为

G−1(y)=1+3y−1, 0<y<1.

因此，取

Uiiid U(0,1),
令

Xi=1+3Ui−1,

i=1,2,…,N,
则重要抽样法的积分公式为

I^=1N∑i=1NeXi23(1+Xi).

Julia程序实现：

g(x) = 2/3*(1+x)
Ginv(p) = sqrt(1 + 3p) - 1
Random.seed!(101)
mean(map((x -> exp(x) / g(x)) ∘ Ginv, rand(1_000_000)))
## 1.7181425608987528

其中

∘表示两个函数的复合运算。
map函数将某个一元函数作用到一个序列的每个元素，
上面的程序如果写成广播格式则会比较繁琐，
难以阅读。

注意，
如果将程序写成下面这样就又犯了使用不同随机数的错误：

mean(exp(Ginv(rand())) / g(Ginv(rand)) for _=1:1_000_000)

这时分子和分母使用了不同的随机数，
结果错误。

条件分布法

产生随机向量的一种方法是条件分布法。
设

X=(X1,X2,…,Xr)的分布密度或分布概率

p(x1,x2,…,xr)可以分解为

p(x1,x2,…,xr)=p(x1)p(x2|x1)p(x3|x1,x2)⋯p(xr|x1,x2,…,xr−1)

则可以先生成

X1， 由已知的

X1的值从条件分布

p(x2|x1)产生

X2，
再从已知的

X1,X2的值从条件分布

p(x3|x1,x2)产生

X3，
如此重复直到产生

Xr。

变换法

设随机向量

X=(X1,X2,…,Xp)T服从多元正态分布

N(μ,Σ)， 联合密度函数为

p(x)=(2π)−p2|Σ|−12exp⁡{−12(x−μ)TΣ−1(x−μ)}, x∈Rp

正定矩阵

Σ有Cholesky分解

Σ=LLT，
其中

L为下三角矩阵。
设

Z=(Z1,Z2,…,Zp)T服从

p元标准正态分布

N(0,I)(

I表示单位阵)，
则

X=μ+LZ
服从

N(μ,Σ)分布。

Julia的LinearAlgebra包提供了chol函数。
chol(A).L返回的矩阵

L使得

LTL=A。

# 泊松过程模拟
function randpproc(T=100.0, λ=1.0)
    local xs, S, nev, X 
    xs = Float64[]
    S = 0.0 # 当前时刻
    nev = 0 # 已发生事件个数
    while true 
        X = -log(rand()) / λ
        S += X
        if S <= T 
            nev += 1
            push!(xs, S)
        else
            break 
        end
    end # while true
    return xs 
end

using Random, Distributions 
using CairoMakie, AlgebraOfGraphics 
CairoMakie.activate!()
Random.seed!(101)
x = randpproc(50, 1)
plt = data((; x=x)) * mapping(:x) * visual(VLines)
draw(plt)
## 估计速率
mean(x[2:end] .- x[begin:end-1])
## 0.8870887289991723

function randnhpproc(T, λ = x -> 1.0, M=1.0)
    local xs 
    xs = randpproc(T, M)
    xs = [t for t in xs if rand() < λ(t) / M]
    return xs 
end 

fl(t) = t <= 50 ? 0.2 : 1.0
Random.seed!(101)
x = randnhpproc(100.0, fl, 1.0)
plt = data((; x=x)) * mapping(:x) * visual(VLines)
dr = draw(plt; axis = (; width=600, height=100))

布朗运动模拟

布朗运动是具有独立增量的高斯过程，
其模拟比较简单。

时间区间

[0,T]上的一维的标准布朗运动为随机过程

{W(t),0≤t≤T}，
满足：

• W(0)=0

  ;

• 以概率1地， W(t)

  作为 t

  的函数，（每条轨道）在 [0,T]

  连续；

• 独立增量；
• W(t)−W(s)∼N(0,t−s)

  ,
  对任意 0≤s<t≤T

  。

W(t)分布为

W(t)∼N(0,t), 0<t≤T.

对标准布朗运动

W(t)，令

X(t)=μt+σW(t),0≤t≤T,

称

X(t)为有漂移

μ和扩散系数

σ2的布朗运动，
简记为BM(

μ,σ2)。
其中

μ，

σ为常数，

σ>0。

X(t)仍是高斯过程，
从0出发，
有独立增量，
轨道连续，
边缘分布为

X(t)∼N(μt,σ2t).

也可以构造

X(0)=x的解，
只要给每个

X(t)加

x即可。

可以定义漂移和扩散系数时变的布朗运动，
存在时变的函数

μ(t)和

σ(t)，
使得增量

X(t)−X(s)(

0≤s<t≤T)服从正态分布，期望为

E[X(t)−X(s)]=∫stμ(u)du,

方差为

Var[X(t)−X(s)]=∫stσ2(u)du.

设

0<t1<t2<⋯<tn≤T，
要生成

(W(t1),…,W(tn))或者

(X(t1),…,X(tn))的抽样，
最简单的方法是利用布朗运动的独立增量性和正态性。
设

Z1,…,Zn是独立同标准正态分布的随机变量（随机数），

t0=0,

W(0)=0,

X(0)=0。

令

(16.1)W(ti+1)=W(ti)+ti+1−tiZi+1, i=0,1,…,n−1.

对

X(t)∼BM(μ,σ2)，令

(16.2)X(ti+1)=X(ti)+μ⋅(ti+1−ti)+σti+1−tiZi+1, i=0,1,…,n−1.

对于漂移和扩散系数时变的

X(t)，令

(16.3)X(ti+1)=X(ti)+∫titi+1μ(u)du+(∫titi+1σ2(u)du)1/2Zi+1, i=0,1,…,n−1.

只要输入的

Z1,…,Zn的联合分布正确，
则输出的

(W(t1),…,W(tn))和

(X(t1),…,X(tn))的联合分布也是正确的，
没有离散化误差。

这只输出了

t1,t2,…,tn这有限个点上的轨道的值，
中间的值未知，
如果用通常的线性插值等方法近似，
就会产生离散化误差，
分布也不准确。

在漂移和扩散系数时变时，
如果将(16.3)中的积分替换成如下的欧拉(Euler)近似：

X(ti+1)=X(ti)+μ(ti)⋅(ti+1−ti)+σ(ti)ti+1−tiZi+1, i=0,1,…,n−1,

则输出的

(X(t1),…,X(tn))的联合分布也是近似的，
带有离散化误差(discretization error)，
而且这种误差是累积的，
时间越往后误差越大。

如下的函数输入时间点的数组，
输出模拟的标准布朗运动轨道在这些时间点的值：

function sim_bm(ts::Vector)
    xs = zeros(length(ts))
    t0 = 0
    x0 = 0.0
    for i in eachindex(ts)
        x1 = x0 + sqrt(ts[i] - t0)*randn()
        xs[i] = x1
        t0 = ts[i]
        x0 = x1 
    end
    return (ts, xs)
end
function sim_bm(T::Real, n::Integer)
    ts = collect(1:n) .* (T/n)
    return sim_bm(ts)
end

测试：

using Random 
using CairoMakie
CairoMakie.activate!()
Random.seed!(101)
n = 100
ts, xs = sim_bm(10.0, 100)
lines(ts, xs)

采样点较少情形

布朗运动是正交增量过程，
所以只要在新的时间点模拟出独立的增量即可。
一般的高斯过程不一定有独立增量。

比如，
设随机过程

X(t),t≥0满足：

X(0)=0,

E[X(t)]=0，

Cov[X(t),X(s)]=exp⁡(−(t−s)22ℓ2),

记右侧的函数为

σ(t,s)，
称其为高斯核函数，
这是一个对称正定核，即：

σ(t,s)=σ(s,t);∀n≥1, t1,…,tn,(σ(ti,tj))n×n≥0.

高斯核函数中的参数

ℓ代表不同时间的状态的相关性随时间距离衰减的平均时间，
各

Xt,Xs之间都是正相关，

ℓ越大，
相关性随时间距离衰减越慢，
不同时刻的状态的相关性越强，
这时轨道也越光滑；

ℓ越小，
轨道的局部变动也越剧烈。

实际上，
因为

σ(s,t)是

s−t的函数，
所以这是一个高斯平稳过程，
即任意的有限维分布都服从多元正态分布，
均值为常数，
协方差函数

σ(s,t)仅依赖于

t−s的随机过程。

为了模拟这个过程在时间点

t1,…,tn上的值，
当

n较小时，
因为

(X(t1),…,X(tn))服从多元正态分布，
均值向量为

0，
协方差阵为

Σ=(σ(ti,tj))n×n,

所以可以对

Σ作Cholesky分解得

Σ=LLT，
其中

L为满秩下三角阵，
令

Z为

n元标准正态分布随机向量，
则

LZ服从多元正态分布

N(0,Σ)。

Julia的Distributions包的MvNormal(mu, Sigma)函数可以用来表示多元正态分布。
为了生成

m条轨道，方法如：

dist = MvNormal(zeros(n), Sigma)
sampmat = rand(dist, m)

结果是一个

n×m矩阵，
每列是随机过程一条轨道。

但是，这样做法有一个问题，
就是虽然在理论上

Σ矩阵是正定矩阵，
但是当采样点个数

n较大时，
这个矩阵很接近于不满秩，
使得Cholesky分解的数值算法出错。
为此，借鉴(Kochenderfer 2019)中的做法，
将

Σ替换成

Σ+ϵI，
其中

ϵ是较小的正常数，如

10−6。

生成高斯核的高斯过程的Julia函数如下：

using Distributions
# 给定表示时间点的向量ts, 需要的轨道条数npath，
# 高斯核参数decay, 保证Cholesky分解可行的扰动量inflate,
# 输出模拟结果为n*m矩阵，每列为一条轨道.
function sim_Gauss(ts::Vector, npath::Integer, 
    decay=1.0; inflate = 1e-6)
    n = length(ts)
    ker(t, s) = exp(-(t-s)^2 / (2*decay^2))
    Sig = [ker(ts[i], ts[j]) for i=1:n, j=1:n]
    dist = MvNormal(zeros(n), Sig + inflate .* I)
    xs = rand(dist, npath)
    return ts, xs
end

# 输入要生成的时间长度T，等分时间点点数n的模拟函数
function sim_Gauss(T::Real, n::Integer, npath=1, 
    decay=1.0; inflate=1e-6)
    ts = collect(1:n) .* (T/n)
    return sim_Gauss(ts, npath, decay; inflate=inflate)
end

测试，模拟100个时间点，取

ℓ=2.0，
模拟5条轨道：

using CairoMakie
CairoMakie.activate!()

ts, xs = sim_Gauss(10.0, 100, 5, 2.0)
fig = lines(ts, xs[:,1])
for i=2:5 ; lines!(ts, xs[:,i]); end
fig 

结果图形：

因为取的相关性衰减时间参数

ℓ=2.0的值较大，
所以每条轨道都很光滑。

下面取

ℓ=0.2:

ts, xs = sim_Gauss(10.0, 100, 5, 0.2)

结果图形：

因为相关性减弱，
所以轨道变得不那么光滑。

平稳高斯列

若

X(t)为平稳高斯列，
均值为0，
协方差函数为

γk,k=0,1,2,…，
则有限维分布完全由

{γk}决定。
为了生成

X(t)的一条长度为

n的轨道，
可以递推地进行抽样。
抽取

X(1)∼N(0,γ0)，
然后，
对

X(t+1),

t=1,…,n−1，
因为

X(t+1)|X(t),…,X(1)是条件正态分布，
有

E(X(t+1)|X(t),…,X(1))=∑j=1tat,jX(t+1−j)=at,1X(t)+at,2X(t−1)+⋯+at,tX(1),Var(X(t+1)|X(t),…,X(1))=σt2,

其中

σ02,σ12,…和

at,j,j=1,…,t,

t=1,2,…可以用如下的Levinson递推公式从

{γk}序列计算得到：

σ02=γ0a1,1=γ1/γ0σk2=σk−12(1−ak,k2)ak+1,k+1=γk+1−ak,1γk−ak,2γk−1−⋯−ak,kγ1γ0−ak,1γ1−ak,2γ2−⋯−ak,kγk(16.4)ak+1,j=ak,j−ak+1,k+1ak,k+1−j,1≤j≤k

这样，就可以递推地计算出

X(t+1)|X(t),…,X(1)的条件期望和条件方差，
从这个条件分布抽样得到

X(t+1)。

在递推过程中，
如果平稳高斯序列

{X(t)}相关性较强，
有可能出现

σk2→0(

k→∞)的情况。
这时，
可以设置一个阈值

ϵ，
当

σp2<ϵ时，
对

t>p就不再用所有的

X(1),…,X(t−1)计算条件分布，
而是仅使用

X(t−1),…,X(t−p)这

p个，
递推预测系数统一使用

ap,j,

j=1,…,p。

下面的函数gauss_levinson函数输入长度为

n的自协方差函数序列

{γk,k=0,1,…,n−1}，
要求的轨道条数npath，
输出一个结果矩阵，
每行为一条轨道的

n个观测值。
也可以输入一个计算

γk的函数，
此函数输入非负整数值

k，
返回

γk的值，
用这个函数规定自协方差函数，
进行模拟。

function gauss_levinson(gam_vec::Vector{Float64}, npath=1; eps=1e-6)
    n = length(gam_vec)
    x = zeros(npath, n)
    a1 = zeros(n)
    a2 = zeros(n)
    S2 = zeros(n)

    S2[1+0] = gam_vec[0+1] # 用0个去预报下一个的方差
    x[:,1] .= rand_norm(zeros(npath), sqrt(S2[1]), npath)

    a1[1] = gam_vec[1+1]/gam_vec[0+1] # a_{1,1}，用X_1预报X_2的系数

    # 下面是Levinson递推计算条件分布的条件期望的线性组合系数，
    # 与条件分布的方差，
    # 并从条件分布抽样
    p = -1
    for k=1:(n-1)
        # 从前k个预报第k+1个的方差
        S2[1+k] = S2[k] * (1 - a1[k]^2)
        #@show S2[1+k]
        # 如果基本可完全线性预报，不需要更多预报系数
        if S2[1+k] < eps 
            p = k 
            break 
        end

        # 从前k个给出条件期望和条件方差，抽取第k+1个
        cond_mean = sum(a1[j] .* x[:,k+1-j] for j=1:k)
        x[:,1+k] .= rand_norm(cond_mean, sqrt(S2[1+k]), npath)

        # 不需要计算用n个预报n+1个的公式
        k == n-1 && break

        # 继续计算从前k+1个预报第k+2个的系数
        a2[1+k] = (gam_vec[1+(k+1)] - 
            sum(a1[1:k] .* gam_vec[(1+k):-1:(1+1)])) /
            S2[1+k]
        #println("a[",k+1,", ", k+1, "]=", a2[1+k])
        a2[1:k] .= a1[1:k] .- a2[1+k] .* a1[k:-1:1]

        # 重复利用a1, a2的空间，用a1表示当前可用的预报系数
        a1, a2 = a2, a1
    end # for k

    if p > 0 # 用p个可完全线性预报
        for k=(p+1):n 
            cond_mean = sum(a1[j] .* x[:,k-j] for j=1:p)
            x[:,k] .= rand_norm(cond_mean, sqrt(S2[1+p]), npath)
        end
    end # if p>0
    #show(round.(S2, digits=4))
    return x
end 

# 输入：自协方差函数gam(k)，以函数形式输入
# n是需要输出的长度
function gauss_levinson(gam::Function, 
    n::Int = 100, npath=1; eps=1e-6)
    gam_vec = gam.(0:(n-1))
    #@show gam_vec
    return gauss_levinson(gam_vec, npath; 
        eps = eps)
end

设某个高斯平稳列的协方差函数

{γk}满足

γk=exp⁡(−k2/(2ℓ2)),

其中

ℓ>0描述了前后相关性持续的时间长度，

ℓ越大相关性越强，
结果的轨道越光滑。
测试，生成5条轨道：

using CairoMakie
CairoMakie.activate!()
using Random
n, npath, ell = 200, 5, 5.0
Random.seed!(101)
x = gauss_levinson(k -> exp(-k^2 / (2*ell^2)), n, npath);
fig = series(1:n, x)

dc <- data.frame(
  x = c(c(13, 27, 33, 40, 43, 61, 125, 135, 161, 330),
        c(2, 3, 6, 7, 9, 12, 24, 35, 35, 74, 112, 122)),
  type = factor(rep(c("腺癌", "鳞癌"), c(10, 12)), 
                levels=c("腺癌", "鳞癌")))
ttc <- t.test(x ~ type, data = dc,
  var.equal=TRUE); ttc

using Combinatorics

# 两样本t统计量计算函数
function tstat(x, y)
    n1 = length(x)
    n2 = length(y)
    mx = mean(x)
    my = mean(y)
    S2p = 1/(n1+n2-2)*(
        sum((x .- mx) .^ 2) +
        sum((y .- my) .^ 2))
    tstat = (mx - my) / sqrt(S2p * (1/n1 + 1/n2))
    return tstat
end 

xv = [13, 27, 33, 40, 43, 61, 125, 135, 161, 330]
yv = [2, 3, 6, 7, 9, 12, 24, 35, 35, 74, 112, 122]

n1 = length(xv)
n2 = length(yv)
mixedxy = [xv; yv]

# 真实样本的统计量值
t0 = tstat(xv, yv)
# 所有的从n1+n2个样本点中取出n1个的罗列，
# 其中每一个是1:(n1+n2)的n1个下标的集合
allcomb = Combinatorics.combinations(1:(n1+n2), n1)
# 计算每一个的置换样本的统计量
tstats = map(allcomb) do ix 
    xv1 = mixedxy[ix]
    yv1 = mixedxy[setdiff(1:(n1+n2), ix)]
    tstat(xv1, yv1)
end
# 双侧检验p值
pvalue = 2*min(mean(tstats .> t0), 
    mean(tstats .< t0))
@show pvalue;
## pvalue = 0.0577472063540175

Random.seed!(101)
M = 10_000
tstats = map(1:M) do _
    rp = Random.randperm(n1+n2)
    xv1 = mixedxy[rp[1:n1]]
    yv1 = mixedxy[rp[n1+1:end]]
    tstat(xv1, yv1)
end
pvalue = 2*min(mean([t0; tstats] .>= t0), 
    mean([t0; tstats] .<= t0))
@show pvalue;
## pvalue = 0.0613938606139386