forked from sijichun/MathStatsCode
-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathrejection.py
More file actions
44 lines (43 loc) · 1.21 KB
/
rejection.py
File metadata and controls
44 lines (43 loc) · 1.21 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
#!/usr/bin/python3
## file: rejection.py
import numpy as np
import numpy.random as nprd
import scipy.special as scisp
## 设定参数
c=3 #抽取大于5的正态分布样本
N=1000 #抽取200个
## 直接使用正态分布进行拒绝抽样
times=0 #产生了多少次正态分布
accept=0 #接受了多少个
sample1=[] #结果
while accept<N:
x=nprd.normal()
if x>=c:
sample1.append(x)
accept=accept+1
times=times+1
print("接受率=",N/times)
## 使用指数分布进行拒绝抽样
times=0 #产生了多少次指数分布
accept=0 #接受了多少个
sample2=[] #结果
# 计算参数
lam=(c+np.sqrt(c**2+4))/2
M=np.exp((lam**2-2*lam*c)/2)/(np.sqrt(2*np.pi)* \
lam*(1-scisp.ndtr(c)))
normal_m=1/np.sqrt(2*np.pi) #正态分布密度函数前面的常数
# 截断正态分布密度函数
f_trunc_normal = lambda x: \
normal_m*np.exp(-0.5*x**2)/(1-scisp.ndtr(c))
# 指数分布密度函数
f_exponential = lambda x: lam*np.exp(-1*lam*(x-c))
while accept<N:
## 产生指数分布
x=c-np.log(nprd.uniform())/lam
## 接受概率
r=f_trunc_normal(x)/(M*f_exponential(x))
if nprd.uniform()<=r:
sample2.append(x)
accept=accept+1
times=times+1
print("接受率=",N/times)