如何生成一个随机变量/随机向量的随机样本?
import random, math
from typing import List
import matplotlib.pyplot as plt
import seaborn as sns
import numpy as np
import pandas as pd
from mpl_toolkits.mplot3d import Axes3D
1. 连续型随机变量
在已知分布函数的表达式的情况下,有
X
:
=
F
−
1
(
U
)
∼
F
X:=F^{-1}(U)\sim F
X:=F−1(U)∼F.
以指数分布为例,
F
(
x
)
=
1
−
e
−
λ
x
(
x
≥
0
)
F(x)=1-e^{-\lambda x} (x \ge 0)
F(x)=1−e−λx(x≥0)
先生成一个
[
0
,
1
]
[0,1]
[0,1]之间的均匀分布随机数
r
r
r,再求出
F
(
x
∗
)
=
r
F(x^{*})=r
F(x∗)=r 的一个近似根,这个根就是我们要的指数分布随机变量的一个样本。
关于连续的单调函数的求根,二分法足够起到良好的效果。
def bisect_exp(lambda_, r: float) -> float:
"""用二分法求指数分布函数 F(x)=r 的根"""
F = lambda x: 1 - math.exp(-lambda_ * x) if x >= 0 else 0
lo, hi = 0, 1
while F(hi) < r:
hi *= 10
while (hi - lo) > 1e-5: # 设置求根误差限
mid = (lo + hi) / 2
if F(mid) > r:
hi = mid
else:
lo = mid
return (lo + hi) / 2
def random_exp(lambda_, size:int =100) -> List[float]:
"""生成长度为size的指数分布随机样本"""
res = []
for _ in range(size):
r = random.random()
res.append(bisect_exp(lambda_, r))
return res
以 E x p ( 2 ) \mathrm{Exp(2)} Exp(2) 看看效果,先画出该分布的密度函数 f ( x ) = 2 e − 2 x ( x ≥ 0 ) f(x)=2e^{-2x}(x\ge 0) f(x)=2e−2x(x≥0)
plt.figure(figsize=(12, 4))
u = np.linspace(0, 4, 10000); v = 2 * np.exp(-2 * u)
plt.subplot(121); plt.plot(u, v)
plt.subplot(122); l = random_exp(2, 10000); sns.distplot(l, kde=False, norm_hist=True)
plt.show()
两图对比,可以看到分布还是很接近的!
离散型随机变量
对于取值有限的离散型分布,其随机样本是容易模拟的。
考虑Bernoulli分布:
X
=
{
0
,
p
=
1
2
1
,
p
=
1
2
X=\left\{ \begin{aligned} 0,\;& p=\frac{1}{2}\\ 1,\;& p=\frac{1}{2} \end{aligned} \right.
X=⎩⎪⎨⎪⎧0,1,p=21p=21
直接生成
[
0
,
1
]
[0,1]
[0,1]之间的均匀分布的随机数,小于0.5记为0,大于0.5记为1,这里不做展示。
以下以 possion 分布为例
def random_possion(lambda_, size=100) -> int:
res = []
for _ in range(size):
r = random.random()
k, s = 0, math.exp(-lambda_)
while s < r:
k += 1
s += math.exp(-lambda_) * lambda_ ** k / math.factorial(k)
res.append(k)
return res
还是以 P ( 2 ) \mathrm{P(2)} P(2) 为例,先计算标准的概率分布
pd.DataFrame({
'prob': [math.exp(-2) * 2 ** k / math.factorial(k) for k in range(10)]
})
prob | |
---|---|
0 | 0.135335 |
1 | 0.270671 |
2 | 0.270671 |
3 | 0.180447 |
4 | 0.090224 |
5 | 0.036089 |
6 | 0.012030 |
7 | 0.003437 |
8 | 0.000859 |
9 | 0.000191 |
rp = random_possion(2, 10000)
pd.Series(rp).value_counts(True, False).sort_index()
0 0.1356
1 0.2661
2 0.2744
3 0.1801
4 0.0892
5 0.0370
6 0.0116
7 0.0046
8 0.0009
9 0.0004
11 0.0001
dtype: float64
从数值上看,与实际情况是接近的!
pd.Series(rp).value_counts(True, False).sort_index().plot.bar();plt.show()
随机向量
对于随机向量
(
X
,
Y
)
(X, Y)
(X,Y),
X
X
X 和
Y
Y
Y 可能不是独立的
考虑一个二元正态分布
(
X
,
Y
)
∼
N
(
μ
1
,
μ
2
,
σ
1
2
,
σ
2
2
,
ρ
)
(X, Y) \sim \mathcal{N}(\mu_1, \mu_2, \sigma^{2}_1, \sigma_2^{2}, \rho)
(X,Y)∼N(μ1,μ2,σ12,σ22,ρ),在概率论中我们学过二元正态分布的标准分解式,条件分布
Y
∣
X
=
x
Y|X=x
Y∣X=x 也是一个正态分布
边际分布
Y
∣
X
=
x
∼
N
(
μ
2
+
ρ
σ
2
σ
1
(
x
−
μ
1
)
,
σ
2
2
(
1
−
ρ
2
)
)
,
X
∣
Y
=
y
∼
N
(
μ
1
+
ρ
σ
1
σ
2
(
y
−
μ
2
)
,
σ
1
2
(
1
−
ρ
2
)
)
Y|X=x \sim \mathcal{N}\left(\mu_2+\rho \frac{\sigma_2}{\sigma_1}(x-\mu_1), \sigma_2^2(1-\rho^2)\right),\quad X|Y=y \sim \mathcal{N}\left(\mu_1+\rho\frac{\sigma_1}{\sigma_2}(y-\mu_2), \sigma_1^2(1-\rho^2)\right)
Y∣X=x∼N(μ2+ρσ1σ2(x−μ1),σ22(1−ρ2)),X∣Y=y∼N(μ1+ρσ2σ1(y−μ2),σ12(1−ρ2))
在这个情况下就可以使用 Gibbs 采样的方式得到一串随机采样,这串随机采样的极限分布就是原二元正态分布
random.normalvariate(mu, sigma) 返回均值为 mu, 标准差为 sigma 的一个随机正态样本
考虑 N ( 1 , 2 , 9 , 4 , 0.5 ) \mathcal{N}(1, 2, 9, 4, 0.5) N(1,2,9,4,0.5)
def random_norm(mu1, mu2, sigma1, sigma2, r, size=100) -> List[List[float]]:
res = []
chain = [mu1, mu2]
for i in range(1000+size): # 这里的 1000 是为了让 Markov 链趋向极限分布
# 计算给定 Y = chain[1] 时 X 的边际分布
y = chain[1]
chain[0] = random.normalvariate(
mu=mu1 + r * sigma1 / sigma2 * (y - mu2),
sigma=sigma1 ** 2 * (1 - r ** 2))
# 计算给定 X = chain[0] 时 Y 的边际分布
x = chain[0]
chain[1] = random.normalvariate(
mu=mu2 + r * sigma2 / sigma1 * (x - mu1),
sigma=sigma2 ** 2 * (1 - r ** 2))
if i >= 1000:
res.append(chain[:])
return res
sample = np.array(random_norm(1, 2, 3, 2, 0.5, size=10000))
画出样本的散点图
plt.scatter(sample[:,0], sample[:, 1], s=3)
plt.axis('equal'); plt.show()
fig = plt.figure()
ax = Axes3D(fig)
x = np.arange(-25, 25, 1)
y = np.arange(-15, 15, 0.5)
xn, yn = x.size, y.size
z = np.zeros((xn, yn))
for i in sample:
z[min(x.searchsorted(i[0]), xn-1)][min(y.searchsorted(i[1]), yn-1)] += 1
x, y = np.meshgrid(x, y)
ax.plot_surface(x, y, z.T, rstride=1, cstride=1, cmap='rainbow')
ax.set_xlabel('x'); ax.set_ylabel('y') ; plt.show()
观察样本散点图,可以看到在某个方向上,随机向量呈现出一种不相关的状态。事实上,由概率论的知识,我们知道
存在正交变换
Q
,
s
t
.
(
X
′
,
Y
′
)
=
Q
(
X
,
Y
)
Q, st. (X^{'}, Y^{'}) = Q (X, Y)
Q,st.(X′,Y′)=Q(X,Y),且
X
′
,
Y
′
X^{'}, Y^{'}
X′,Y′ 是不相关的。而正态变量的不相关性等价于独立性,我们可以生成
(
X
′
,
Y
′
)
(X^{'}, Y^{'})
(X′,Y′) 的随机样本再变换回去。这是根据多元正态分布的特点得到的模拟方法。
Markov 链的一个轨道与其极限分布的关系
事实上,对于一个比较“良好”的 离散Markov过程,通过对一个轨道进行采样,计算各个状态经过的频次比,是可以得到其平稳分布的近似的!
假定一个状态转移矩阵
[
P
11
P
12
P
13
P
21
P
22
P
23
P
31
P
32
P
33
]
=
[
0.2
0.4
0.4
0.3
0.1
0.6
0.7
0.2
0.1
]
\begin{matrix} \begin{bmatrix} P_{11} & P_{12} & P_{13} \\ P_{21} & P_{22} & P_{23} \\ P_{31} & P_{32} & P_{33} \end{bmatrix} =\begin{bmatrix} 0.2 & 0.4 & 0.4 \\ 0.3 & 0.1 & 0.6 \\ 0.7 & 0.2 & 0.1 \end{bmatrix} \end{matrix}
⎣⎡P11P21P31P12P22P32P13P23P33⎦⎤=⎣⎡0.20.30.70.40.10.20.40.60.1⎦⎤,有三个状态
{
1
,
2
,
3
}
\{1,2,3\}
{1,2,3}
初始状态为
1
1
1,考虑画出这条Markov过程一个轨道
m = np.array([[0.2, 0.4, 0.4],
[0.3, 0.1, 0.6],
[0.7, 0.2, 0.1]])
eiv = np.abs(np.real(np.linalg.eig(m.T)[1][:, 0]))
eiv /= sum(eiv)
eiv
array([0.39884393, 0.25433526, 0.34682081])
eiv 是状态转移矩阵的特征值 1 1 1 的左特征向量,代表这个马氏过程的平稳分布!
cumsum = np.cumsum(m, axis=1)
def transfer(cumsum: np.ndarray, state: int) -> int:
"""返回从状态 state 随机转移到的下一个状态"""
return cumsum[state-1].searchsorted(random.random()) + 1
现在记录一个长度 1000000 1000000 1000000 的轨道
state = 1
record = [1]
for _ in range(1000000):
state = transfer(cumsum, state)
record.append(state)
pd.Series(record).value_counts(True, False)
1 0.399224
2 0.254174
3 0.346603
dtype: float64
可以看到,三个状态经过的频次刚好是平稳分布的较好近似!进一步,如果要估计“用频次估计平稳分布”的好坏,可以继续研究这样子做的方差,进而得到相应平稳分布估计量的区间估计!