python实现四位一并法_(番外篇1) 贝叶斯超参数的确定方法(附自适应搜索法、牛顿迭代法和二分查找法的python实现)...

最新推荐文章于 2024-06-28 01:30:54 发布

weixin_39543478

最新推荐文章于 2024-06-28 01:30:54 发布

阅读量284

点赞数

文章标签： python实现四位一并法

(番外篇1) 贝叶斯超参数的确定方法(附自适应搜索法、牛顿迭代法和二分查找法的python实现)

背景说明

在贝叶斯统计中，总体分布的参数θ\thetaθ不再被视为一个未知的常数处理，而是认为参数θ\thetaθ也有一个分布，事实上，由于θ\thetaθ是未知的不确定的东西，因此用概率分布的语言来描述是很有优势的。

参数θ\thetaθ一般需要先综合已有信息，确定好先验分布π(θ)\pi(\theta)π(θ)，然后利用新进的样本信息进行更新，从而得到后验分布π(θ∣x)\pi(\theta|x)π(θ∣x)，这个过程中如果π(θ)\pi(\theta)π(θ)与π(θ∣x)\pi(\theta|x)π(θ∣x)仍然是同一种分布，则称π(θ)\pi(\theta)π(θ)是总体分布参数θ\thetaθ的共轭先验分布，共轭先验分布具有很好的性质，体现在计算简单和具备可解释性上。

常用的共轭先验分布有

总体分布

参数

共轭先验分布

二项分布

发生概率

贝塔分布

正态分布

均值

正态分布

方差

倒伽马分布

泊松分布

均值

伽马分布

指数分布

均值倒数

伽马分布我们知道，二项分布在现实中是广泛存在的，因此当我们选择其先验分布为贝塔分布Be(α,β)Be(\alpha,\beta)Be(α,β)时(共轭)，我们还需要给其两个参数赋予初始值，这就是超参数确定问题。

确定超参数的几个方法

利用先验矩

利用先验分位数

利用先验矩和先验分位数

其他方法

一个实例：二项分布参数p的超参数确定

思路分析

根据上文所讲的，二项分布的共轭先验分布为Be(α,β)Be(\alpha,\beta)Be(α,β)，我们将采用3. 利用先验分位数来确定。

我们根据先验信息，知道了θ\thetaθ的两个分位数，比如为上下四分位数θL\theta_{L}θL和θU\theta_{U}θU，因此可以得到两个方程式：

∫0θLΓ(α+β)Γ(α)Γ(β)θα−1(1−θ)β−1dθ=0.25∫0θuΓ(α+β)Γ(α)Γ(β)θα−1(1−θ)β−1dθ=0.75

\begin{aligned}

&\int_{0}^{\theta_L}\frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}\theta^{\alpha-1}(1-\theta)^{\beta-1}d\theta = 0.25 \\

&\int_{0}^{\theta_{u}}\frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}\theta^{\alpha-1}(1-\theta)^{\beta-1}d\theta = 0.75

\end{aligned}

∫0θLΓ(α)Γ(β)Γ(α+β)θα−1(1−θ)β−1dθ=0.25∫0θuΓ(α)Γ(β)Γ(α+β)θα−1(1−θ)β−1dθ=0.75

由于我们关注的变量是α\alphaα和β\betaβ，因此上式可以抽象为

f1(α,β)=0f2(α,β)=0

\begin{aligned}

f_1(\alpha,\beta) = 0 \\

f_2(\alpha,\beta) = 0

\end{aligned}

f1(α,β)=0f2(α,β)=0这是一个非线性方程组的零点问题。理论上，利用牛顿迭代的二阶形式，在一定条件下，是可以求解此类问题的。但是这涉及到α\alphaα和β\betaβ的求导，是比较复杂的。

可以考虑从另一个角度来求解这个问题。首先这是一个超参数的初始值确定问题，因此结果并不要求十分的精准，我们可以在损失一定精度的前提下得到α\alphaα和β\betaβ的近似值。

因此我们可以在一个α−β平面\alpha-\beta平面α−β平面内对不同的α\alphaα和β\betaβ进行多轮计算，通过计算的精确程度来控制计算的时间。比如我们大概可以确定α\alphaα和β\betaβ的范围分别在[1,100][1,100][1,100]和[1,100][1,100][1,100]之间，那么就可以某个步长遍历这个α−β平面\alpha-\beta平面α−β平面，在遍历的过程中找到使得∥f1∥+∥f2∥\|f_1\|+\|f_2\|∥f1∥+∥f2∥最小的α\alphaα和β\betaβ，从而将得到的α\alphaα和β\betaβ视为非线性方程的近似解。由于f1f_1f1和f2f_2f2是关于α\alphaα和β\betaβ的连续函数，因此可以保证结果的可靠性。

