幂律分布拟合曲线

plot (on log-log axes) the empirical distribution along with the fitted power-law distribution

function h=plplot（x xmin alpha）
% reshape input vector
x = reshape（xnumel（x）1）;
% initialize storage for output handles
h = zeros（21）;

% select method （discrete or continuous） for plotting
if     isempty（setdiff（xfloor（x））） f_dattype = ‘INTS‘;
elseif isreal（x）    f_dattype = ‘REAL‘;
else                 f_dattype = ‘UNKN‘;
end;
if strcmp（f_dattype‘INTS‘） && min（x） > 50
    f_dattype = ‘REAL‘;
end;

% estimate xmin and alpha accordingly
switch f_dattype
    
    case ‘REAL‘
        n = length（x）;
        c = [sort（x） （n:-1:1）‘./n];
        q = sort（x（x>=xmin））;
        cf = [q （q./xmin）.^（1-alpha）];
        cf（:2） = cf（:2） .* c（find（c（:1）>=xmin1‘first‘）2）;

        figure;
        h（1） = loglog（c（:1）c（:2）‘bo‘‘MarkerSize‘8‘MarkerFaceColor‘[1 1 1]）; hold on;
        h（2） = loglog（cf（:1）cf（:2）‘k--‘‘LineWidth‘2）; hold off;
        xr  = [10.^floor（log10（min（x））） 10.^ceil（log10（max（x）））];
        xrt = （round（log10（xr（1）））:2:round（log10（xr（2））））;
        if length（xrt）<4 xrt = （round（log10（xr（1）））:1:round（log10（xr（2））））; end;
        yr  = [10.^floor（log10（1/n）） 1];
        yrt = （round（log10（yr（1）））:2:round（log10（yr（2））））;
        if length（yrt）<4 yrt = （round（log10（yr（1）））:1:round（log10（yr（2））））; end;
        set（gca‘XLim‘xr‘XTick‘10.^xrt）;
        set（gca‘YLim‘yr‘YTick‘10.