python 計算概率密度、累計分佈、逆函式的例子
計算概率分佈的相關引數時，一般使用 scipy 包，常用的函式包括以下幾個：
pdf：連續隨機分佈的概率密度函式
pmf：離散隨機分佈的概率密度函式
cdf：累計分佈函式
百分位函式（累計分佈函式的逆函式）
生存函式的逆函式（1 - cdf 的逆函式）
函式裡面不僅能跟一個數據，還能跟一個數組。下面用正態分佈舉例說明：
>>> import scipy.stats as st >>> st.norm.cdf(0) # 標準正態分佈在 0 處的累計分佈概率值 0.5 >>> st.norm.cdf([-1,1])# 標準正態分佈分別在 -1， 0， 1 處的累計分佈概率值 array([0.15865525,0.5,0.84134475]) >>> st.norm.pdf(0) # 標準正態分佈在 0 處的概率密度值 0.3989422804014327 >>> st.norm.ppf(0.975)# 標準正態分佈在 0.975 處的逆函式值 1.959963984540054 >>> st.norm.lsf(0.975)# 標準正態分佈在 0.025 處的生存函式的逆函式值 1.959963984540054
對於非標準正態分佈，通過更改引數 loc 與 scale 來改變均值與標準差：
>>> st.norm.cdf(0,loc=2,scale=1) # 均值為 2，標準差為 1 的正態分佈在 0 處的累計分佈概率值 0.022750131948179195
對於其他隨機分佈，可能更改的引數不一樣，具體需要查官方文件。下面我們舉一些常用分佈的例子：
>>> st.binom.pmf(4,n=100,p=0.05) # 引數值 n=100,p=0.05 的二項分佈在 4 處的概率密度值 0.17814264156968956 >>> st.geom.pmf(4,p=0.05) # 引數值 p=0.05 的幾何分佈在 4 處的概率密度值 0.04286875 >>> st.poisson.pmf(2,mu=3) # 引數值 mu=3 的泊松分佈在 2 處的概率密度值 0.22404180765538775 >>> st.chi2.ppf(0.95,df=10) # 自由度為 10 的卡方分佈在 0.95 處的逆函式值 18.307038053275146 >>> st.t.ppf(0.975,df=10) # 自由度為 10 的 t 分佈在 0.975 處的逆函式值 2.2281388519649385 >>> st.f.ppf(0.95,dfn=2,dfd=12) # 自由度為 2,12 的 F 分佈在 0.95 處的逆函式值 3.8852938346523933
補充拓展：給定概率密度,生成隨機數 python實現
實現的方法可以不止一種：
rejection sampling
invert the cdf
Metropolis Algorithm (MCMC)
本篇介紹根據累積概率分佈函式的逆函式(2：invert the CDF)生成的方法。
自己的理解不一定正確，有錯誤望指正。
目標：
已知 y=pdf(x），現想由給定的pdf,生成對應分佈的x
PDF是概率分佈函式，對其積分或者求和可以得到CDF（累積概率分佈函式），PDF積分或求和的結果始終為1
步驟（具體解釋後面會說）：
1、根據pdf得到cdf
2、由cdf得到inverse of the cdf
3、對於給定的均勻分佈[0,1),帶入inverse cdf，得到的結果即是我們需要的x
求cdf逆函式的具體方法：
對於上面的第二步，可以分成兩類：
1、當CDF的逆函式好求時，直接根據公式求取，
2、反之當CDF的逆函式不好求時，用數值模擬方法
自己的理解：為什麼需要根據cdf的逆去獲得x？
原因一：
因為cdf是單調函式因此一定存在逆函式（cdf是s型函式，而pdf則不一定，例如正態分佈，不單調，對於給定的y，可能存在兩個對應的x，就不可逆）
原因二：
這僅是我自己的直觀理解，根據下圖所示（左上為pdf，右上為cdf）
由步驟3可知，我們首先生成[0，1)的均勻隨機數，此隨機數作為cdf的y，去對映到cdf的x（若用cdf的逆函式表示則是由x對映到y），可以參考上圖的右上，既然cdf的y是均勻隨機的，那麼對於cdf中同樣範圍的x，斜率大的部分將會有更大的機會被對映，因為對應的y範圍更大（而y是隨即均勻分佈的），那麼，cdf的斜率也就等同於pdf的值，這正好符合若x的pdf較大，那麼有更大的概率出現（即重複很多次後，該x會出現的次數最多）
程式碼實現——方法一，公式法
import numpy as np
import math
import random
import matplotlib.pyplot as plt
import collections
count_dict = dict()
bin_count = 20
def 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,其逆函式很好求，因此直接用公式法
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) + 1
draw_pdf(count_dict)
結果：
程式碼實現——方法二，數值法
數值模擬cdf的關鍵是建立lookup table，
table的size越大則結果越真實（即區間劃分的個數）
import numpy as np
import math
import random
import matplotlib.pyplot as plt
import collections
lookup_table_size = 40
CDFlookup_table = np.zeros((lookup_table_size))
count_dict = dict()
bin_count = 20
def 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_size
def 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)),align='center')
value_list = [(key + 0.5) / bin_count for key in D.keys()]
plt.xticks(range(len(D)),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) + 1
draw_pdf(count_dict)
真實函式與模擬結果
擴充套件：生成伯努利、正太分佈
import numpy as np
import 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.title('counting bits')
plt.show()
# if dice_side=2,is binomial distribution
# if dice_side>2,is multinomial distribution
dice_side = 2
# if N becomes larger,then multinomial distribution will more like normal distribution
N = 100
counting = 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
# normalization
counting /= np.sum(counting)
plot_bar_x()
以上這篇python 計算概率密度、累計分佈、逆函式的例子就是小編分享給大家的全部內容了，希望能給大家一個參考，也希望大家多多支援我們。