python代码实现

import numpy as np

from scipy.integrate import quad

def gamma(alpha):

return quad(lambda x:x**(alpha-1)*np.exp(-x),0,np.inf)[0]

def Fbeta_q(alpha,beta,theta,q):

'''

本函数返回beta分布函数theta对应的概率值与理想的概率值q之间的差值绝对值

'''

k = gamma(alpha+beta)/(gamma(alpha)*gamma(beta))

integValue = quad(lambda x:x**(alpha-1)*(1-x)**(beta-1),0,theta)[0]

return abs(k*integValue-q)

def betaSolve(theta1,q1,theta2,q2,lowAlpha = 0.1,highAlpha = 20,lowBeta = 0.1,highBeta = 20,step = 0.2):

'''

本函数接受两对分位数，返回计算近似到的beta分布的两个参数值及相应的误差

参数的意义如下：

theta1 q1 theta2 q2 为两对分位数

lowAlpha highAlpha lowBeta highBeta 为两个参数的搜索范围，默认均为(0,20)

step 为搜索步长，默认为 0.2 (超参数的初始确定并不需要十分精确)

'''

AlphaPoints = np.arange(lowAlpha+0.1,highAlpha,step)

BetaPoints = np.arange(lowBeta+0.1,highBeta,step)

xAlpha,yBeta = np.meshgrid(AlphaPoints,BetaPoints)

minError = np.inf

minAlpha = np.NaN

minBeta = np.NaN

for i in range(len(AlphaPoints)):

for j in range(len(BetaPoints)):

currentError = Fbeta_q(alpha = xAlpha[i,j], beta = yBeta[i,j],theta = theta1, q = q1)\

+ Fbeta_q(alpha = xAlpha[i,j], beta = yBeta[i,j],theta = theta2, q = q2)

#print(xAlpha[i,j],yBeta[i,j],currentError,minError,(currentError < minError))

if(currentError < minError):

minAlpha = xAlpha[i,j]

minBeta = yBeta[i,j]

minError = currentError

return minAlpha,minBeta,minError

数据测试

我们已知的几组数据(上下四分位数)如下

Be(α,β)Be(\alpha,\beta)Be(α,β)

θL\theta_LθL

θU\theta_UθU

(1,1)

0.25

0.75

(2,8)

0.107

0.274

(6,11)

0.224

0.368

(9,11)

0.372

0.526

测试代码如下所示：

import betaQuantile

import betaParameters

if __name__ == '__main__':

alpha,beta,error = betaParameters.betaSolve(theta1 = 0.25,q1 = 0.25,theta2 = 0.75,q2 = 0.75,step = 0.1,lowAlpha = 0.1,highAlpha = 5,lowBeta = 0.1,highBeta = 5)

print('theta1 = 0.250','theta1 = 0.750','alpha = ',alpha,'beta = ',beta,'error = ',error)

alpha,beta,error = betaParameters.betaSolve(theta1 = 0.107,q1 = 0.25,theta2 = 0.274,q2 = 0.75,step = 0.1,lowAlpha = 0.1,highAlpha = 10,lowBeta = 0.1,highBeta = 10)

print('theta1 = 0.107','theta1 = 0.274','alpha = ',alpha,'beta = ',beta,'error = ',error)

alpha,beta,error = betaParameters.betaSolve(theta1 = 0.224,q1 = 0.25,theta2 = 0.368,q2 = 0.75,step = 0.2,lowAlpha = 0.1,highAlpha = 20,lowBeta = 0.1,highBeta = 20)

print('theta1 = 0.224','theta1 = 0.368','alpha = ',alpha,'beta = ',beta,'error = ',error)

alpha,beta,error = betaParameters.betaSolve(theta1 = 0.372,q1 = 0.25,theta2 = 0.526,q2 = 0.75,step = 0.2,lowAlpha = 0.1,highAlpha = 20,lowBeta = 0.1,highBeta = 20)

print('theta1 = 0.372','theta1 = 0.526','alpha = ',alpha,'beta = ',beta,'error = ',error)

测试结果如下所示：

