필터 지우기
필터 지우기

Best fit line on PDF plot with logarithmic axes

조회 수: 1 (최근 30일)
Kevin
Kevin 2013년 8월 29일
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

답변 (0개)

카테고리

Help CenterFile Exchange에서 Line Plots에 대해 자세히 알아보기

Community Treasure Hunt

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

Start Hunting!

Translated by