beta rejection is finished

Author: sijichun

code_in_notes/check_random.py

@@ -10,6 +10,6 @@
 # 画图
 import matplotlib.pyplot as plt
 # 设定图像大小
-plt.rcParams['figure.figsize'] = (10.0, 10.0)
+plt.rcParams['figure.figsize'] = (8.0, 5.0)
 plt.scatter(x0,x1,color='blue') ## 画出散点图
 plt.savefig("check_random_py.eps")

code_in_notes/importance_sampling_beta.py

@@ -6,23 +6,45 @@
 # 设定参数
 alpha=4
 beta=2
+# beta分布密度函数
+f=lambda x: 1/scisp.beta(alpha,beta)* \
+    x**(alpha-1) * (1-x)**(beta-1)
 # 计算密度函数最大值
-x_star=(1-alpha)/(2-alpha-beta)
-M=1/scisp.beta(alpha,beta)*x_star**(alpha-1)*x_star**(beta-1)
+M=f((1-alpha)/(2-alpha-beta))
 print(M)
-# 随机抽两个均匀分布，一个为(0,1)，一个为(0,M)，抽300个
-N=300
-u1=nprd.random(N)
-u2=nprd.random(N)*M
+# 随机抽两个均匀分布，一个为(0,1)，一个为(0,M)，抽500个
+N=500
+x=nprd.random(N)
+u=nprd.random(N)*M
+# 挑出使得u<beta密度函数的x
+accepted=[i for i in range(N) if u[i]<=f(x[i])]
+rand_beta=x[accepted] #生成的Beta分布随机数
+rand_U_selected=u[accepted]
 # 画图
 import matplotlib.pyplot as plt
-# 设定图像大小
-plt.rcParams['figure.figsize'] = (10.0, 10.0)
+# 设定图像大小和坐标范围
+plt.rcParams['figure.figsize'] = (8.0, 5.0)
+plt.xlim(0,1)
+plt.ylim(0,M+0.1)
 # 横线和竖线
-x=
-# 标题
-plt.xlabel('$u_1$')
-plt.ylabel("$u_2$")
-plt.title('Relationship of x and y')
-plt.scatter(u1,u2,color='blue') ## 画出散点图
+x_grid=np.linspace(0,1,100)#(0,1)均匀的100个点
+y_M=np.ones(100)*M
+beta_dens=f(x_grid)
+y1=np.linspace(0,f(0.4),20)
+y2=np.linspace(f(0.4),M,20)
+x_hline=np.ones(20)*0.4
+plt.xlabel(r'$X$')
+plt.ylabel(r'$U$')
+plt.title('Sampling from Beta dist.') # 标题
+## 画出散点图
+plt.scatter(x,u,color='blue',s=0.8)
+plt.scatter(rand_beta,rand_U_selected,color='black',s=0.8)
+plt.plot(x_grid,y_M,color='red') ## 最大值
+plt.plot(x_grid,beta_dens,color='orange') ## 密度函数
+plt.plot(x_hline,y1,color='grey') ## 接受区域
+plt.text(0.42,0.3,"Acceptance",fontsize=15,
+    horizontalalignment="left")
+plt.plot(x_hline,y2,color='pink') ## 拒绝区域
+plt.text(0.38,1.5,"Rejection",fontsize=15,
+    horizontalalignment="right")
 plt.savefig("importance_sampling_beta.eps")

code_in_notes/rejection.py

@@ -24,19 +24,17 @@
 sample2=[] #结果
 # 计算参数
 lam=(c+np.sqrt(c**2+4))/2
-M=np.exp((lam**2-2*lam*c)/2)/(np.sqrt(2*np.pi)*
+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) #正态分布密度函数前面的常数
 while accept<N:
     ## 产生指数分布
     x=-1*np.log(nprd.uniform())/lam+c
     ## 接受概率
-    r=normal_m*np.exp(-1*(x-c)**2/2)/(M*
+    r=normal_m*np.exp(-1*(x-c)**2/2)/(M* \
         lam*np.exp(-1*lam*(x-c)))
     if nprd.uniform()<r:
         sample2.append(x)
         accept=accept+1
     times=times+1
 print("接受率=",N/times)
-print(sample1)
-print(sample2)