Линия наилучшего соответствия с логарифмическими осями (MATLAB)

Я пытаюсь построить наиболее подходящую линию на функции плотности вероятности с логарифмическими осями. Ось Y (PDF) составляет от 10 ^ -12 до 10 ^ -28, а ось X - от 10 ^ 10 до 10 ^ 20. Я пробовал полифит, но безуспешно. Любые идеи? Прилагаю мой код.

Спасибо, Кевин

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

person user2639937    schedule 21.08.2013    source источник
comment
Часть бункера может быть уменьшена на histc(log10(W_SH),11:0.3:20.3).   -  person marsei    schedule 22.08.2013
comment
Будет ли polyfit с log10(bin) и log10(prob_dens_mean) лучшим вариантом для поиска подходящих коэффициентов? Сначала я думаю о линейном графике, а затем о журнале.   -  person marsei    schedule 22.08.2013
comment
Если обе оси правильно отображаются при преобразовании журнала, то есть ли веская причина, по которой вы не использовали polyfit (log (x), log (y), 1)?   -  person    schedule 22.08.2013
comment
Привет, я попробовал приведенный выше код как с best_fit = polyfit (log10 (bin), log10 (prob_dens_mean), 1), так и с best_fit = polyfit (log (bin), log (prob_dens_mean), 1), и ни один из графиков не подходит лучше всего линия. Вместо этого я вижу, что best_fit - это NaN. Любые идеи?   -  person user2639937    schedule 23.08.2013


Ответы (1)


У меня нет ваших данных, но вот пример с использованием случайных случайных данных с нормальным распределением

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)
person outdoor_guy    schedule 27.01.2014