dggbeta.m

这个matlab程序用于实现广义高斯分布的参数估计

function [f, g] = dggbeta(beta, absx)
% DGGBETA Derivative of generalized Gaussian pdf with respective to beta
%	[F, G] = DGGBETA(BETA, ABSX)
%	Return partial derivative and second order partial derivative
%	This is an utility function for GGMLE
%
%	ABSX is the distance of the data from the mean
%
% Author: Minh N. Do, Dec. 1999


% Function is undefined for negative beta
if beta <= 0
    f = NaN;
    g = NaN;
    return;
end
    
% Fast computation of \Psi(1/beta) and \Psi(1, 1/beta) based on interpolation
% A table of \Psi(1/beta) for beta = 0.1:0.1:10;
% t1 = mfun('Psi', 10 ./ (1:100));
t1 = [ 2.251753e+00 1.506118e+00 1.046538e+00 7.031566e-01 4.227843e-01 ...
       1.817656e-01 -3.248542e-02 -2.274535e-01 -4.079433e-01 -5.772157e-01 ...
      -7.375804e-01 -8.907294e-01 -1.037936e+00 -1.180181e+00 -1.318234e+00 ...
      -1.452709e+00 -1.584101e+00 -1.712817e+00 -1.839193e+00 -1.963510e+00 ...
      -2.086005e+00 -2.206880e+00 -2.326306e+00 -2.444431e+00 -2.561385e+00 ...
      -2.677277e+00 -2.792208e+00 -2.906261e+00 -3.019514e+00 -3.132034e+00 ...
      -3.243880e+00 -3.355106e+00 -3.465759e+00 -3.575884e+00 -3.685518e+00 ...
      -3.794697e+00 -3.903452e+00 -4.011813e+00 -4.119806e+00 -4.227454e+00 ...
      -4.334779e+00 -4.441802e+00 -4.548542e+00 -4.655014e+00 -4.761235e+00 ...
      -4.867219e+00 -4.972980e+00 -5.078529e+00 -5.183879e+00 -5.289040e+00 ...
      -5.394021e+00 -5.498833e+00 -5.603483e+00 -5.707980e+00 -5.812331e+00 ...
      -5.916543e+00 -6.020622e+00 -6.124576e+00 -6.228409e+00 -6.332128e+00 ...
      -6.435736e+00 -6.539239e+00 -6.642642e+00 -6.745949e+00 -6.849164e+00 ...
      -6.952290e+00 -7.055331e+00 -7.158291e+00 -7.261173e+00 -7.363980e+00 ...
      -7.466715e+00 -7.569380e+00 -7.671979e+00 -7.774513e+00 -7.876985e+00 ...
      -7.979398e+00 -8.081753e+00 -8.184053e+00 -8.286298e+00 -8.388493e+00 ...
      -8.490637e+00 -8.592733e+00 -8.694782e+00 -8.796787e+00 -8.898747e+00 ...
      -9.000666e+00 -9.102543e+00 -9.204381e+00 -9.306181e+00 -9.407943e+00 ...
      -9.509670e+00 -9.611361e+00 -9.713019e+00 -9.814644e+00 -9.916237e+00 ...
      -1.001780e+01 -1.011933e+01 -1.022083e+01 -1.032231e+01 -1.042375e+01];

% A table of \Psi(1, 1/beta) for beta = 0.1:0.1:10;
% t2 = mfun('Psi', 1, 10 ./ (1:100));
t2 = [1.051663e-01 2.213230e-01 3.494237e-01 4.903578e-01 6.449341e-01 ...
      8.138754e-01 9.978204e-01 1.197329e+00 1.412891e+00 1.644934e+00 ...
      1.893830e+00 2.159905e+00 2.443442e+00 2.744693e+00 3.063875e+00 ...
      3.401183e+00 3.756787e+00 4.130837e+00 4.523469e+00 4.934802e+00 ...
      5.364943e+00 5.813987e+00 6.282020e+00 6.769120e+00 7.275357e+00 ...
      7.800792e+00 8.345485e+00 8.909485e+00 9.492842e+00 1.009560e+01 ...
      1.071779e+01 1.135946e+01 1.202063e+01 1.270134e+01 1.340162e+01 ...
      1.412149e+01 1.486098e+01 1.562010e+01 1.639887e+01 1.719733e+01 ...
      1.801548e+01 1.885334e+01 1.971092e+01 2.058824e+01 2.148531e+01 ...
      2.240215e+01 2.333877e+01 2.429517e+01 2.527137e+01 2.626738e+01 ...
      2.728320e+01 2.831884e+01 2.937432e+01 3.044964e+01 3.154480e+01 ...
      3.265981e+01 3.379468e+01 3.494942e+01 3.612403e+01 3.731851e+01 ...
      3.853288e+01 3.976713e+01 4.102126e+01 4.229530e+01 4.358923e+01 ...
      4.490306e+01 4.623679e+01 4.759044e+01 4.896400e+01 5.035747e+01 ...
      5.177086e+01 5.320418e+01 5.465741e+01 5.613057e+01 5.762366e+01 ...
      5.913669e+01 6.066964e+01 6.222253e+01 6.379536e+01 6.538813e+01 ...
      6.700084e+01 6.863350e+01 7.028610e+01 7.195864e+01 7.365114e+01 ...
      7.536358e+01 7.709598e+01 7.884833e+01 8.062063e+01 8.241289e+01 ...
      8.422511e+01 8.605728e+01 8.790941e+01 8.978151e+01 9.167356e+01 ...
      9.358558e+01 9.551756e+01 9.746951e+01 9.944142e+01 1.014333e+02];

if beta < 0.1 | beta > 10
    p1 = mfun('Psi', 1/beta);
    p2 = mfun('Psi', 1, 1/beta);
else
    p1 = interp1(t1, beta*10, '*cubic');
    p2 = interp1(t2, beta*10, '*cubic');
end

n = length(absx);
x1 = absx .^ beta;
x2 = log(absx + eps);

s1 = sum(x1);
s2 = sum(x1 .* x2);
s3 = sum(x1 .* x2 .* x2);
l1 = log(beta * s1 / n + eps);

if s1 == 0
    f = NaN;
    g = NaN;
    return;
end    

f = 1 + (p1 + l1) / beta - s2 / s1;
g = (-p1 - p2/beta + 1 + beta*s2/s1 - l1) / beta^2 - s3/s1 + (s2/s1)^2;