广义高斯分布(GGD)和非对称广义高斯分布(AGGD)的形状参数快速估计

2023-12-27 09:58

本文主要是介绍广义高斯分布(GGD)和非对称广义高斯分布(AGGD)的形状参数快速估计,希望对大家解决编程问题提供一定的参考价值,需要的开发者们随着小编来一起学习吧!

0 引言

广义高斯分布(generalized Gaussian distribution,GGD)和非对称广义高斯分布( asymmetric generalized Gaussian distribution,AGGD)被经常使用与图像/视频信号的统计分析,其形状参数常被用为图像的特征进行分类或回归。如在图像质量评价任务中,Anish Mittal等人提出的BRISQUE模型利用GGD拟合归一化后图像(MSCN)的分布,利用GGD参数 ( α , σ 2 ) (α, σ ^2) (α,σ2)作为一组特征;利用AGGD拟合图像四个方向上的MSCN的内积,每个方向利用AGGD参数 ( η , ν , σ l 2 , σ r 2 ) (η, ν, σ_l ^2, σ_r ^2) (η,ν,σl2,σr2)作为另外四组特征。本文介绍BRISQUE采用的GGD和AGGD形状参数估计方法,并给出MATLAB、Python实现。

f ( x ; α , σ 2 ) = α 2 β Γ ( 1 / α ) exp ⁡ ( − ( ∣ x ∣ β ) α ) f\left(x ; \alpha, \sigma^{2}\right)=\frac{\alpha}{2 \beta \Gamma(1 / \alpha)} \exp \left(-\left(\frac{|x|}{\beta}\right)^{\alpha}\right) f(x;α,σ2)=2βΓ(1/α)αexp((βx)α)
where
β = σ Γ ( 1 / α ) Γ ( 3 / α ) \beta=\sigma \sqrt{\frac{\Gamma(1 / \alpha)}{\Gamma(3 / \alpha)}} β=σΓ(3/α)Γ(1/α)
and Γ ( ⋅ ) \Gamma(\cdot) Γ() is the gamma function:
Γ ( a ) = ∫ 0 ∞ t a − 1 e − t d t a > 0 \Gamma(a)=\int_{0}^{\infty} t^{a-1} e^{-t} d t \quad a>0 Γ(a)=0ta1etdta>0

1 零均值广义高斯分布(GGD)的参数估计

广义高斯分布定义如下:
相同方差不同形状参数的GGD

  GGD共有两个参数:参数 α α α控制着分布的“形状”,即衰减的速度; σ σ σ控制着方差(标准差)。这两个参数的估计方法参考的是论文《Estimation of Shape Parameter for Generalized Gaussian Distributions in Subband Decompositions of Video》的方法。推导过程如下(注意里面的 γ γ γ即为本文的形状参数 α α α):
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
我的数学功底不强,所以我只能对照和论文理解个大概。欢迎大佬们在评论区指正!
方法的核心是下面的比函数:

ρ ( α ) = σ X 2 E 2 [ ∣ X ∣ ] = Γ ( 1 α ) ⋅ Γ ( 3 α ) Γ 2 ( 2 α ) \rho(\alpha)=\frac{\sigma_{X}^{2}}{E^{2}[|X|]}=\frac{\Gamma\left(\frac{1}{\alpha}\right) \cdot \Gamma\left(\frac{3}{\alpha}\right)}{\Gamma^{2}\left(\frac{2}{\alpha}\right)} ρ(α)=E2[X]σX2=Γ2(α2)Γ(α1)Γ(α3)
1.首先计算均值和方差的估计值:
方差估计值即为需要估计的方差 μ ^ X = 1 M ∑ i = 1 11 x i , σ ^ N 2 = 1 M ∑ i = 1 M ( x i − μ ^ N ) 2 \hat{\mu}_{X}=\frac{1}{M} \sum_{i=1}^{11} x_{i}, \quad \hat{\sigma}_{N}^{2}=\frac{1}{M} \sum_{i=1}^{M}\left(x_{i}-\hat{\mu}_{N}\right)^{2} μ^X=M1i=111xi,σ^N2=M1i=1M(xiμ^N)2

需要注意的是,对于零均值广义高斯分布,则 μ ^ X = 0 \hat{\mu}_{X}=0 μ^X=0,方差估计可简化为下式:
σ X 2 ^ = 1 M ∑ i = 1 M x i 2 \widehat{\sigma_{X}^{2}}=\frac{1}{M} \sum_{i=1}^{M} x_{i}^{2} σX2 =M1i=1Mxi2

