跳转至

舍选法生成随机变量

针对离散变量，使用判别方法生成随机变量。针对连续变量，使用舍选法生成随机变量。

unnamed-chunk-4-1

分布形状：\(\Gamma(2,1)\) 随机数

unnamed-chunk-8-3

舍选法的理论推导

Python
import random
import numpy as np
import math
import matplotlib.pyplot as plt

plt.rcParams["axes.unicode_minus"] = False
plt.rcParams["font.sans-serif"] = ["SimHei"]
# 字号大小
plt.rcParams["font.size"] = 16
# 渲染公式
plt.rcParams["text.usetex"] = True

离散变量：生成 1 个二项分布随机数的算法

  1. 生成在 0 到 1 上均匀分布的随机数\(U \sim U(0,1)\)，且令\(t=1\)；

  2. 判断\(U\)是否小于\(0.7\)，若是，则\(X_t=1\)，否则\(X_t=0\)，并令\(t=t+1\)；

  3. 若\(t>n\)，则退出循环，返回\(\sum_{t=1}^{n} X_t\)，否则重复步骤 1 和 2；

Python
def random_binomial(n, p):
    success = 0
    for i in range(n):
        if random.uniform(0, 1) <= p:
            success += 1
    return success

生成 500 个服从 b(100, 0.7) 的二项分布随机数

Python
number_of_random_variavles = 500

n = 100
p = 0.7

random_binomial_list = []
for i in range(number_of_random_variavles):
    random_binomial_list.append(random_binomial(n, p))

绘制随机数的柱状图

Python
fig = plt.figure(figsize=(10, 6))
ax = fig.add_subplot(111)
ax.bar(range(0, n + 1), [random_binomial_list.count(i) for i in range(0, n + 1)])
Text Only
## <BarContainer object of 101 artists>
Python
ax.set_xlabel(r"$Binomial(100, 0.7)$")
ax.set_ylabel("Count", usetex=False)

unnamed-chunk-4-1

连续变量：舍选法

用舍选法生成 500 个取值大于 5 的\(\Gamma(2, 1)\)随机数，其密度函数为：

\[ p(x)=\frac{x e^{-x} 1\{x \geq 5\}}{\int_5^{\infty} u e^{-u} d u} \]

定义函数，生成指数分布随机数

Python
def rondom_exponential(lambd):
    return -math.log(1 - random.uniform(0, 1)) / lambd + 5

定义函数，生成\(\Gamma\)分布随机数

Python
def rondom_gamma():
    iteration = 0
    while True:
        iteration += 1
        x = rondom_exponential(0.5)
        u = random.uniform(0, 1)
        h_x = (1 / 5) * math.exp((5 - x) / 2) * x
        if u <= h_x:
            return x, iteration

生成取值大于 5 的\(\Gamma(2,1)\)随机数

Python
number_of_random_variavles = 500

alpha = 2
beta = 1

rondom_gamma_list = []
total_iteration = 0
for i in range(number_of_random_variavles):
    rondom_gamma_list.append(rondom_gamma()[0])
    total_iteration += rondom_gamma()[1]

绘制随机数的柱状图

Python
fig = plt.figure(figsize=(10, 6))
ax = fig.add_subplot(111)
ax.hist(rondom_gamma_list, bins=50)
Text Only
## (array([76., 44., 41., 37., 39., 27., 29., 28., 19., 21., 12., 17., 13.,
##         8.,  9.,  8.,  6.,  7., 10.,  5.,  5.,  3.,  4.,  1.,  3.,  3.,
##         3.,  6.,  1.,  2.,  2.,  2.,  1.,  2.,  1.,  0.,  1.,  1.,  0.,
##         0.,  0.,  0.,  0.,  0.,  1.,  1.,  0.,  0.,  0.,  1.]), array([ 5.00354213,  5.15295126,  5.3023604 ,  5.45176953,  5.60117866,
##         5.7505878 ,  5.89999693,  6.04940606,  6.19881519,  6.34822433,
##         6.49763346,  6.64704259,  6.79645172,  6.94586086,  7.09526999,
##         7.24467912,  7.39408825,  7.54349739,  7.69290652,  7.84231565,
##         7.99172479,  8.14113392,  8.29054305,  8.43995218,  8.58936132,
##         8.73877045,  8.88817958,  9.03758871,  9.18699785,  9.33640698,
##         9.48581611,  9.63522524,  9.78463438,  9.93404351, 10.08345264,
##        10.23286177, 10.38227091, 10.53168004, 10.68108917, 10.83049831,
##        10.97990744, 11.12931657, 11.2787257 , 11.42813484, 11.57754397,
##        11.7269531 , 11.87636223, 12.02577137, 12.1751805 , 12.32458963,
##        12.47399876]), <BarContainer object of 50 artists>)
Python
ax.set_xlabel(r"$Gamma(2, 1)$")
# x 轴刻度
ax.set_xticks(np.arange(5, 15, 1))
ax.set_ylabel("Count", usetex=False)

unnamed-chunk-8-3

输出总迭代次数和实际接受率

Python
print("总迭代次数：", total_iteration)
Text Only
## 总迭代次数： 789
Python
print("实际接受率：{:.2%}".format(number_of_random_variavles / total_iteration))
Text Only
## 实际接受率：63.37%

评论