%----------------- Statistical Data Processing -----------------
%
% Technical University of Sofia - MUSSPI lab course
% December 2016
%---------------------------------------------------------------

%% Part I

close all;
clear all;

filename = 'til';
extension = 'jpg';

A = imread(filename,extension);
A1 = A(:,:,1);
A11 = double(A1);
A11 = reshape(A11,size(A1,1)*size(A1,2),1);
X = A11;

% data 1
X=sort(X);
n=length(X);
xmin=X(1);
xmax=X(n);
k=round(n^0.5);     % number of intervals
d=(xmax-xmin)/k;    % width of a single interval
del=(xmax-xmin)/20; % a bit more width for ease of observation
xl=xmin-del;        
xr=xmax+del;
fprintf('Number of intervals k=%d\n Interval h=%14.7f\n',k,d)
hist(X,100);
title('Histogram');
xlim([xl xr]);

Mx=mean(X);
Sx=std(X);
Dx=Sx^2;
Ax=sum((X-Mx).^3)/(n-1)/Sx^3;
Ex=sum((X-Mx).^4)/(n-1)/Dx^2-3;
fprintf('Expected value: Mx=%14.7f\n',Mx);
fprintf('Standard deviation: Sx=%14.7f\n',Sx);
fprintf('Dispersion: Dx=%14.7f\n',Dx);
fprintf('Asymmetry: Ax=%14.7f\n',Ax);
fprintf('Excess: Ex=%14.7f\n',Ex);
pause(20);
%% Part II

clear all;
close all;

load('danni.mat');
X=sort(X);
n=length(X);
xmin=X(1);
xmax=X(n);
k=round(n^0.5);     % number of intervals
d=(xmax-xmin)/k;    % width of a single interval
del=(xmax-xmin)/20; % a bit more width for ease of observation
xl=xmin-del;        
xr=xmax+del;
[nj,xm]=hist(X,k);   
delta=xm(2)-xm(1);  % width of the interval (difference between two centers)
hist(X,k);
title('Histogram');
xlim([xl xr]);


xfv=[xm-delta/2; xm+delta/2];  %x-axis for the emperical f(x)
xfv=reshape(xfv,prod(size(xfv)),1); % reshape into a column
xfv=[xl;xfv(1);xfv;xfv(end);xr];  % add the corner elements
fv=nj/(n*delta);  % values of the empirical f(x) in one row
fv=[fv;fv];   % 2 rows
fv=[0;0;reshape(fv,prod(size(fv)),1);0;0];   % +corner elements
xft=linspace(xl,xr,1000);  % x-axis for the theoretical f(x)


ft=[normpdf(xft,mean(X),std(X));exppdf(xft,mean(X));unifpdf(xft,xmin,xmax); raylpdf(xft,mean(X)*sqrt(2/pi))];
col='bgmy'; %plot colors
plot(xfv,fv,'-r', xft,ft(1,:),col(1), xft,ft(2,:),col(2), xft,ft(3,:),col(3), xft,ft(4,:),col(4))
legend('Data','Normal pdf','Exp. pdf','Cont. Uniform pdf','Rayleigh pdf');
xlabel('x');
ylabel('f(x)');
xlim([xl xr]);
ylim([0 1.4*max(fv)]);