2.计算绝对值的修正平均值的估计值:
E ^ [ ∣ X ∣ ] = ( 1 / M ) ∑ i = 1 M ∣ x i − μ ^ X ∣ \hat{E}[|X|]=(1 / M) \sum_{i=1}^{M}\left|x_{i}-\hat{\mu}_{X}\right| E^[X]=(1/M)i=1Mxiμ^X

3.设 R ^ \hat{R} R^ ρ ( α ) ρ(α) ρ(α)的估计值,由式(4)得:
R ^ = σ X 2 E 2 [ ∣ X ∣ ] \widehat{\mathrm{R}}=\frac{\sigma_{X}^{2}}{E^{2}[|X|]} R =E2[X]σX2

4.根据 (4)式等式两端相等的原则,利用查找匹配的方法确定 α α α。具体过程如下:
1)给出候选的 α α α:设定 α α α的查找区间,如[0.2:0.001:10]即 α α α在(0.2,10)区间范围每隔0.001取值。
2)将所有 α α α代入式(4)右边计算 ρ ( α ) ρ(α) ρ(α)
3)计算所有 ρ ( α ) ρ(α) ρ(α)与步骤3得出估计值 R ^ \hat{R} R^的距离,选取距离最小对应的 α α α即为所求。

BRISQUE算法作者给出了零均值GGD参数估计的MATLAB代码(注意其中的 γ γ γ即为形状参数(本文的 α α α)):

function [gamparam sigma] = estimateggdparam(vec)
gam                              = 0.2:0.001:10;%给定候选的形状参数γ
r_gam                            = (gamma(1./gam).*gamma(3./gam))./((gamma(2./gam)).^2);%根据式(4)计算比值
sigma_sq                         = mean((vec).^2);%σ^2的零均值估计,非零均值需要计算均值然后按照(3)式
sigma                            = sqrt(sigma_sq);
E                                = mean(abs(vec));%计算步骤2
rho                              = sigma_sq/E^2;%计算步骤3的比值
[min_difference, array_position] = min(abs(rho - r_gam));
gamparam                         = gam(array_position);  

我进行了对应的Python实现:

import numpy as np
from scipy.special import gamma
def estimate_GGD_parameters(vec):gam =np.arange(0.2,10.0,0.001)#产生候选的γr_gam = (gamma(1/gam)*gamma(3/gam))/((gamma(2/gam))**2)#根据候选的γ计算r(γ)sigma_sq=np.mean((vec)**2)#σ^2的零均值估计,非零均值需要计算均值然后按照(3)式sigma=np.sqrt(sigma_sq)E=np.mean(np.abs(vec-u_vec))r=sigma_sq/(E**2)#根据sigma和E计算r(γ)diff=np.abs(r-r_gam)gamma_param=gam[np.argmin(diff, axis=0)]return sigma,gamma_param

2 非对称广义高斯分布(AGGD)的参数估计

AGGD的定义如下:
在这里插入图片描述
其中:
在这里插入图片描述
AGGD共有三个参数, α 、 σ l 、 σ r α、σ_l、σ_r ασlσr α α α控制着分布的“形状”, σ l 、 σ r σ_l、σ_r σlσr分别控制左右两边的扩散程度。AGGD参数的估计方法参考论文《MULTISCALE SKEWED HEAVY TAILED MODEL FOR TEXTURE ANALYSIS》的二阶估计法。推导过程如下:
在这里插入图片描述
在这里插入图片描述
我的数学功底不强,所以我只能对照和论文理解个大概。欢迎大佬们在评论区指正!

方法的核心是以下的比:
r = m 1 2 μ 2 = ( γ 2 + 1 ) 2 ( γ 3 + 1 ) ( γ + 1 ) × ρ ( α ) ( 8 ) r=\frac{m_{1}^{2}}{\mu_{2}}=\frac{\left(\gamma^{2}+1\right)^{2}}{\left(\gamma^{3}+1\right)(\gamma+1)} \times \rho(\alpha) \ \ \ (8) r=μ2m12=(γ3+1)(γ+1)(γ2+1)2×ρ(α)   (8)

1. 首先估计 σ σ σ
σ ^ l = 1 N l − 1 ∑ k = 1 , x k < 0 N l x k 2 ∣ σ ^ r = 1 N r − 1 ∑ k = 1 , x k ≥ 0 N r x k 2 \begin{array}{l} \hat{\sigma}_{l}=\sqrt{\frac{1}{N_{l}-1} \sum_{k=1, x_{k}<0}^{N_{l}} x_{k}^{2} \mid} \\ \hat{\sigma}_{r}=\sqrt{\frac{1}{N_{r}-1} \sum_{k=1, x_{k} \geq 0}^{N_{r}} x_{k}^{2}} \end{array} σ^l=Nl11k=1,xk<0Nlxk2 σ^r=Nr11k=1,xk0Nrxk2

