如何生成一个随机变量/随机向量的随机样本？连续型随机变量离散型随机变量随机向量Markov 链的一个轨道与其极限分布的关系
如何生成一个随机变量/随机向量的随机样本？
代码语言：javascript复制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
连续型随机变量
def bisect_exp(lambda_, r: float) -> float:
"""用二分法求指数分布函数 F(x)=r 的根"""
F = lambda x: - math.exp(-lambda_ * x) if x >= else
lo, hi = ,
while F(hi) < r:
hi *=
while (hi - lo) > 1e-5: # 设置求根误差限
mid = (lo hi) /
if F(mid) > r:
hi = mid
else:
lo = mid
return (lo hi) /
def random_exp(lambda_, size:int =) -> List[float]:
"""生成长度为size的指数分布随机样本"""
res = []
for _ in range(size):
r = random.random()
res.append(bisect_exp(lambda_, r))
return res
plt.figure(figsize=(, ))
u = np.linspace(, , ); v = * np.exp(-2 * u)
plt.subplot(); plt.plot(u, v)
plt.subplot(); l = random_exp(, ); sns.distplot(l, kde=False, norm_hist=True)
plt.show()
两图对比，可以看到分布还是很接近的！
离散型随机变量
直接生成之间的均匀分布的随机数，小于0.5记为0，大于0.5记为1，这里不做展示。 以下以 possion 分布为例
代码语言：javascript复制def random_possion(lambda_, size=) -> int:
res = []
for _ in range(size):
r = random.random()
k, s = , math.exp(-lambda_)
while s < r:
k =
s = math.exp(-lambda_) * lambda_ ** k / math.factorial(k)
res.append(k)
return res
pd.DataFrame({
'prob': [math.exp(-2) * ** k / math.factorial(k) for k in range()]
})
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(, )
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
从数值上看，与实际情况是接近的！
代码语言：javascript复制pd.Series(rp).value_counts(True, False).sort_index().plot.bar();plt.show()
随机向量
random.normalvariate(mu, sigma) 返回均值为 mu, 标准差为 sigma 的一个随机正态样本
考虑
代码语言：javascript复制def random_norm(mu1, mu2, sigma1, sigma2, r, size=) -> List[List[float]]:
res = []
chain = [mu1, mu2]
for i in range( size): # 这里的 1000 是为了让 Markov 链趋向极限分布
# 计算给定 Y = chain[1] 时 X 的边际分布
y = chain[]
chain[] = random.normalvariate(
mu=mu1 r * sigma1 / sigma2 * (y - mu2),
sigma=sigma1 ** * ( - r ** ))
# 计算给定 X = chain[0] 时 Y 的边际分布
x = chain[]
chain[] = random.normalvariate(
mu=mu2 r * sigma2 / sigma1 * (x - mu1),
sigma=sigma2 ** * ( - r ** ))
if i >= :
res.append(chain[:])
return res
sample = np.array(random_norm(, , , , 0.5, size=))
画出样本的散点图
代码语言：javascript复制plt.scatter(sample[:,], sample[:, ], s=)
plt.axis('equal'); plt.show()
fig = plt.figure()
ax = Axes3D(fig)
x = np.arange(-25, , )
y = np.arange(-15, , 0.5)
xn, yn = x.size, y.size
z = np.zeros((xn, yn))
for i in sample:
z[min(x.searchsorted(i[]), xn-1)][min(y.searchsorted(i[]), yn-1)] =
x, y = np.meshgrid(x, y)
ax.plot_surface(x, y, z.T, rstride=, cstride=, cmap='rainbow')
ax.set_xlabel('x'); ax.set_ylabel('y') ; plt.show()
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)[][:, ]))
eiv /= sum(eiv)
eiv
array([0.39884393, 0.25433526, 0.34682081])
eiv 是状态转移矩阵的特征值 的左特征向量，代表这个马氏过程的平稳分布！
代码语言：javascript复制cumsum = np.cumsum(m, axis=)
def transfer(cumsum: np.ndarray, state: int) -> int:
"""返回从状态 state 随机转移到的下一个状态"""
return cumsum[state-1].searchsorted(random.random())
现在记录一个长度 的轨道
代码语言：javascript复制state =
record = []
for _ in range():
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
可以看到，三个状态经过的频次刚好是平稳分布的较好近似！进一步，如果要估计“用频次估计平稳分布”的好坏，可以继续研究这样子做的方差，进而得到相应平稳分布估计量的区间估计！
