題目
用舍選法在計算機中產生概率密度函數為
\(f(x)=\left\{\begin{array}{l} \frac{12}{(3+2 \sqrt{3}) \pi}\left(\frac{\pi}{4}+\frac{2 \sqrt{3}}{3} \sqrt{1-x^{2}}\right), 0 \leq x \leq 1 \\ 0 \end{array}\right.\)
的100個隨機數，具體要求：
（1）[0，1]均勻分布隨機數用線性同余法產生，參數由自己確定，不能用計算語言已有的函數。
（2）對用線性同余法產生的（0，1）均勻隨機數進行均勻性和獨立性檢驗，檢驗樣本為100個。
（３）計算舍選法的理論上舍選效率和實際仿真的舍選效率。
程序與分析
函數定義
# 線性同余法產生n個服從(0,1)均勻分布的隨機數
import time
import numpy as np
from matplotlib import pyplot as plt
from matplotlib import rcParams
# 設置字體為楷體
rcParams['font.sans-serif'] = ['KaiTi']
# 函數avg_random(n)產生n個服從(0,1)均勻分布的隨機數
def avg_random(n):
a = 32310901 # 乘子
c = 1729 # 增量
M = pow(2, 32)
theats = []
x0 = time.time() # 以時間戳為種子
xi = x0
for i in range(0, n):
xi = (a * x0 + c) % M
theati = xi / M
theats.append(theati)
x0 = xi
theats = np.array(theats)
return theats
# 1. 對產生的隨機數進行粗略檢驗，如果偏差超過10%，提示生成的隨機數均勻性很差，不進一步進行檢驗
# 粗略檢驗函數rough_test()檢驗通過返回True，失敗返回False
def rough_test(theats):
mean = np.average(theats)
sec = np.average(theats * theats)
if np.abs(mean - 0.5) / 0.5 < 0.1 * 0.5 or np.abs(sec - (1 / 3)) < 0.1 * (1 / 3):
return True
else:
return False
# 2. 粗略檢驗通過后，進行頻率檢驗
# 頻率檢驗函數frequency_test()，參數"sec_num"為區間個數，默認值為"10"，檢驗通過(置信度取5%)返回True，失敗返回False
def frequency_test(theats, sec_num=10):
N = len(theats) # 抽樣值的個數
m = N / sec_num # 理論頻數
sec_val = np.zeros(sec_num)
# 計算每個區間的實際頻數
for theat in theats:
for i in range(1, 11):
if theat < (i / sec_num):
sec_val[i - 1] += 1
break
C = np.sum(np.square(sec_val - m) / m)
# 采用5%置信度的9自由度卡方分布上分位數為16.919
if C < 16.919:
return True
else:
return False
# 3. 獨立性檢驗--相關系數檢驗，返回距離為distance（默認值為5）的兩個隨機數的相關系數
def corr_coeffient_test(theats, distance=5):
distance = np.abs(int(distance)) # 均值限制為正整數
mean = np.average(theats) # 樣本均值
std = np.var(theats) # 樣本方差
sumval = 0
for i in range(0, len(theats) - distance):
temp = theats[i] * theats[i + distance]
sumval += temp
R = (sumval / (len(theats) - distance) - pow(mean, 2)) / (pow(std, 2))
return R
# 4. 獨立性檢驗--聯立表檢驗，參數"k"為正方形每個維度被分割的個數（默認為4），參數"e"為數據偏離位數（默認為5），val_alpha為5%置信度k^2-1自由度卡方分布上分位數(默認25)；隨機數在獨立性上95%可信則返回True，反之返回False
def simu_table_test(theats, k=4, e=5, val_alpha=25):
sec_val = np.zeros((k, k))
for i, j in zip(range(0, len(theats)), np.roll(range(0, len(theats)), -e)):
for x in range(1, k + 1):
if theats[i] < (x / k):
for y in range(1, k + 1):
if theats[j] < (y / k):
sec_val[x - 1][y - 1] += 1
break
break
C = 0
for i in range(0, k):
for j in range(0, k):
c = pow((sec_val[i][j] - np.sum(sec_val, 1)[i] * np.sum(sec_val, 0)[j] / len(theats)), 2) / (np.sum(sec_val, 1)[i] * np.sum(sec_val, 0)[j] / len(theats))
C += c
if C < val_alpha:
return True
else:
return False
# 密度函數
def fuc(x):
return 12 / ((3 + 2 * np.sqrt(3)) * np.pi) * (np.pi / 4 + 2 * np.sqrt(3) / 3 * np.sqrt(1 - pow(x, 2)))
# 函數rejection_method(R1,R2,n)產生n個服從題目密度函數分布的隨機數，R1、R2為兩個服從(0,1)均勻分布的隨機數，同時返回理論效率'theory_efficience'和仿真效率'emulational_efficience'
def rejection_method(randow_num_list, n):
a = 0#函數下界
b = 1#函數上界
M = fuc(0)#函數最大值
X = []#隨機數初始化
num = 0#隨機數數量初始化
step=0#
failnum = 0
addstep=0
while num <n:
if not step%10:
addstep+=1
R1 = randow_num_list[step % (len(randow_num_list))]
R2 = randow_num_list[(step+addstep) % (len(randow_num_list))]
x = a + (b - a) * R1
if M * R2 <= fuc(x):
num += 1
X.append(x)
else:
failnum += 1
step+=1
return {'obje_randow_list': np.array(X), 'theory_efficience': 1 / (M * (a + b)), 'emulational_efficience': n / (n + failnum)}
生成100個服從(0,1)均勻分布的隨機數並進行粗略檢驗
\(\frac{1}{N} \sum_{i=1}^{N} R_{i} \approx \frac{1}{2}\)
\(\frac{1}{N} \sum_{i=1}^{N} R_{i}^{2} \approx \frac{1}{3}\)
\(R_i\)為隨機數，N為隨機數個數,如果兩個都明顯不成立，就可否定隨機性不夠要求。
randow_num_list = avg_random(100)
if rough_test(randow_num_list):
print("粗略檢驗通過！\n隨機數組為:\n ", randow_num_list)
else:
print("粗略檢驗失敗！\n隨機數組為：\n", randow_num_list)
輸出結果：
粗略檢驗通過！
隨機數組為:
[0.85159487 0.38553153 0.09123632 0.80387689 0.66494775 0.86175315
0.74260715 0.99250419 0.76666027 0.18396305 0.75009426 0.5288669
0.11249721 0.23132686 0.37451616 0.54591484 0.46852607 0.31343018
0.36082217 0.4449084 0.11814833 0.89293472 0.48887375 0.33682473
0.50663072 0.99573468 0.56172571 0.67126407 0.82313645 0.30861287
0.83140746 0.05285592 0.52025646 0.97192002 0.69488506 0.26584835
0.71075032 0.31535233 0.97442643 0.90242584 0.87177382 0.60438133
0.45030455 0.70864482 0.57518814 0.09044136 0.74346058 0.10024626
0.0408114 0.13287688 0.85784692 0.03787268 0.51931731 0.20832441
0.3880822 0.50524611 0.04926207 0.8398687 0.49066355 0.29365035
0.31442561 0.63639672 0.56949927 0.62473007 0.41733016 0.6282564
0.31666614 0.45719453 0.05920186 0.28278621 0.21508996 0.52478958
0.19995624 0.21906851 0.01380596 0.86616696 0.03639796 0.88335764
0.27265405 0.98767091 0.86073049 0.7010772 0.01781981 0.25839967
0.25399892 0.11907349 0.70671217 0.91646761 0.12017356 0.95450994
0.29517175 0.27535272 0.36116569 0.77153018 0.3022805 0.32689434
0.75088046 0.36458859 0.78798403 0.93310726]
粗略檢驗通過，進行頻率檢驗
對於（0，1）均勻分布數分成n組\(\left(0, \quad \frac{1}{n}\right),\left(\frac{1}{n}, \frac{2}{n}\right), \cdots\left(\frac{n-1}{n}, 1\right)\),每組的理論頻數\(E_{i}=\frac{N}{n}\),\(N\)為隨機數個數，實際頻數為\(Q_i\)，統計量\(\bar{c}=\sum_{i=1}^{n} \frac{\left(Q_{i}-E_{i}\right)^{2}}{E_{i}^{2}}\)服從\(\chi_{n-1}^2\)分布，對於一定數量\(\frac{N}{n}\geqslant5\)，若采用5%的置信度，查\(\chi_{n-1}^2\)分布表，得\(\chi_{0.05(n-1)}^2\)，若\(\bar{c}<\chi_{0.05(n-1)}^2\)，認為這批隨機數在統計性能上是95%可信。
if frequency_test(randow_num_list):
print("頻率檢驗通過！")
else:
print("頻率檢驗失敗！")
輸出結果：
頻率檢驗通過！
獨立性檢驗--相關系數檢驗
給定\(N\)個服從相同分布的隨機數\(r_1,r_2,\dots,r_n\)，我們計算前后距離為\(j\)的兩個隨機樣本的相關系數\(R_j\)
R = corr_coeffient_test(randow_num_list, 1)
print("距離為1的相關系數為：", R)
R = corr_coeffient_test(randow_num_list)
print("距離為5的相關系數為：", R)
R = corr_coeffient_test(randow_num_list, 10)
print("距離為10的相關系數為：", R)
R = corr_coeffient_test(randow_num_list, 20)
print("距離為20的相關系數為：", R)
R = corr_coeffient_test(randow_num_list, 50)
print("距離為50的相關系數為：", R)
輸出結果：
距離為1的相關系數為： -1.0571511721357674
距離為5的相關系數為： -2.2119618825288287
距離為10的相關系數為： -1.1470855033688878
距離為20的相關系數為： -0.9568683998367279
距離為50的相關系數為： 0.4553730525624434
獨立性檢驗--聯立表檢驗
if simu_table_test(randow_num_list):
print("這100個隨機數在獨立性上95%是可信的！")
else:
print("沒有把握確定這100個隨機數在獨立性上95%是可信的！")
輸出結果：
這100個隨機數在獨立性上95%是可信的！
舍選法生成100個隨機數並進行效率求解
num = 100#舍選法隨機數的個數
objection = rejection_method(randow_num_list, num)
print("生成的",num,"個服從要求密度函數的隨機數列為：\n", objection['obje_randow_list'])
print("理論效率為：", objection['theory_efficience'])
print("仿真效率為：", objection['emulational_efficience'])
輸出結果：
生成的 100 個服從要求密度函數的隨機數列為：
[0.85159487 0.38553153 0.09123632 0.80387689 0.76666027 0.18396305
0.75009426 0.5288669 0.11249721 0.23132686 0.37451616 0.54591484
0.46852607 0.31343018 0.36082217 0.4449084 0.11814833 0.89293472
0.33682473 0.50663072 0.56172571 0.67126407 0.82313645 0.30861287
0.83140746 0.05285592 0.52025646 0.97192002 0.26584835 0.31535233
0.97442643 0.87177382 0.60438133 0.45030455 0.70864482 0.57518814
0.09044136 0.74346058 0.10024626 0.0408114 0.13287688 0.85784692
0.03787268 0.51931731 0.20832441 0.3880822 0.50524611 0.04926207
0.8398687 0.49066355 0.29365035 0.31442561 0.63639672 0.56949927
0.62473007 0.41733016 0.6282564 0.31666614 0.45719453 0.05920186
0.28278621 0.21508996 0.19995624 0.21906851 0.01380596 0.86616696
0.03639796 0.88335764 0.27265405 0.7010772 0.01781981 0.25839967
0.25399892 0.11907349 0.70671217 0.12017356 0.29517175 0.27535272
0.36116569 0.3022805 0.32689434 0.75088046 0.78798403 0.93310726
0.85159487 0.38553153 0.09123632 0.80387689 0.66494775 0.86175315
0.74260715 0.99250419 0.76666027 0.18396305 0.75009426 0.5288669
0.11249721 0.37451616 0.54591484 0.46852607]
理論效率為： 0.8722740982016987
仿真效率為： 0.8547008547008547
查看舍選法生成隨機數的擬合效果
randow_num_list_test=avg_random(100000)
obje_randow_list_test=rejection_method(randow_num_list_test, 10000)['obje_randow_list']
bins = np.arange(0, 1.1,0.1)
weights = np.ones_like(obje_randow_list_test) / len(obje_randow_list_test)
plt.hist(obje_randow_list_test, weights=weights, bins=bins)
samples = np.arange(0.05,1,0.1)
density = fuc(samples)
density_1 = density / sum(density)
plt.plot(samples, density_1)
plt.title('概率密度函數擬合圖')
plt.show()
概率分布擬合圖