%%%%%%%%%%%%%%%%%%%%%%%%
%%%%  STATISTICS  %%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%


=====================================
====  descriptive statistics  =======
=====================================

N = 100;
x = linspace(0,10,N);
data = round(x.^2 - 0.01*x.^4 +6*randn(1,N));
plot(x,data,'o','linewidth',2);

hist(data,10)            % histogram

data_mean = mean(data);
data_med = median(data);
data_mod = mode(data);

data_range = range(data);  % = max(data)-min(data)
data_iqr = iqr(data);    % interquartile range (a measure of dispersion)

data_std = std(data);    % standard deviation
data_var = var(data);    % data variance
help var  % higher moments: skewness, ...

data_quant = quantile(data);    % data min, quartiles, max
p = [0.25 0.5 0.75];
data_quart = quantile(data, p); % data quartiles

 % min, 1. quartile, med, 3. quart., max, mean, std, skewness, kurt.:
statistics(data);


=====================================
======    polynomial fitting    =====
=====================================

N = 100;
x = linspace(-4,10,N);
f=@(x)  -0.08*x.^3 + 0.8*x.^2 + 2*x;
plot(x,f(x),'linewidth',2)
hold on

data = f(x) +6*randn(1,N);
plot(x,data,'o','linewidth',2);

p = polyfit(x, data, 1);     % polynomial of 1-st order
y = polyval(p, x);           % evaluation at x
plot(x,y,'k','linewidth',2);
p = polyfit(x, data, 2);     % 2-nd order
y = polyval(p, x);
plot(x,y,'m','linewidth',2);
p = polyfit(x, data, 3)      % 3-rd order
y = polyval(p, x);
plot(x,y,'r','linewidth',2);


=====================================
======    correlation    ============
=====================================

N = 100;
x = linspace(-4,10,N);
f=@(x) -0.08*x.^3 + 0.8*x.^2 + 2*x;
data = round(f(x) +6*randn(1,N));
plot(x,data,'o','linewidth',2);
hold on
corr(x,data)

data1 = round(f(x) +6*randn(1,N));
plot(x,data1,'ro','linewidth',2);

corr(data,data1)
corr(data,data)
corr(data,-data)


=====================================
random numbers generation, histogram
=====================================

y = randn(10000,1);  % standard normal distribution
hist(y,9)
hist(y,100)

y = rand(10000,1);   % uniform dist.
hist(y,100)

y = randi(20, 10000,1);  % integers in <1,20>, uniform. dist.
hist(y,20)

y = rande(10000,1);  % exponential dist.
hist(y,100)

y = randp(2.3, 10000,1); % Poisson dist. with mean 2.3
hist(y,10)

----------------------

v=[1 2 3 4 5];                   % sampled values
p = [1 2 4 2 1];                 % relative probabilities
y = discrete_rnd (v, p, 1,1000); % sample of given size
hist(y,5)

---------------------

v = randn(1000,1);              % given data
hist(v,20)
   % more data generated from the same distr.:
y = empirical_rnd (v, 1,10000);
figure;
hist(y,20)

------------------