2013-08-21 49 views
2

我想绘制对数坐标轴的概率密度函数的最佳拟合线。 Y轴(PDF)为10^-12至10^-28,而X轴为10^10至10^20。我尝试了polyfit,没有运气。有任何想法吗?附件是我的代码。对数轴的最佳拟合线(MATLAB)

感谢, 凯文

clc; 
clear all; 

load Aug2005_basin_variables.mat 

% Initialize 

j_len = length(W_SH); 
prob_dens_all = zeros(j_len,30); 
ii = 1 : j_len; 
count(1:30) = 0; 
bin(1:30) = 0; 

for i = 1 : 30 
    bin(i) = 10^(11 + (0.3*i)); 
end 


% Bin the Watts 

for i = 1 : j_len 
    if((log10(W_SH(i)) >= 11) && (log10(W_SH(i)) < 11.3)) 
     count(1) = count(1) + 1; 
    end 
    if((log10(W_SH(i)) >= 11.3) && (log10(W_SH(i)) < 11.6)) 
     count(2) = count(2) + 1; 
    end 
    if((log10(W_SH(i)) >= 11.6) && (log10(W_SH(i)) < 11.9)) 
     count(3) = count(3) + 1; 
    end 
    if((log10(W_SH(i)) >= 11.9) && (log10(W_SH(i)) < 12.2)) 
     count(4) = count(4) + 1; 
    end 
    if((log10(W_SH(i)) >= 12.2) && (log10(W_SH(i)) < 12.5)) 
     count(5) = count(5) + 1; 
    end 
    if((log10(W_SH(i)) >= 12.5) && (log10(W_SH(i)) < 12.8)) 
     count(6) = count(6) + 1; 
    end 
    if((log10(W_SH(i)) >= 12.8) && (log10(W_SH(i)) < 13.1)) 
     count(7) = count(7) + 1; 
    end 
    if((log10(W_SH(i)) >= 13.1) && (log10(W_SH(i)) < 13.4)) 
     count(8) = count(8) + 1; 
    end 
    if((log10(W_SH(i)) >= 13.4) && (log10(W_SH(i)) < 13.7)) 
     count(9) = count(9) + 1; 
    end 
    if((log10(W_SH(i)) >= 13.7) && (log10(W_SH(i)) < 14.0)) 
     count(10) = count(10) + 1; 
    end 
    if((log10(W_SH(i)) >= 14.0) && (log10(W_SH(i)) < 14.3)) 
     count(11) = count(11) + 1; 
    end 
    if((log10(W_SH(i)) >= 14.3) && (log10(W_SH(i)) < 14.6)) 
     count(12) = count(12) + 1; 
    end 
    if((log10(W_SH(i)) >= 14.6) && (log10(W_SH(i)) < 14.9)) 
     count(13) = count(13) + 1; 
    end 
    if((log10(W_SH(i)) >= 14.9) && (log10(W_SH(i)) < 15.2)) 
     count(14) = count(14) + 1; 
    end 
    if((log10(W_SH(i)) >= 15.2) && (log10(W_SH(i)) < 15.5)) 
     count(15) = count(15) + 1; 
    end 
    if((log10(W_SH(i)) >= 15.5) && (log10(W_SH(i)) < 15.8)) 
     count(16) = count(16) + 1; 
    end 
    if((log10(W_SH(i)) >= 15.8) && (log10(W_SH(i)) < 16.1)) 
     count(17) = count(17) + 1; 
    end 
    if((log10(W_SH(i)) >= 16.1) && (log10(W_SH(i)) < 16.4)) 
     count(18) = count(18) + 1; 
    end 
    if((log10(W_SH(i)) >= 16.4) && (log10(W_SH(i)) < 16.7)) 
     count(19) = count(19) + 1; 
    end 
    if((log10(W_SH(i)) >= 16.7) && (log10(W_SH(i)) < 17.0)) 
     count(20) = count(20) + 1; 
    end 
    if((log10(W_SH(i)) >= 17.3) && (log10(W_SH(i)) < 17.6)) 
     count(21) = count(21) + 1; 
    end 
    if((log10(W_SH(i)) >= 17.6) && (log10(W_SH(i)) < 17.9)) 
     count(22) = count(22) + 1; 
    end 
    if((log10(W_SH(i)) >= 17.9) && (log10(W_SH(i)) < 18.2)) 
     count(23) = count(23) + 1; 
    end 
    if((log10(W_SH(i)) >= 18.2) && (log10(W_SH(i)) < 18.5)) 
     count(24) = count(24) + 1; 
    end 
    if((log10(W_SH(i)) >= 18.5) && (log10(W_SH(i)) < 18.8)) 
     count(25) = count(25) + 1; 
    end 
    if((log10(W_SH(i)) >= 18.8) && (log10(W_SH(i)) < 19.1)) 
     count(26) = count(26) + 1; 
    end 
    if((log10(W_SH(i)) >= 19.1) && (log10(W_SH(i)) < 19.4)) 
     count(27) = count(27) + 1; 
    end 
    if((log10(W_SH(i)) >= 19.4) && (log10(W_SH(i)) < 19.7)) 
     count(28) = count(28) + 1; 
    end 
    if((log10(W_SH(i)) >= 19.7) && (log10(W_SH(i)) < 20.0)) 
     count(29) = count(29) + 1; 
    end 
    if((log10(W_SH(i)) >= 20.0) && (log10(W_SH(i)) < 20.3)) 
     count(30) = count(30) + 1; 
    end 
end 


for i=1:30 
    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,30); 


for i = 1 : 30 
    prob_dens_mean(1,i) = mean(prob_dens_all(:,i)); 
    %prob_dens_std(1,i) = std(prob_dens_all(:,i)); 
end 

% Plot 

best_fit = polyfit(bin,log10(prob_dens_mean),11) 

h = figure; 
loglog(bin,prob_dens_mean,'ro','MarkerSize',10) 
hold on; 
plot(best_fit,'b') 
t = title('Event Power Distribution, SHem, August 2005'); 
set(t, 'FontWeight', 'bold', 'FontSize', 12) 
set(gca, 'FontWeight', 'bold', 'FontSize', 12) 
xlabel('Event Power (W)'); 
ylabel('Probability Density'); 
print -dpng SHem_Wattage_PDF_AUG2005.png 
+1

纸槽部分可以通过'减少的例子histc(LOG10(W_SH),11:0.3:20.3)'。 – marsei

+0

'log10(bin)''''log10(prob_dens_mean)''''polyfit''是找到一些拟合系数的更好选择吗?我首先想到了一个线性图,然后是一个loglog。 – marsei

+0

如果两个轴在日志转换中都正确,那么为什么不使用polyfit(log(x),log(y),1)是否有很好的理由? – 2013-08-21 23:49:57

回答

0

我没有数据,但这里是用一些随机正态分布的随机数据

x=randn(1000,1)+5; % create some data, keep numbers positive by adding 5 
[n,xb]=hist(x); % Create the histogram 
n = n/sum(n); % convert counts to a pdf 
p=polyfit(log(xb), log(n), 3); % Do a 3rd order fit 
loglog(xb,n, '*-', xb, exp(polyval(p, log(xb))), 'r') 
grid on 
legend('PDF', 'Fit', 0)