Best fit line on PDF plot with logarithmic axes

Question

Kevin 2013-8-29

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/85950-best-fit-line-on-pdf-plot-with-logarithmic-axes

Want to plot a best fit line on a PDF with logarithmic axes. The x-axis is 10^0 to 10^12 Gigawatts. The y-axis is 10^-2 to 10^16, and is the probability distribution. So far I have tried:

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');

But this plots the best fit line as a horizontal line, tailing off at the far end of the distribution, while the actual distro has a negative slope.

Thanks for your help!

Code below:

% 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
%     if((log10(W_U_T_Jan(i)) >= 17.3) && (log10(W_U_T_Jan(i)) < 17.6))
%         count(21) = count(21) + 1;
%     end
%     if((log10(W_U_T_Jan(i)) >= 17.6) && (log10(W_U_T_Jan(i)) < 17.9))
%         count(22) = count(22) + 1;
%     end
%     if((log10(W_U_T_Jan(i)) >= 17.9) && (log10(W_U_T_Jan(i)) < 18.2))
%         count(23) = count(23) + 1;
%     end
%     if((log10(W_U_T_Jan(i)) >= 18.2) && (log10(W_U_T_Jan(i)) < 18.5))
%         count(24) = count(24) + 1;
%     end
%     if((log10(W_U_T_Jan(i)) >= 18.5) && (log10(W_U_T_Jan(i)) < 18.8))
%         count(25) = count(25) + 1;
%     end
%     if((log10(W_U_T_Jan(i)) >= 18.8) && (log10(W_U_T_Jan(i)) < 19.1))
%         count(26) = count(26) + 1;
%     end
%     if((log10(W_U_T_Jan(i)) >= 19.1) && (log10(W_U_T_Jan(i)) < 19.4))
%         count(27) = count(27) + 1;
%     end
%     if((log10(W_U_T_Jan(i)) >= 19.4) && (log10(W_U_T_Jan(i)) < 19.7))
%         count(28) = count(28) + 1;
%     end
%     if((log10(W_U_T_Jan(i)) >= 19.7) && (log10(W_U_T_Jan(i)) < 20.0))
%         count(29) = count(29) + 1;
%     end
%     if((log10(W_U_T_Jan(i)) >= 20.0) && (log10(W_U_T_Jan(i)) < 20.3))
%         count(30) = count(30) + 1;
%     end
end
%histo = histc(log10(W_U_T_Jan),11:0.3:20.3);
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);
%bin_counts = histc(log10(W_U_T_Jan),1.25:0.5:11.75)
% 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)
%loglog(bin,histo,'ro','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