0 引言

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

$$f\left(x;\alpha ,{\sigma }^2\right)=\frac{\alpha }{2\beta \Gamma (1/\alpha )}\exp \left(-{\left(\frac{|x|}{\beta }\right)}^{\alpha }\right)$$
where
$$\beta =\sigma \sqrt{\frac{\Gamma (1/\alpha )}{\Gamma (3/\alpha )}}$$
and $\Gamma (\cdot )$ is the gamma function:
$$\Gamma (a)={\int }_0^{\infty }{t}^{a-1}{e}^{-t}dta>0$$

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

广义高斯分布定义如下：

  GGD共有两个参数：参数$\alpha$控制着分布的“形状”，即衰减的速度；$\sigma$控制着方差（标准差）。这两个参数的估计方法参考的是论文《Estimation of Shape Parameter for Generalized Gaussian Distributions in Subband Decompositions of Video》的方法。推导过程如下（注意里面的$\gamma$即为本文的形状参数$\alpha$）：



我的数学功底不强，所以我只能对照和论文理解个大概。欢迎大佬们在评论区指正！
方法的核心是下面的比函数：

$$\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)}$$
1.首先计算均值和方差的估计值：
方差估计值即为需要估计的方差$${\hat{\mu }}_X=\frac{1}{M}\sum _{i=1}^{11}{x}_i,{\hat{\sigma }}_N^2=\frac{1}{M}\sum _{i=1}^M{\left(x_i-{\hat{\mu }}_N\right)}^2$$

需要注意的是，对于零均值广义高斯分布，则${\hat{\mu }}_X=0$,方差估计可简化为下式：
$$\widehat{{\sigma }_X^2}=\frac{1}{M}\sum _{i=1}^M{x}_i^2$$

2.计算绝对值的修正平均值的估计值：
$$\hat{E}[|X|]=(1/M)\sum _{i=1}^M\left|x_i-{\hat{\mu }}_X\right|$$

3.设$\hat{R}$为$\rho (\alpha )$的估计值，由式(4)得：
$$\hat{\mathrm{R}}=\frac{{\sigma }_X^2}{E^2[|X|]}$$

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

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

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共有三个参数，$\alpha 、{\sigma }_l、{\sigma }_r$。$\alpha$控制着分布的“形状”，${\sigma }_l、{\sigma }_r$分别控制左右两边的扩散程度。AGGD参数的估计方法参考论文《MULTISCALE SKEWED HEAVY TAILED MODEL FOR TEXTURE ANALYSIS》的二阶估计法。推导过程如下：


我的数学功底不强，所以我只能对照和论文理解个大概。欢迎大佬们在评论区指正！

方法的核心是以下的比：
$$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)$$

1. 首先估计$\sigma$：
$$\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\ge 0}^{N_r}{x}_k^2}\end{array}$$

2. 计算$\gamma$估计量：
$$\hat{\gamma }=\frac{{\hat{\beta }}_l}{{\hat{\beta }}_r}=\frac{{\hat{\sigma }}_l}{{\hat{\sigma }}_r}$$

3.计算$u$的二阶估计和$m$的一阶估计，$u、m$的k阶估计如下：
$${\mu }_k=E\left(X^k\right)\text{ and }{m}_k=E\left(|X{|}^k\right)$$

然后根据(8)式左边估计$\hat{r}$：
$$\hat{r}=\frac{m_1^2}{{\mu }_2}$$

4. 设$\hat{R}$为$\rho (\alpha )$的估计量，则根据式(8)：
$$\hat{R}=\hat{r}\frac{\left({\hat{\gamma }}^3+1\right)(\hat{\gamma }+1)}{{\left({\hat{\gamma }}^2+1\right)}^2}$$

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

3）计算所得的$\rho (\alpha )$与步骤3估计的$\hat{R}$的距离，选取距离最小对应的$\alpha$即为所求

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

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))#估计σ_l
    sigma_r=np.sqrt(np.mean(vec[vec>0]**2))
    gamma_=sigma_l/sigma_r#步骤2
    #步骤3
    u2=np.mean(vec**2)
    m1=np.mean(np.abs(vec))
    r_=m1**2/u2
    R_=r_*(gamma_**3+1)*(gamma_+1)/((gamma_**2+1)**2)#步骤4
    diff=(R_-r_alpha)**2
    alpha_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.