MC采样时是否使用while true——后附代码

理论篇

1.理论基础：while true:直到循环内的条件成立才停止,可能会循环多次，直到条件成立才退出循环。

2.应用：while true在采样中的应用如下：

2.1使用 while True:直到条件成立才停止

如下述代码。会循环10000+次，结果会产生10000个成功的采样点

def acceptreject(split_value):
    global t
    global power
    global sum_
    global c
    while True:
        x=random.uniform(0,1)
        y=random.uniform(0,1)
        if y*c<=math.pow((x-0.4),4)/sum_:
            return x
samples=[]
for i in range(10000):
    samples.append(AcceptReject(t))

2.2不使用 while True:只运行一次

如下述代码。在条件不成立的情况下可能返回None。（一般需要判断返回值是否为空 if data is not None:）会循环10000次，产生<10000个结果。

def acceptreject1(split_value):
    global t
    global power
    global sum_
    global c
    x=random.uniform(0,1)
    y=random.uniform(0,1)
    if y*c<=math.pow((x-0.4),7)/sum_:
        return x
    else:
        return 
samples=[]
for i in range(10000):
    si=acceptreject1(t)
    if si is not None:
        samples.append(si)

3.实验对比

方法一：设置拒绝接受域时不使用while true采样效果图

方法二：设置拒绝接受域时使用while true采样效果图

方法一与方法二对比采样效果图

实验平台：jupyter

#f(x)=(x-0.4)^0.4/0.176
import matplotlib.pyplot as plt
import numpy as np
import math
import random

#使采样可以复现
import os
seed = 13
random.seed(seed)
np.random.seed(seed)
os.environ['PYTHONHASHSEED'] = str(seed) # 为了禁止hash随机化，使得实验可复现。

def acceptreject(split_value):
    global t
    global power
    global sum_
    global c
    while True:
        x=random.uniform(0,1)
        y=random.uniform(0,1)
        if y*c<=math.pow((x-0.4),4)/sum_:
            return x

def acceptreject1(split_value):
    global t
    global power
    global sum_
    global c
    x=random.uniform(0,1)
    y=random.uniform(0,1)
    if y*c<=math.pow((x-0.4),4)/sum_:
        return x
    else:
        return         
        
t=0.4  
power=4 
sum_=(math.pow((1-t),power+1)-math.pow((0-t),power+1))/(power+1)
c=math.pow((1-t),power)/sum_

x=np.linspace(0,1,100)
cc=[c for i in x]
fx=[math.pow((xi-t),power)/sum_ for xi in x]


#调用acceptreject(split_value)采样
samples=[acceptreject(t) for i in range(10000)]
print("使用while true时有效采样个数:{0}".format(len(samples)))
#print(samples1)

#调用acceptreject1(split_value)采样
samples1=[]
for i in range(10000):
    si=acceptreject1(t)
    if si is not None:
        samples1.append(si)
print("不使用while true时有效采样个数:{0}".format(len(samples1)))
#print(samples1)

plt.plot(x,cc,'--',label='c:f(x)_max=(1-0.4)^4/0.0176')
plt.plot(x,fx,label='f(x)_max=(x-0.4)^4/f_0_1_(x-0.4)^4')
plt.hist(samples1,bins=50,density=True,label='samples< distribution density',color='red') 
plt.hist(samples,bins=50,density=True, label='samples  distribution density')   
plt.legend(loc="upper right",title='function',color='green')
plt.show()

参考文章
1.While(true)无限循环
2.python代码if not x: 和if x is not None:和if not x is None:使用
3.马尔可夫链蒙特卡罗算法（MCMC）
4.MCMC(一)蒙特卡罗方法