2013-08-29 134 views
0

想要繪製具有對數座標軸的PDF上的最佳擬合線。 x軸爲10 至10 千兆瓦。 y軸是10 -2至10 ,並且是概率分佈。到目前爲止,我曾嘗試:使用對數座標軸在PDF圖上繪製最佳擬合線

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 

回答

1

我以爲你是因爲混合線性空間和日誌空間的問題。嘗試安裝和日誌空間繪製這樣的:

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

此外,我建議使用polyval,而不是計算linePointsY自己。

相關問題