python 計算概率密度、累計分布、逆函數(shù)的例子
計算概率分布的相關(guān)參數(shù)時,一般使用 scipy 包,常用的函數(shù)包括以下幾個:
pdf:連續(xù)隨機分布的概率密度函數(shù)
pmf:離散隨機分布的概率密度函數(shù)
cdf:累計分布函數(shù)
百分位函數(shù)(累計分布函數(shù)的逆函數(shù))
生存函數(shù)的逆函數(shù)(1 - cdf 的逆函數(shù))
函數(shù)里面不僅能跟一個數(shù)據(jù),還能跟一個數(shù)組。下面用正態(tài)分布舉例說明:
>>> import scipy.stats as st>>> st.norm.cdf(0) # 標(biāo)準正態(tài)分布在 0 處的累計分布概率值0.5>>> st.norm.cdf([-1, 0, 1])# 標(biāo)準正態(tài)分布分別在 -1, 0, 1 處的累計分布概率值array([0.15865525, 0.5, 0.84134475])>>> st.norm.pdf(0) # 標(biāo)準正態(tài)分布在 0 處的概率密度值0.3989422804014327>>> st.norm.ppf(0.975)# 標(biāo)準正態(tài)分布在 0.975 處的逆函數(shù)值1.959963984540054>>> st.norm.lsf(0.975)# 標(biāo)準正態(tài)分布在 0.025 處的生存函數(shù)的逆函數(shù)值1.959963984540054
對于非標(biāo)準正態(tài)分布,通過更改參數(shù) loc 與 scale 來改變均值與標(biāo)準差:
>>> st.norm.cdf(0, loc=2, scale=1) # 均值為 2,標(biāo)準差為 1 的正態(tài)分布在 0 處的累計分布概率值0.022750131948179195
對于其他隨機分布,可能更改的參數(shù)不一樣,具體需要查官方文檔。下面我們舉一些常用分布的例子:
>>> st.binom.pmf(4, n=100, p=0.05) # 參數(shù)值 n=100, p=0.05 的二項分布在 4 處的概率密度值0.17814264156968956>>> st.geom.pmf(4, p=0.05) # 參數(shù)值 p=0.05 的幾何分布在 4 處的概率密度值0.04286875>>> st.poisson.pmf(2, mu=3) # 參數(shù)值 mu=3 的泊松分布在 2 處的概率密度值0.22404180765538775>>> st.chi2.ppf(0.95, df=10) # 自由度為 10 的卡方分布在 0.95 處的逆函數(shù)值18.307038053275146>>> st.t.ppf(0.975, df=10) # 自由度為 10 的 t 分布在 0.975 處的逆函數(shù)值2.2281388519649385>>> st.f.ppf(0.95, dfn=2, dfd=12) # 自由度為 2, 12 的 F 分布在 0.95 處的逆函數(shù)值3.8852938346523933
補充拓展:給定概率密度,生成隨機數(shù) python實現(xiàn)
實現(xiàn)的方法可以不止一種:
rejection sampling
invert the cdf
Metropolis Algorithm (MCMC)
本篇介紹根據(jù)累積概率分布函數(shù)的逆函數(shù)(2:invert the CDF)生成的方法。
自己的理解不一定正確,有錯誤望指正。
目標(biāo):
已知 y=pdf(x),現(xiàn)想由給定的pdf, 生成對應(yīng)分布的x
PDF是概率分布函數(shù),對其積分或者求和可以得到CDF(累積概率分布函數(shù)),PDF積分或求和的結(jié)果始終為1
步驟(具體解釋后面會說):
1、根據(jù)pdf得到cdf
2、由cdf得到inverse of the cdf
3、對于給定的均勻分布[0,1),帶入inverse cdf,得到的結(jié)果即是我們需要的x
求cdf逆函數(shù)的具體方法:
對于上面的第二步,可以分成兩類:
1、當(dāng)CDF的逆函數(shù)好求時,直接根據(jù)公式求取,
2、反之當(dāng)CDF的逆函數(shù)不好求時,用數(shù)值模擬方法
自己的理解:為什么需要根據(jù)cdf的逆去獲得x?
原因一:
因為cdf是單調(diào)函數(shù)因此一定存在逆函數(shù)(cdf是s型函數(shù),而pdf則不一定,例如正態(tài)分布,不單調(diào),對于給定的y,可能存在兩個對應(yīng)的x,就不可逆)
原因二:
這僅是我自己的直觀理解,根據(jù)下圖所示(左上為pdf,右上為cdf)
由步驟3可知,我們首先生成[0,1)的均勻隨機數(shù),此隨機數(shù)作為cdf的y,去映射到cdf的x(若用cdf的逆函數(shù)表示則是由x映射到y(tǒng)),可以參考上圖的右上,既然cdf的y是均勻隨機的,那么對于cdf中同樣范圍的x,斜率大的部分將會有更大的機會被映射,因為對應(yīng)的y范圍更大(而y是隨即均勻分布的),那么,cdf的斜率也就等同于pdf的值,這正好符合若x的pdf較大,那么有更大的概率出現(xiàn)(即重復(fù)很多次后,該x會出現(xiàn)的次數(shù)最多)
代碼實現(xiàn)——方法一,公式法
import numpy as npimport mathimport randomimport matplotlib.pyplot as pltimport collectionscount_dict = dict()bin_count = 20def inverseCDF(): ''' return the x value in PDF ''' uniform_random = random.random() return inverse_cdf(uniform_random) def pdf(x): return 2 * x # cdf = x^2, 其逆函數(shù)很好求,因此直接用公式法def inverse_cdf(x): return math.sqrt(x)def draw_pdf(D):global bin_count D = collections.OrderedDict(sorted(D.items())) plt.bar(range(len(D)), list(D.values()), align=’center’) # 因為映射bin的時候采用的floor操作,因此加上0.5 value_list = [(key + 0.5) / bin_count for key in D.keys()] plt.xticks(range(len(D)), value_list) plt.xlabel(’x’, fontsize=5) plt.ylabel(’counts’, fontsize=5) plt.title(’counting bits’) plt.show()for i in range(90000): x = inverseCDF() # 用bin去映射,否則不好操作 bin = math.floor(x * bin_count) # type(bin): int count_dict[bin] = count_dict.get(bin, 0) + 1draw_pdf(count_dict)
結(jié)果:
代碼實現(xiàn)——方法二,數(shù)值法
數(shù)值模擬cdf的關(guān)鍵是創(chuàng)建lookup table,
table的size越大則結(jié)果越真實(即區(qū)間劃分的個數(shù))
import numpy as npimport mathimport randomimport matplotlib.pyplot as pltimport collectionslookup_table_size = 40CDFlookup_table = np.zeros((lookup_table_size))count_dict = dict()bin_count = 20def inverse_cdf_numerically(y): global lookup_table_size global CDFlookup_table value = 0.0 for i in range(lookup_table_size): x = i * 1.0 / (lookup_table_size - 1) value += pdf2(x) CDFlookup_table[i] = value CDFlookup_table /= value # normalize the cdf if y < CDFlookup_table[0]: t = y / CDFlookup_table[0] return t / lookup_table_size index = -1 for j in range(lookup_table_size): if CDFlookup_table[j] >= y: index = j break # linear interpolation t = (y - CDFlookup_table[index - 1]) / (CDFlookup_table[index] - CDFlookup_table[index - 1]) fractional_index = index + t # 因為index從0開始,所以不是 (index-1)+t return fractional_index / lookup_table_sizedef inverseCDF(): ''' return the x value in PDF ''' uniform_random = random.random() return inverse_cdf_numerically(uniform_random)def pdf2(x): return (x * x * x - 10.0 * x * x + 5.0 * x + 11.0) / (10.417)def draw_pdf(D): global bin_count D = collections.OrderedDict(sorted(D.items())) plt.bar(range(len(D)), list(D.values()), align=’center’) value_list = [(key + 0.5) / bin_count for key in D.keys()] plt.xticks(range(len(D)), value_list) plt.xlabel(’x’, fontsize=5) plt.ylabel(’counts’, fontsize=5) plt.title(’counting bits’) plt.show()for i in range(90000): x = inverseCDF() bin = math.floor(x * bin_count) # type(bin): int count_dict[bin] = count_dict.get(bin, 0) + 1draw_pdf(count_dict)
真實函數(shù)與模擬結(jié)果
擴展:生成伯努利、正太分布
import numpy as npimport matplotlib.pyplot as plt'''reference:https://blog.demofox.org/2017/07/25/counting-bits-the-normal-distribution/'''def plot_bar_x(): # this is for plotting purpose index = np.arange(counting.shape[0]) plt.bar(index, counting) plt.xlabel(’x’, fontsize=5) plt.ylabel(’counts’, fontsize=5) plt.title(’counting bits’) plt.show()# if dice_side=2, is binomial distribution# if dice_side>2 , is multinomial distributiondice_side = 2# if N becomes larger, then multinomial distribution will more like normal distributionN = 100counting = np.zeros(((dice_side - 1) * N + 1))for i in range(30000): sum = 0 for j in range(N): dice_result = np.random.randint(0, dice_side) sum += dice_result counting[sum] += 1# normalizationcounting /= np.sum(counting)plot_bar_x()
以上這篇python 計算概率密度、累計分布、逆函數(shù)的例子就是小編分享給大家的全部內(nèi)容了,希望能給大家一個參考,也希望大家多多支持好吧啦網(wǎng)。
相關(guān)文章:
