在Python中使用逆變換方法生成隨機(jī)變量

作者：deephub 2021-03-18 10:14:06

在仿真理論中，生成隨機(jī)變量是最重要的“構(gòu)建塊”之一，而這些隨機(jī)變量大多是由均勻分布的隨機(jī)變量生成的。其中一種可以用來產(chǎn)生隨機(jī)變量的方法是逆變換法。在本文中，我將向您展示如何使用Python中的逆變換方法生成隨機(jī)變量(包括離散和連續(xù)的情況)。

目標(biāo)

[[388329]]

概念

給定隨機(jī)變量U，其中U在(0,1)中均勻分布。假設(shè)我們要生成隨機(jī)變量X，其中累積分布函數(shù)(CDF)為：

逆變換方法的思想是通過如下使用其逆CDF從任何概率分布中生成一個(gè)隨機(jī)數(shù)。

對(duì)于離散隨機(jī)變量，步驟略有不同。假設(shè)我們想生成一個(gè)離散隨機(jī)變量X的值，它具有一個(gè)概率質(zhì)量函數(shù)(PMF)

為了生成X的值，需要生成一個(gè)隨機(jī)變量U，U在(0,1)中均勻分布，并且定義

通過以上步驟，我們可以按如下方法創(chuàng)建逆變換方法的算法。

連續(xù)隨機(jī)數(shù)代碼實(shí)現(xiàn)

首先，我們實(shí)現(xiàn)此方法以生成連續(xù)隨機(jī)變量。假設(shè)我們要模擬一個(gè)隨機(jī)變量X，該變量遵循均值λ(即X〜EXP(λ))的指數(shù)分布。我們知道指數(shù)分布的概率分布函數(shù)(PDF)是

CDF如下

然后，我們可以使用以下的方法寫出逆CDF

在Python中，我們可以通過如下編寫這些代碼行來簡單地實(shí)現(xiàn)它。

### Generate exponential distributed random variables given the mean  
### and number of random variables  
def exponential_inverse_trans(n=1,mean=1):  
U=uniform.rvs(size=n)  
X=-mean*np.log(1-U)  
actual=expon.rvs(size=n,scale=mean)  
 
plt.figure(figsize=(12,9))  
plt.hist(X, bins=50, alpha=0.5, label="Generated r.v.")  
plt.hist(actual, bins=50, alpha=0.5, label="Actual r.v.")  
plt.title("Generated vs Actual %i Exponential Random Variables" %n)  
plt.legend()  
plt.show()  
return X

我們可以通過運(yùn)行以下示例來嘗試上面的代碼。請(qǐng)注意，由于我們要生成隨機(jī)變量，因此結(jié)果可能會(huì)有所不同。

cont_example1=exponential_inverse_trans(n=100,mean=4)  
cont_example2=exponential_inverse_trans(n=500,mean=4)  
cont_example3=exponential_inverse_trans(n=1000,mean=4)

在Python中使用逆變換方法生成隨機(jī)變量

看起來很有趣。如果將其與實(shí)際變量進(jìn)行比較，我們可以看到生成的隨機(jī)變量具有非常相似的結(jié)果。可以調(diào)整均值(請(qǐng)注意，我為expon.rvs()函數(shù)定義的均值是指數(shù)分布中的比例參數(shù))和/或生成的隨機(jī)變量的數(shù)量，以查看不同的結(jié)果。

離散隨機(jī)數(shù)實(shí)現(xiàn)代碼

對(duì)于離散隨機(jī)變量情況，假設(shè)我們要模擬遵循以下分布的離散隨機(jī)變量情況X

首先，我們編寫函數(shù)以使用這些代碼行為一個(gè)樣本生成離散隨機(jī)變量。

### Generate arbitary discrete distributed random variables given  
### the probability vector  
def discrete_inverse_trans(prob_vec):  
U=uniform.rvs(size=1)  
if U<=prob_vec[0]:  
return 1  
else:  
for i in range(1,len(prob_vec)+1):  
if sum(prob_vec[0:i])<U and sum(prob_vec[0:i+1])>U:  
return i+1

然后，我們創(chuàng)建一個(gè)函數(shù)以使用這些代碼行生成許多隨機(jī)變量樣本。

def discrete_samples(prob_vec,n=1):  
sample=[]  
for i in range(0,n):  
sample.append(discrete_inverse_trans(prob_vec))  
return np.array(sample)

最后，我們創(chuàng)建一個(gè)函數(shù)來模擬結(jié)果，并通過這些代碼行將其與實(shí)際結(jié)果進(jìn)行比較。

def discrete_simulate(prob_vec,numbers,n=1):  
sample_disc=discrete_samples(prob_vec,n)  
unique, counts=np.unique(sample_disc,return_counts=True)  
 
fig=plt.figure()  
ax=fig.add_axes([0,0,1,1])  
prob=counts/n  
ax.bar(numbers,prob)  
ax.set_title("Simulation of Generating %i Discrete Random Variables" %n)  
plt.show()  
 
data={'X':unique,'Number of samples':counts,'Empirical Probability':prob,'Actual Probability':prob_vec}  
df=pd.DataFrame(data=data)  
return df

我們可以在下面運(yùn)行一些示例以查看結(jié)果。同樣，請(qǐng)注意，由于我們要生成隨機(jī)變量，因此結(jié)果可能會(huì)有所不同。

prob_vec=np.array([0.1,0.3,0.5,0.05,0.05])  
numbers=np.array([1,2,3,4,5])  
 
dis_example1=discrete_simulate(prob_vec, numbers, n=100)  
dis_example2=discrete_simulate(prob_vec, numbers, n=500)  
dis_example3=discrete_simulate(prob_vec, numbers, n=1000)

在Python中使用逆變換方法生成隨機(jī)變量

In[11]: dis_example1  
Out[11]:  
X Number of samples Empirical Probability Actual Probability  
0 1 8 0.08 0.10  
1 2 35 0.35 0.30  
2 3 50 0.50 0.50  
3 4 5 0.05 0.05  
4 5 2 0.02 0.05In[12]: dis_example2  
Out[12]:  
X Number of samples Empirical Probability Actual Probability  
0 1 53 0.106 0.10  
1 2 159 0.318 0.30  
2 3 234 0.468 0.50  
3 4 30 0.060 0.05  
4 5 24 0.048 0.05In[13]: dis_example3  
Out[13]:  
X Number of samples Empirical Probability Actual Probability  
0 1 108 0.108 0.10  
1 2 290 0.290 0.30  
2 3 491 0.491 0.50  
3 4 51 0.051 0.05  
4 5 60 0.060 0.05