Python+Scipy實現(xiàn)自定義任意的概率分布

更新時間：2022年08月10日 15:40:34 作者：nejssd

Scipy自帶了多種常見的分布，如正態(tài)分布、均勻分布、二項分布、多項分布、伽馬分布等等，還可以自定義任意的概率分布。本文將為大家介紹如何利用Scipy自定義任意的概率分布，感興趣的可以了解下

連續(xù)變量分布

考慮連續(xù)變量x滿足如下概率密度分布函數(shù):

其在實數(shù)域積分為1。我們可以通過scipy.stats中的rv_continuous類去實現(xiàn)這個分布，代碼如下:

from scipy.stats import rv_continuous
import matplotlib.pyplot as plt
import numpy as np
class MyDistribution(rv_continuous):
    def _pdf(self, x):#概率密度分布函數(shù)
        return 2*sqrt(0.1)*exp(-0.1*x**2)*cos(x)**2/(sqrt(pi)*(exp(-10) + 1))
distribution = MyDistribution()
xlist=np.linspace(-8,8,300)
ylist=distribution.pdf(xlist)
samples=distribution.rvs(size=200);#取200次樣

fig,ax=plt.subplots(figsize=(8,6))
ax.plot(xlist,ylist,lw=3,color='red',label="$\mathrm{ideal}$");
ax.hist(samples,color='blue',density=True, bins=np.arange(-8,8,0.25), histtype='barstacked', rwidth=0.9,label=r"$\mathrm{samples}$")
ax.legend(fontsize=20);
ax.set_xlabel(r"$x$",size=25)
ax.set_ylabel(r"$\mathrm{PDF}$",size=20)
ax.set_xlim(-8,8);
ax.tick_params(axis='both',direction='in',width=1.3,length=3,top=1,right=1,labelsize=20,pad=2)
fig.tight_layout();
fig.show();

運行結(jié)果如下：

增加采樣次數(shù)，分布直方圖逐漸趨于理想的概率分布函數(shù)P(x)。

離散變量分布

考慮連續(xù)變量x滿足泊松分布，則可以用scipy.stats中的rv_discrete類去實現(xiàn)這個分布，代碼如下:

from scipy.stats import rv_discrete
import matplotlib.pyplot as plt
import numpy as np
from scipy.special import factorial
class MyDistribution(rv_discrete):
    def _pmf(self, k, mu):
        return exp(-mu)*mu**k/factorial(k)
distribution = MyDistribution()
mu=2
samples=distribution.rvs(size=500,mu=mu);#取500次樣
klist = np.arange(0,10,1)
plist = distribution.pmf(klist,mu)
fig, ax = plt.subplots()
ax.plot(klist, plist, 'ro', ms=12, mec='r',label="$\mathrm{ideal}$");
ax.hist(samples,color='blue',density=True, bins=klist, histtype='barstacked', rwidth=0.8,label=r"$\mathrm{samples}$",align="left")
ax.legend(fontsize=20);
fig.show();

運行結(jié)果如下：

可以修改上述MyDistribution類中的pmf函數(shù)，實現(xiàn)任意想要的離散分布。

二項分布Binomial Distribution

是n個獨立的成功/失敗試驗中成功的次數(shù)的離散概率分布，其中每次試驗的成功概率為p。這樣的單次成功/失敗試驗又稱為伯努利試驗。實際上，當(dāng)n=1時，二項分布就是伯努利分布。

'''1、定義隨機變量'''
# 比如5次擲硬幣實驗，正面朝上的次數(shù)
n2=5
x2=np.arange(1,n2+1,1)
x2

array([1, 2, 3, 4, 5])

'''2、求對應(yīng)的概率質(zhì)量函數(shù) (PMF)'''
p2=0.5
pList2=stats.binom.pmf(x2,n2,p2)
# 返回一個列表，列表中每個元素表示隨機變量中對應(yīng)值的概率
pList2

array([0.15625, 0.3125 , 0.3125 , 0.15625, 0.03125])

'''3、繪圖'''
fig=plt.figure()
# plot在此的作用為顯示兩個標(biāo)記點
plt.plot(x2,pList2,marker='o',linestyle='None')
'''
vlines用于繪制豎直線(vertical lines),
參數(shù)說明：vline(x坐標(biāo)值, y坐標(biāo)最小值, y坐標(biāo)值最大值)
'''
plt.vlines(x2, 0, pList2)
plt.xlabel('隨機變量：拋硬幣5次')
plt.ylabel('概率')
plt.title('二項分布：n=%d,p2=%0.2f' % (n2,p2))
plt.show()

幾何分布Geometric Distribution

在n次伯努利試驗中，試驗k次才得到第一次成功的機率。詳細地說，是：前k-1次皆失敗，第k次成功的概率。幾何分布是帕斯卡分布當(dāng)r=1時的特例。

'''1、定義隨機變量'''
# 比如射箭1次中靶的概率為90%，射5次箭
k=5
x3=np.arange(1,k+1,1)
x3

array([1, 2, 3, 4, 5])

'''2、求對應(yīng)的概率質(zhì)量函數(shù) (PMF)'''
p3=0.7
pList3=stats.geom.pmf(x3,p3)
# 返回一個列表，表示在第i次射擊中，第一次射中的概率
pList3

array([0.7    , 0.21   , 0.063  , 0.0189 , 0.00567])

'''3、繪圖'''
fig=plt.figure()
# plot在此的作用為顯示兩個標(biāo)記點
plt.plot(x3,pList3,marker='o',linestyle='None')
'''
vlines用于繪制豎直線(vertical lines),
參數(shù)說明：vline(x坐標(biāo)值, y坐標(biāo)最小值, y坐標(biāo)值最大值)
'''
plt.vlines(x3, 0, pList3)
plt.xlabel('隨機變量：射擊5次')
plt.ylabel('概率')
plt.title('幾何分布：n=%d,p=%0.2f' % (k,p3))
plt.show()

泊松分布Poisson Distribution

描述在某單位時間內(nèi)，事件發(fā)生n次的概率

'''1、定義隨機變量'''
# 某機器每季度發(fā)生故障平均為1次，那么在一年中機器發(fā)生10次的概率為
mu=4 # 平均值
k=10 # 要求發(fā)生10次的概率
x4=np.arange(1,k+1,1)
x4

array([ 1,  2,  3,  4,  5,  6,  7,  8,  9, 10])

'''2、求對應(yīng)的概率質(zhì)量函數(shù) (PMF)'''
pList4=stats.poisson.pmf(x4,mu) # 一年的平均值為4
# 返回一個列表，表示1年中發(fā)生i次故障的概率
pList4

array([0.07326256, 0.14652511, 0.19536681, 0.19536681, 0.15629345,
       0.10419563, 0.05954036, 0.02977018, 0.01323119, 0.00529248])

'''3、繪圖'''
fig=plt.figure()
# plot在此的作用為顯示兩個標(biāo)記點
plt.plot(x4,pList4,marker='o',linestyle='None')
'''
vlines用于繪制豎直線(vertical lines),
參數(shù)說明：vline(x坐標(biāo)值, y坐標(biāo)最小值, y坐標(biāo)值最大值)
'''
plt.vlines(x4, 0, pList4)
plt.xlabel('隨機變量：發(fā)生k次故障')
plt.ylabel('概率')
plt.title('泊松分布：n=%d' % k)
plt.show()