Составьте гистограмму каждого столбца матрицы, но без нулевых элементов - PullRequest
Новая
       18

Составьте гистограмму каждого столбца матрицы, но без нулевых элементов

0 голосов
/ 12 мая 2019

Я запускаю симуляцию частиц в коробке. Когда частица покидает ящик, ее кинетическая энергия становится равной нулю (на время> t убегает). Поэтому я хочу составить гистограмму того, как Wkinet (который является функцией nP = числа частиц, ntM = временных шагов) эволюционирует во времени, но я не хочу принимать во внимание нулевые значения каждого столбца. Есть ли способ кодировать его, чтобы он мог найти оптимальное количество бинов?

Вот что я пробовал: 

nbins = 1000;
for j = 1:ntM/5
    Wkinet(Wkinet==0) = NaN;
    y = Wkinet(:,j).*erg2eV;
end
histfit(y,nbins)

here

1 Ответ

1 голос
/ 13 мая 2019

Логическое индексирование часто бывает довольно быстрым и интуитивно понятным, когда вы освоите синтаксис. 

myTolerance=1e-7; % in erg units.
nbins=1000;
for j=1:ntM/5
    %Wkinet(Wkinet==0)=NaN;
    % y=Wkinet(:,j).*erg2eV; % An extra assigment is costly and probably not needed.
    H = histfit(Wkinet(abs(Wkinet(:,j))>myTolerance, j) * erg2ev, nbins);
    % Select from column j all rows in column j whose absolute values are greater than the tolerance. 
    % Assumption; erg2ev is just a scalar, otherwise select its entries with erg2ev(abs(Wkinet(:,j))>myTolerance)
    H(1).delete; Remove bins, only keep the fit.
    set(gca, 'YScale', 'log'); % Make logarithmic Y
    set(gca, 'XScale', 'log'); % Make logarithmic X
    pause
end

Если вам нужно явное ограничение оси, используйте 

xlim([lowerBound upperBound]); ylim(etc...

... или иногда полезно использовать команду оси для точного управления, например, 

ax=axis; ax(3)=min( 8ax(3) maxAllowedY]); axis(ax);

«Пауза» (для интерактивного использования) может быть заменена командой печати для сохранения графиков на диск. Э.Г. 

print(sprintf('My_plot_%02d',j),'-dpng');

Или сохранить фигуру: 

savefig(sprintf('My_fig_%02d',j));

Если вы уверены, что количество графиков меньше, чем, скажем, 16, вы можете поместить в цикл команду subplot. Заменить паузу на 

subplot(4,4,j);

Конечная нота; если вы намереваетесь построить нормальное распределение, соответствующее вашим ненулевым данным, вы можете получить лучшие результаты, заменив функцию histfit с помощью 

myFit = fitdist(Wkinet(Wkinet(:,j)>myTolerance, j) * erg2ev), 'Normal');
maxEv = max(Wkinet(Wkinet(:,j)>myTolerance, j) * erg2ev);
myX   = [myTolerance; maxEv/100; maxEv]; % Alter for different plot X-axis
myY   = pdf(myFit, myX);
plot(myX, myY);

Я проверил, и - это разница между fitdist иstistist, вероятно, вызванная дискретизацией бина.

...