0

想要用对数轴在 PDF 上绘制最佳拟合线。x 轴为 10 0到 10 12吉瓦。y轴是10 -2到10 16，是概率分布。到目前为止，我已经尝试过：

pPoly = polyfit(bin,prob_dens_mean,1);
linePointsX = logspace(min(log10(bin)), max(log10(bin)), 200);
linePointsY = pPoly(1)*linePointsX+pPoly(2);

h = figure; loglog(bin,prob_dens_mean,'b*','MarkerSize',10) hold on; plot(linePointsX,linePointsY,'-r');

但这将最佳拟合线绘制为一条水平线，在分布的远端拖尾，而实际的发行版具有负斜率。

下面的代码：

% Purpose:  To characterize the clustered latent heat distribution
% Input:
% W_U_T_Jan = latent heat for each cluster, AM (J)

% Output:
% 1PDF chart per month

clc;
clear all;

load Jan2005_basin_variables.mat

% Initialize

j_len = length(W_U_T_Jan);
prob_dens_all = zeros(j_len,20);
ii = 1 : j_len;
count(1:20) = 0;
bin(1:20) = 0;

for i = 1 : 20
    bin(i) = 10^(0.6*i);
end

%histo = histc(log10(W_U_T_Jan),10:0.5:20.0);
% Bin the Watts

for i = 1 : j_len
    if((log10(W_U_T_Jan(i)) >= 1.25) && (log10(W_U_T_Jan(i)) < 1.75))
        count(1) = count(1) + 1;
    end
    if((log10(W_U_T_Jan(i)) >= 1.75) && (log10(W_U_T_Jan(i)) < 2.25))
        count(2) = count(2) + 1;
    end
    if((log10(W_U_T_Jan(i)) >= 2.25) && (log10(W_U_T_Jan(i)) < 2.75))
        count(3) = count(3) + 1;
    end
    if((log10(W_U_T_Jan(i)) >= 2.75) && (log10(W_U_T_Jan(i)) < 3.25))
        count(4) = count(4) + 1;
    end
    if((log10(W_U_T_Jan(i)) >= 3.25) && (log10(W_U_T_Jan(i)) < 3.75))
        count(5) = count(5) + 1;
    end
    if((log10(W_U_T_Jan(i)) >= 3.75) && (log10(W_U_T_Jan(i)) < 4.25))
        count(6) = count(6) + 1;
    end
    if((log10(W_U_T_Jan(i)) >= 4.25) && (log10(W_U_T_Jan(i)) < 4.75))
        count(7) = count(7) + 1;
    end
    if((log10(W_U_T_Jan(i)) >= 5.25) && (log10(W_U_T_Jan(i)) < 5.75))
        count(8) = count(8) + 1;
    end
    if((log10(W_U_T_Jan(i)) >= 5.75) && (log10(W_U_T_Jan(i)) < 6.25))
        count(9) = count(9) + 1;
    end
    if((log10(W_U_T_Jan(i)) >= 6.25) && (log10(W_U_T_Jan(i)) < 6.75))
        count(10) = count(10) + 1;
    end
    if((log10(W_U_T_Jan(i)) >= 6.75) && (log10(W_U_T_Jan(i)) < 7.25))
        count(11) = count(11) + 1;
    end
    if((log10(W_U_T_Jan(i)) >= 7.25) && (log10(W_U_T_Jan(i)) < 7.75))
        count(12) = count(12) + 1;
    end
    if((log10(W_U_T_Jan(i)) >= 7.75) && (log10(W_U_T_Jan(i)) < 8.25))
        count(13) = count(13) + 1;
    end
    if((log10(W_U_T_Jan(i)) >= 8.25) && (log10(W_U_T_Jan(i)) < 8.75))
        count(14) = count(14) + 1;
    end
    if((log10(W_U_T_Jan(i)) >= 8.75) && (log10(W_U_T_Jan(i)) < 9.25))
        count(15) = count(15) + 1;
    end
    if((log10(W_U_T_Jan(i)) >= 9.25) && (log10(W_U_T_Jan(i)) < 9.75))
        count(16) = count(16) + 1;
    end
    if((log10(W_U_T_Jan(i)) >= 9.75) && (log10(W_U_T_Jan(i)) < 10.25))
        count(17) = count(17) + 1;
    end
    if((log10(W_U_T_Jan(i)) >= 10.25) && (log10(W_U_T_Jan(i)) < 10.75))
        count(18) = count(18) + 1;
    end
    if((log10(W_U_T_Jan(i)) >= 10.75) && (log10(W_U_T_Jan(i)) < 11.25))
        count(19) = count(19) + 1;
    end
    if((log10(W_U_T_Jan(i)) >= 11.25) && (log10(W_U_T_Jan(i)) < 11.75))
        count(20) = count(20) + 1;
    end
end



for i=1:20
    prob(i) = count(i)/sum(count);
    prob_dens(i) = prob(i)/bin(i);
end

% Check
sum(prob_dens.*bin);
prob_dens_all(i,:) = prob_dens(:);

%end

prob_dens_mean = zeros(1,20);

for i = 1 : 20
  prob_dens_mean(1,i) = mean(prob_dens_all(:,i));
end

% Plot

%best_fit = polyfit(log(bin),log(prob_dens_mean),1);



% pPoly = polyfit(bin,prob_dens_mean,1)
% linePointsX = [min(log10(bin)) max(log10(bin))]
% linePointsY = polyval(pPoly,[min(log10(bin)),max(log10(bin))])

pPoly = polyfit(bin,prob_dens_mean,1);
linePointsX = logspace(min(log10(bin)), max(log10(bin)), 200);
linePointsY = pPoly(1)*linePointsX+pPoly(2);

h = figure;
loglog(bin,prob_dens_mean,'b*','MarkerSize',10)
hold on;
plot(linePointsX,linePointsY,'-r');

%loglog(best_fit,'ro')
t = title('Event Power Distribution, Tropics, Jan 2005');
set(t, 'FontWeight', 'bold', 'FontSize', 12)
set(gca, 'FontWeight', 'bold', 'FontSize', 12)
set(gca,'XLim',[10^0 10^12],'YLim',[10^(-16) 10^(-2)]);
xlabel('Event Power (GW)');
ylabel('Probability Density');
print -dpng Tropics_Wattage_PDF_Jan2005.png
4

1 回答 1

1

I think you are having problems because of mixing linear space and log space. Try fitting and plotting in log space like this:

pPoly = polyfit(log10(bin),log10(prob_dens_mean),1);
linePointsX = logspace(min(bin), max(bin), 200);
linePointsY = polyval(pPoly,log10(linePointsX));

figure
loglog(bin,prob_dens_mean,'b*','MarkerSize',10)
hold on
loglog(linePointsX,linePointsY,'-r');

Also I recommend using polyval instead of computing linePointsY yourself.

于 2013-08-29T19:08:21.570 回答