用Python生成随机样本

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()

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()]
})

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

pd.Series(rp).value_counts(True, False).sort_index().plot.bar();plt.show()

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=))

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()

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])

cumsum = np.cumsum(m, axis=)
def transfer(cumsum: np.ndarray, state: int) -> int:
    """返回从状态 state 随机转移到的下一个状态"""
    return cumsum[state-1].searchsorted(random.random()) + 

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