2. 计算 γ γ γ估计量:
γ ^ = β ^ l β ^ r = σ ^ l σ ^ r \hat{\gamma}=\frac{\hat{\beta}_{l}}{\hat{\beta}_{r}}=\frac{\hat{\sigma}_{l}}{\hat{\sigma}_{r}} γ^=β^rβ^l=σ^rσ^l

3.计算 u u u的二阶估计和 m m m的一阶估计, u 、 m u、m um的k阶估计如下:
μ k = E ( X k ) and  m k = E ( ∣ X ∣ k ) \mu_{k}=E\left(X^{k}\right) \text { and } m_{k}=E\left(|X|^{k}\right) μk=E(Xk) and mk=E(Xk)

然后根据(8)式左边估计 r ^ \hat{r} r^
r ^ = m 1 2 μ 2 \hat{r}=\frac{m_{1}^{2}}{\mu_{2}} r^=μ2m12

4. 设 R ^ \hat{R} R^ ρ ( α ) ρ(α) ρ(α)的估计量,则根据式(8):
R ^ = r ^ ( γ ^ 3 + 1 ) ( γ ^ + 1 ) ( γ ^ 2 + 1 ) 2 \hat{R}=\hat{r} \frac{\left(\hat{\gamma}^{3}+1\right)(\hat{\gamma}+1)}{\left(\hat{\gamma}^{2}+1\right)^{2}} R^=r^(γ^2+1)2(γ^3+1)(γ^+1)

5. 和GDD类似,利用查找匹配的方法确定 α α α。具体过程如下:
1)给出候选的 α α α:设定 α α α的查找区间,如[0.2:0.001:10]即 α α α在(0.2,10)区间范围每隔0.001取值。
2)将所有 α α α代入下式计算 ρ ( α ) ρ(α) ρ(α)
ρ ( α ) = Γ ( 2 / α ) 2 Γ ( 1 / α ) Γ ( 3 / α ) \rho(\alpha)=\frac{\Gamma(2 / \alpha)^{2}}{\Gamma(1 / \alpha) \Gamma(3 / \alpha)} ρ(α)=Γ(1/α)Γ(3/α)Γ(2/α)2

3)计算所得的 ρ ( α ) ρ(α) ρ(α)与步骤3估计的 R ^ \hat{R} R^的距离,选取距离最小对应的 α α α即为所求

BRISQUE算法作者给出了零均值AGGD参数估计的MATLAB代码(注意其中的 γ γ γ即为形状参数(本文的 α α α)):

function [alpha leftstd rightstd] = estimateaggdparam(vec)
gam   = 0.2:0.001:10;%给定候选的形状参数γ
r_gam = ((gamma(2./gam)).^2)./(gamma(1./gam).*gamma(3./gam));%步骤5.3计算所有可能的ρ(α)
leftstd            = sqrt(mean((vec(vec<0)).^2));%估计σ_l
rightstd           = sqrt(mean((vec(vec>0)).^2));
gammahat           = leftstd/rightstd;%步骤2
rhat               = (mean(abs(vec)))^2/mean((vec).^2);%步骤3
rhatnorm           = (rhat*(gammahat^3 +1)*(gammahat+1))/((gammahat^2 +1)^2);%步骤4
[min_difference, array_position] = min((r_gam - rhatnorm).^2);
alpha              = gam(array_position);

我进行了对应了Python实现:

def estimate_AGGD_parameters(vec):alpha =np.arange(0.2,10.0,0.001)#产生候选的αr_alpha=((gamma(2/alpha))**2)/(gamma(1/alpha)*gamma(3/alpha))#根据候选的α计算r(α)sigma_l=np.sqrt(np.mean(vec[vec<0]**2))#估计σ_lsigma_r=np.sqrt(np.mean(vec[vec>0]**2))gamma_=sigma_l/sigma_r#步骤2#步骤3u2=np.mean(vec**2)m1=np.mean(np.abs(vec))r_=m1**2/u2R_=r_*(gamma_**3+1)*(gamma_+1)/((gamma_**2+1)**2)#步骤4diff=(R_-r_alpha)**2alpha_param=alpha[np.argmin(diff, axis=0)]return sigma_l,sigma_r,alpha_param

