%% Partial data from the Lab of Dr Todd McDevitt, BME
%Area of grown cell [in \mu m^2] after time T=2 under
%the static condition.
clear all
close all force
% 
% Preamble
     disp('Univariate Data: Cells from the Lab of Todd McDevitt')
     lw = 3; 
     set(0, 'DefaultAxesFontSize', 16);
     fs = 15;
     msize = 10;
% 
% load('C:\your path\cellarea.dat')
load('C:\STAT\Descriptive\Descriptivedat\cellarea.dat')
% 
figure(1)
hist(cellarea, 100)
%print -depsc 'C:\Courses\bmestatu\cells0.eps'

%%
%  
sortedca = sort(cellarea);
sortedca(1)
   %ans=7000
sortedca(end)
   %ans = 1260000
   %----------------------------
% eliminate maximal element that looks like outlier  
car =cellarea(cellarea ~= max(cellarea)); 
car = car/1000; % scale measurements to nice numbers
                % now areas are in 1000's \mu m^2
               figure(2)
hist(car, 100)           
n = length(car)
    % n = 462
    % 
%-------------------------------------
%%
% 
%Location Measures
location = [geomean(car) harmmean(car) mean(car) ...
    median(car) mode(car) trimmean(car,20)]
%location = 18.8485   15.4211   24.8701   17.0000   10.0000   19.7298

%
%also Winsorized Mean
 alpha=0.20;
  sa = sort(car); 
  sa(1:floor( n*alpha/2 )) = sa(floor( n*alpha/2 ) + 1);
  sa(end-floor( n*alpha/2 ):end) = sa(end-floor( n*alpha/2 ) - 1);
  winsmean = mean(sa)
% winsmean = 21.9632
 
%%
%Variability Measures 
 var(car) % sample variance
   %ans = 588.9592
std(car) %s.  stdev, sum divided by (n-1)
   % ans = 24.2685
 std(car,1) %s. stdev, sum divided by n
   % ans = 24.2422
 sqrt(var(car)) %should be equal to std
   % ans = 24.2685
mad(car)  % mean absolute deviation from the mean;
         % MAD is usually median absolute deviation from the median
         % ans = 15.3328     
realmad = 1/0.6745 * median( abs(car - median(car)))
         % real mad in MATLAB is 1/0.6745 * mad(car,1)
         %realmad = 10.3781
         
range(car) %Range, span of data, Max - Min
         % ans = 212
iqr(car) % inter-quartile range Q3-Q1
         % ans = 19
%--------------------------------------------

%Sample Quantiles (Sample Percentiles)
% For 5%, 10%, 25%, 50%, 75%, 90%, 95% percentiles
% are:
prctile(car, [5, 10, 25, 50, 75, 90, 95] )
 %ans =  7     8    11    17    30    51    67

quantile(car, [0.05, 0.10, 0.25, 0.50, 0.75, 0.90, 0.95] )
 %ans =  7     8    11    17    30    51    67


%1st quartile Q1 --> prctile(car, 25)
%3rd quartile Q3 --> prctile(car, 75)
% median --> prctile(car, 50)
 
prctile(car, 75)- prctile(car, 25) %should be equal to iqr
     % ans = 19
     
 % z-scores, standardize the sample
 za = zscore(car);
 mean(za)
   % ans = -5.8155e-17
 var(za)
   % ans = 1.0000
  
%%
% Moments of Higher Orders
%kth (row) moment: mean(car.^k)
mean(car.^3) % third  moment
      % ans = 1.1161e+05
%kth central moment mean( (car-mean(car)).^k )
mean( (car-mean(car)).^3 )
      % ans = 5.2383e+04
% sample skewness mean(car.^3)/std(car,1)^3
mean( (car-mean(car)).^3 )/std(car,1)^3
   %ans = 3.6769
skewness(car)
   %ans = 3.6769
%
% sample kurtosis 
mean( (car-mean(car)).^4 )/std(car,1)^4
    % ans = 22.8297
kurtosis(car)
    % ans = 22.8297
    
% %robust estimators   
% sk = (prctile(car, 75)+prctile(car, 25) - ...
%  2 * median(car))/(prctile(car, 75) - prctile(car, 25))
%    %0.3684 
% 
% %robust kurtosis 
% kur = (prctile(car,7/8*100)-prctile(car,5/8*100)+prctile(car,3/8*100)- ...
%   prctile(car,1/8*100))/(prctile(car,6/8*100)-prctile(car,2/8*100))
%   %1.4211     
    
% sample CV [coefficient of variation]
std(car)/mean(car)
   % ans = 0.9758

%%
%----------- Some Graphical Summaries of the Sample
%
figure(3)
boxplot(car)
%print -depsc 'C:\Springer\Descriptive\Descriptiveeps\cells1.eps'
%%
figure(4)
hist(car, 80)

%print -depsc 'C:\Springer\Descriptive\Descriptiveeps\cells2.eps'
%%
[f,x,u] = ksdensity(car );
figure(5)
plot(x,f, 'linewidth',2)
hold on
[f,x] = ksdensity(car,'width',u/3);
plot(x,f,'r','linewidth',2);
[f,x] = ksdensity(car,'width',u*3);
plot(x,f,'g', 'linewidth',2);
legend('default width','default/3','3 * default')
plot(car, zeros(size(car)), 'b+');
hold off
%print -depsc 'C:\Springer\Descriptive\Descriptiveeps\cells3.eps'
title('Density estimate for Cell Areas') 
%%
figure(6)
 xx = min(car)-1:0.01:max(car)+1;
 yy = empiricalcdf(xx, car);
 plot(xx, yy, 'k-','linewidth',2)
 xlabel('x');  ylabel('F_n(x)')
%print -depsc 'C:\Springer\Descriptive\Descriptiveeps\cells4a.eps'



figure(7)
[f,xf] = ecdf(car);
stairs(xf,f)
%print -depsc 'C:\Springer\Descriptive\Descriptiveeps\cells4.eps'
title('Empirical Cumulative Distribution Function (CDF) for Cell Areas') 
%%

figure(8)
qqplot(car);
%% 
%print -depsc 'C:\Springer\Descriptive\Descriptiveeps\cells5.eps'
%%
figure(9)
%pie plot proportions
n1 = sum( car <= 15 );
n2 = sum( (car > 15 ) & (car <= 30) );
n3 = sum( car > 30 );
% part of a circle n1/n, n2/n and n3/n
explode = [0 1 0];
pie([n1, n2, n3], explode)
colormap summer
%print -depsc 'C:\Springer\Descriptive\Descriptiveeps\cellspie1.eps'

figure(10)
pie3([n1, n2, n3], explode)
colormap summer
%print -depsc 'C:\Springer\Descriptive\Descriptiveeps\cellspie2.eps'



