Clear Filters
Clear Filters

Best fit line on PDF plot with logarithmic axes

1 view (last 30 days)
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

Answers (0)

Categories

Find more on Line Plots in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!