真实(α,β)真实(\alpha,\beta)真实(α,β)

近似(α,β)近似(\alpha,\beta)近似(α,β)

(1,1)

(1.0,1.0)

(2,8)

(2.0,8.0)

(6,11)

(5.4,12.6)

(9,11)

(8.8,10.8)

扩展：贝塔分布的上下四分位查询表格

对于任意已知的α\alphaα和β\betaβ，事实上我们都可以计算它的分位数。因此，我们可以对整数组合的α\alphaα和β\betaβ进行计算上下四分位数，从而可以通过查表的方式快速获得其整数近似值。

我们的问题实际上是求解一个非线性方程的零点

∫0θΓ(α+β)Γ(α)Γ(β)θα−1(1−θ)β−1dθ=p

\begin{aligned}

&\int_{0}^{\theta}\frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}\theta^{\alpha-1}(1-\theta)^{\beta-1}d\theta = p

\end{aligned}

∫0θΓ(α)Γ(β)Γ(α+β)θα−1(1−θ)β−1dθ=p

由于我们关注的变量是θ\thetaθ，因此上式可以抽象为

f(θ)=FBe(θ)−p=0f(\theta) = F_{Be}(\theta) - p = 0f(θ)=FBe(θ)−p=0该方程必然有唯一实解，我们有两种方式可以求解该方程：

1. 间接法(迭代法、牛顿法)

2. 直接法(二分法)

利用牛顿迭代的话，需要对f(θ)f(\theta)f(θ)求导，这是简单的，其为