参考文献:
[1]Mittal A , Moorthy A K , Bovik A C . No-Reference Image Quality Assessment in the Spatial Domain[J]. IEEE Transactions on Image Processing A Publication of the IEEE Signal Processing Society, 2012, 21(12):4695.
[2]Sharifi K , Leon-Garcia A . Estimation of shape parameter for generalized Gaussian distribution in subband decomposition of video[J]. IEEE Transactions on Circuits and Systems for Video Technology, 1995, 5(1):52-56.
[3]Lasmar N E , Stitou Y , Berthoumieu Y . Multiscale skewed heavy tailed model for texture analysis[C]// IEEE International Conference on Image Processing. IEEE, 2010.

这篇关于广义高斯分布(GGD)和非对称广义高斯分布(AGGD)的形状参数快速估计的文章就介绍到这儿,希望我们推荐的文章对编程师们有所帮助!



http://www.chinasem.cn/article/542741

相关文章

使用Python快速实现链接转word文档

《使用Python快速实现链接转word文档》这篇文章主要为大家详细介绍了如何使用Python快速实现链接转word文档功能,文中的示例代码讲解详细,感兴趣的小伙伴可以跟随小编一起学习一下... 演示代码展示from newspaper import Articlefrom docx import

Java通过反射获取方法参数名的方式小结

《Java通过反射获取方法参数名的方式小结》这篇文章主要为大家详细介绍了Java如何通过反射获取方法参数名的方式,文中的示例代码讲解详细,感兴趣的小伙伴可以跟随小编一起学习一下... 目录1、前言2、解决方式方式2.1: 添加编译参数配置 -parameters方式2.2: 使用Spring的内部工具类 -

Python调用另一个py文件并传递参数常见的方法及其应用场景

《Python调用另一个py文件并传递参数常见的方法及其应用场景》:本文主要介绍在Python中调用另一个py文件并传递参数的几种常见方法,包括使用import语句、exec函数、subproce... 目录前言1. 使用import语句1.1 基本用法1.2 导入特定函数1.3 处理文件路径2. 使用ex

MySQL中时区参数time_zone解读

《MySQL中时区参数time_zone解读》MySQL时区参数time_zone用于控制系统函数和字段的DEFAULTCURRENT_TIMESTAMP属性,修改时区可能会影响timestamp类型... 目录前言1.时区参数影响2.如何设置3.字段类型选择总结前言mysql 时区参数 time_zon

Python如何使用seleniumwire接管Chrome查看控制台中参数

《Python如何使用seleniumwire接管Chrome查看控制台中参数》文章介绍了如何使用Python的seleniumwire库来接管Chrome浏览器,并通过控制台查看接口参数,本文给大家... 1、cmd打开控制台,启动谷歌并制定端口号,找不到文件的加环境变量chrome.exe --rem

shell脚本快速检查192.168.1网段ip是否在用的方法

《shell脚本快速检查192.168.1网段ip是否在用的方法》该Shell脚本通过并发ping命令检查192.168.1网段中哪些IP地址正在使用,脚本定义了网络段、超时时间和并行扫描数量,并使用... 目录脚本:检查 192.168.1 网段 IP 是否在用脚本说明使用方法示例输出优化建议总结检查 1

Linux中Curl参数详解实践应用

《Linux中Curl参数详解实践应用》在现代网络开发和运维工作中,curl命令是一个不可或缺的工具,它是一个利用URL语法在命令行下工作的文件传输工具,支持多种协议,如HTTP、HTTPS、FTP等... 目录引言一、基础请求参数1. -X 或 --request2. -d 或 --data3. -H 或

Rust中的Option枚举快速入门教程

《Rust中的Option枚举快速入门教程》Rust中的Option枚举用于表示可能不存在的值,提供了多种方法来处理这些值,避免了空指针异常,文章介绍了Option的定义、常见方法、使用场景以及注意事... 目录引言Option介绍Option的常见方法Option使用场景场景一:函数返回可能不存在的值场景

详解Spring Boot接收参数的19种方式

《详解SpringBoot接收参数的19种方式》SpringBoot提供了多种注解来接收不同类型的参数,本文给大家介绍SpringBoot接收参数的19种方式,感兴趣的朋友跟随小编一起看看吧... 目录SpringBoot接受参数相关@PathVariable注解@RequestHeader注解@Reque

Java向kettle8.0传递参数的方式总结

《Java向kettle8.0传递参数的方式总结》介绍了如何在Kettle中传递参数到转换和作业中,包括设置全局properties、使用TransMeta和JobMeta的parameterValu... 目录1.传递参数到转换中2.传递参数到作业中总结1.传递参数到转换中1.1. 通过设置Trans的