CODEWITHSUNDEEP
matlabplotgraphhistogram
wololo
wololo

Reputation: 41

MATLAB Histogram Problems

I need to calculate the probability density of the function 1/2-x where x is a randomly generated number between 0 and 1.

My code looks like this:

n=10^6;
b=10^2;

x=(1./(2-(rand(1,n))));

a=histfit(x,b,'kernel'); % Requires Statistics and Machine Learning Toolbox
xlim([0.5 1.0])

And I get a decent graph that looks like this:

enter image description here

As it may be apparent, there are a couple of problems with this:

  1. MATLAB draws a fit that differs from my histogram, because it counts in the empty space outside the [0.5 1] range of the function as well. This results in a distorted fit towards the edges. (The reason you don't see said empty space is because I entered an xlim there)

  2. I don't know how I could divide every value in the Y-axis by 10^6, which would give me my probability density.

Thanks in advance.

Upvotes: 4

Views: 380

Answers (2)

rozsasarpi
rozsasarpi

Reputation: 1641

1. Better result can be obtained by using ksdensity and specifying the support of the distribution.

2. By using hist you have access to the counts and centers, thus the normalization to get density is straightforward.

Code to demonstrate the suggestions:

rng(125)
n=10^6;
b=10^2;

x=(1./(2-(rand(1,n))));

subplot(1,2,1)
a = histfit(x,b,'kernel');
title('Original')
xlim([0.5 1.0])

[f,c] = hist(x,b);
% normalization to get density
f = f/trapz(c,f);

% kernel density
pts = linspace(0.5, 1, 100);
[fk,xk] = ksdensity(x, pts, 'Support', [0.5, 1]);

subplot(1,2,2)
bar(c,f)
hold on
plot(xk,fk, 'red', 'LineWidth', 2)
title('Improved')
xlim([0.5 1.0])

Comparing the results: enter image description here

EDIT: If you do not like the endings:

pts = linspace(0.5, 1, 500);
[fk,xk] = ksdensity(x, pts, 'Support', [0.5, 1]);
bar(c,f)
hold on
plot(xk(2:end-1),fk(2:end-1), 'red', 'LineWidth', 2)
title('Improved_2')
xlim([0.5 1.0])

enter image description here

Upvotes: 2

GameOfThrows
GameOfThrows

Reputation: 4510

To solve both of your problems, I suggest using hist (Note if you have a version above 2010b, you should use histogram instead) instead of histfit to first get the values of your histogram and then doing a regression and plotting them:

n=10^6;
b=10^2;
x=(1./(2-(rand(1,n))));

[counts,centers]=hist(x,b);
density = counts./trapz(centers, counts);
%// Thanks to @Arpi for this correction 
polyFitting = polyfit(centers,density,3)
polyPlot = polyval(polyFitting,centers)
figure
bar(centers,density)
hold on
plot(centers,polyPlot,'r','LineWidth',3)

You can also up the resolution by adjusting b, which is set to 100 currently. Also try different regressions to see which one you prefer.

result looks like

Upvotes: 3

Related Questions