f′(θ)=∂FBe(θ)∂θ=fBe(θ)=Γ(α+β)Γ(α)Γ(β)θα−1(1−θ)β−1f^{'}(\theta) = \frac{\partial{F_{Be}(\theta)}} {\partial{\theta}} = f_{Be}(\theta) = \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}\theta^{\alpha-1}(1-\theta)^{\beta-1}f′(θ)=∂θ∂FBe(θ)=fBe(θ)=Γ(α)Γ(β)Γ(α+β)θα−1(1−θ)β−1

因此迭代式为

θ:=θ−FBe(θ)fBe′(θ)\theta: = \theta - \frac{F_{Be}(\theta)}{f_{Be}^{'}(\theta)}θ:=θ−fBe′(θ)FBe(θ)

事实上，在这里用牛顿法是不太合适的，根据牛顿法的收敛条件，这里并不能保证其收敛到[0,1][0,1][0,1]的零点(后面我们将看到有时会收敛到负值)，因为FBe(θ)F_{Be}(\theta)FBe(θ)其实在[0,1][0,1][0,1]区间上并不是凸函数或凹函数。更合适的方法是利用二分法直接获得零点，因为我们已经知道其在[0,1][0,1][0,1]区间上存在唯一的零点。二分法思想较简单，这里略去不写。

求上下四分位点的python实现

import numpy as np

from scipy.integrate import quad

def gamma(alpha):

'''

本函数定义了伽马函数

'''

return quad(lambda x:x**(alpha-1)*np.exp(-x),0,np.inf)[0]

def fbeta(alpha,beta,theta):

'''

本函数返回beta分布密度函数theta对应的概率密度值

'''

k = gamma(alpha+beta)/(gamma(alpha)*gamma(beta))

value = theta**(alpha-1)*(1-theta)**(beta-1)

return k*value

def Fbeta(alpha,beta,theta):

'''

本函数返回beta分布函数 P{x <= theta} 的概率值

'''

k = gamma(alpha+beta)/(gamma(alpha)*gamma(beta))

integValue = quad(lambda x:x**(alpha-1)*(1-x)**(beta-1),0,theta)[0]

return k*integValue

def betaQuantileSolve(alpha,beta,p,theta1 = 0.5,epsilon = 0.0001,method = 'Bisection'):

'''

本函数对根据特定参数的贝塔分布，计算其p分位数，返回其值。

参数意义如下：

alpha beta 贝塔分布的参数

p 要计算的分位数q，使得 P(x <= q) = p

theta1 = 0.5 epsilon = 0.0001 初始值事后判断的收敛阀值

method = 'Bisection' 所用的方法，method有两个选择，一是'Bisection'，而是'Newton'

'''

if(method == 'Newton'):

theta0 = np.inf

while(abs(theta1-theta0) > epsilon):

theta0 = theta1

theta1 = theta0 - (Fbeta(alpha,beta,theta0)-p)/fbeta(alpha,beta,theta0)

return theta1

elif(method == 'Bisection'):

high = 1

low = 0

middle = (high+low)/2

while(high - low > epsilon):

direction = Fbeta(alpha,beta,middle)-p

if(direction > 0):

high = middle

elif(direction < 0):

low = middle

else:

return middle

middle = (high+low)/2

return middle

数据测试

Be(α,β)Be(\alpha,\beta)Be(α,β)

θL\theta_LθL

θU\theta_UθU

(1,1)

0.25

0.75

(2,8)

0.107

0.274

(6,11)

0.224

0.368

(9,11)

0.372

0.526

测试代码如下所示

import betaQuantile

import betaParameters

if __name__ == '__main__':

#测试计算上下四分位数的程序(betaQuantile)

print('测试计算上下四分位数的程序(betaQuantile)')

theta_L = betaQuantile.betaQuantileSolve(alpha = 1,beta = 1,p = 0.25,method = 'Newton')

theta_U = betaQuantile.betaQuantileSolve(alpha = 1,beta = 1,p = 0.75,method = 'Newton')

print('alpha = 1','beta = 1','theta_L = ',theta_L,'theta_U = ',theta_U)

theta_L = betaQuantile.betaQuantileSolve(alpha = 2,beta = 8,p = 0.25,method = 'Newton')

theta_U = betaQuantile.betaQuantileSolve(alpha = 2,beta = 8,p = 0.75,method = 'Newton')

print('alpha = 2','beta = 8','theta_L = ',theta_L,'theta_U = ',theta_U)

theta_L = betaQuantile.betaQuantileSolve(alpha = 6,beta = 11,p = 0.25,method = 'Newton')

theta_U = betaQuantile.betaQuantileSolve(alpha = 6,beta = 11,p = 0.75,method = 'Newton')

print('alpha = 6','beta = 11','theta_L = ',theta_L,'theta_U = ',theta_U)

theta_L = betaQuantile.betaQuantileSolve(alpha = 9,beta = 11,p = 0.25,method = 'Newton')

theta_U = betaQuantile.betaQuantileSolve(alpha = 9,beta = 11,p = 0.75,method = 'Newton')

print('alpha = 9','beta = 11','theta_L = ',theta_L,'theta_U = ',theta_U)

测试结果：

从结果中可以看到，牛顿迭代无法保证其收敛到[0,1][0,1][0,1]区间上的零点，有些出现了负值。因此，可以换二分法来计算，测试代码如下：

import betaQuantile

import betaParameters

if __name__ == '__main__':

#测试计算上下四分位数的程序(betaQuantile)(Bisection)

print('测试计算上下四分位数的程序(betaQuantile)(Bisection)')

theta_L = betaQuantile.betaQuantileSolve(alpha = 1,beta = 1,p = 0.25,method = 'Bisection')

theta_U = betaQuantile.betaQuantileSolve(alpha = 1,beta = 1,p = 0.75,method = 'Bisection')

print('alpha = 1','beta = 1','theta_L = ',theta_L,'theta_U = ',theta_U)

theta_L = betaQuantile.betaQuantileSolve(alpha = 2,beta = 8,p = 0.25,method = 'Bisection')

theta_U = betaQuantile.betaQuantileSolve(alpha = 2,beta = 8,p = 0.75,method = 'Bisection')

print('alpha = 2','beta = 8','theta_L = ',theta_L,'theta_U = ',theta_U)

theta_L = betaQuantile.betaQuantileSolve(alpha = 6,beta = 11,p = 0.25,method = 'Bisection')

theta_U = betaQuantile.betaQuantileSolve(alpha = 6,beta = 11,p = 0.75,method = 'Bisection')

print('alpha = 6','beta = 11','theta_L = ',theta_L,'theta_U = ',theta_U)

theta_L = betaQuantile.betaQuantileSolve(alpha = 9,beta = 11,p = 0.25,method = 'Bisection')

theta_U = betaQuantile.betaQuantileSolve(alpha = 9,beta = 11,p = 0.75,method = 'Bisection')

print('alpha = 9','beta = 11','theta_L = ',theta_L,'theta_U = ',theta_U)

测试结果

可以看到，二分法的结果即是正确的。

weixin_39543478

关注

0
点赞
踩
0

收藏

觉得还不错? 一键收藏
0
评论
复制链接

分享到 QQ

分享到新浪微博

扫一扫