clear all
close all


 disp('Auditory Cortex Spikes')
 lw = 3; 
 set(0, 'DefaultAxesFontSize', 16);
 fs = 15;
 msize = 10;

% The data set comes from experiments in the lab of Dr. Robert Liu
% (BME-Emory: http://www.biology.emory.edu/research/Liu/index.html) 
% and concerns single unit electrophysiology in the auditory cortex of 
% non-anesthetized female mice. The motivating question is the exploration 
% of auditory neural differences between female parents vs female virgins
% and its relationship to cortical response.
% 
% All the data represents single neuron responses ( i.e. spikes) from the 
% same animal. However, the data corresponds to multiple single neuron
% recordings in the same animal. Multiple trials are performed on the single
% neurons from one animal and thus to reiterate, the data reflects multiple
% single neurons from the same animal (one mother and one naive).
% 
%      The recordings are made from a region in the auditory cortex of the 
% mouse with a single tungsten electrode. A sound stimulus is presented 
% at the time of 200 ms during each sweep (time shown is 0-611 and 200
% is the point at which a stimulus is presented).  Each sweep is 611 ms
% long and the duration of the tone is 10-70 ms.
% 
%     The data collected consist of several thousand trials and with varying 
% ranges of frequencies. Through descriptive statistics we may be able to
% infer inherent differences that exist in the data between mothers and naive.

% %
 clear all
 close all
    disp('Univariate Data: Auditory Cortex  Spikes')
    lw = 2; 
    set(0, 'DefaultAxesFontSize', 14);
    fs = 15;
    msize = 10;

load 'C:/STAT/Descriptive/Descriptivedat/spikes.dat'

% average number of spikes (over = 6 x 1850 trials) in the
% intervals of 1 ms, 2ms, ..., 611ms
%
ms = spikes(:,1);
qm = spikes(:,2);  %cumulative -- mother
qn = spikes(:,3);   %cumulative -- naive
figure(1)
   plot(ms, qm)
   hold on
   plot(ms, qn, 'r-')
   legend('Mother','Naive','Location','NorthWest')
   axis tight
%-------------------------------------------------
dm = diff(qm,1);  dn=diff(qn,1);  %first differences
dm=dm(1:512); dn=dn(1:512);  %cut  sample size to a nice power of 2
figure(2)
   plot(dm)
   hold on
   plot(dn,'r-')
   legend('Mother','Naive','Location','NorthEast')
   axis tight
%--------------------------------------------------
%--------------------------------------------------
data = dn;
%---------
% measures of location
%---------
format long
me   = mean(data)
gme = geomean(data)
hme = harmmean(data)
med = median(data)
%  trimmean(data, 100) gives error, but
tme1 = trimmean(data, 99.99)  % gives the median
tme2 = trimmean(data, 0)   %gives the ordinary mean -- 0 trimming
[mo fre] = mode(data)

%------------------
% measures of variation
%------------------
stde = std(data)         %by default flag=0, i.e., std(data,0)
stde1 = std(data, 1)  %flag 1 stands for division by n
vari = var(data, 1)
ran = range(data)
mad0 = mad(data)
mad1 = mad(data,1)
realmad = 1/0.6745 * mad(data,1)   % compare to std(data)
%-----------------------
% missing values NaN
%-----------------------
datanan=data;
datanan([5 7 101 418])= NaN;
men1 = mean(datanan)
men2 = nanmean(datanan)
nmax = nanmax(datanan)
nsum = nansum(datanan)
%-------------------
qts = quantile(data,[0.1 0.25 0.5 0.75 0.9])
prts = prctile(data,[10 25 50 75 90])
%
figure(3)
% boxplot of data
boxplot(data,'notch','on')
%print -depsc 'spikes1.eps'
%---------------------------
% Sample Moments of Higher Order
%kth (row) sample moment: mean(data.^k)
 mo3=mean(data.^3) % third moment
 %kth central moment mean( (data-mean(data)).^k )
 cmo3=mean( (data-mean(data)).^3 )
%-----------
% sample skewness mean(data.^3)/std(data,1)^3
skdef = mean( (data-mean(data)).^3 )/std(data,1)^3
skmlab=skewness(data)
% sample kurtosis
kudef=mean( (data-mean(data)).^4 )/std(data,1)^4
kumlab= kurtosis(data)

% sample CV [coefficient of variation]
cv = std(data)/mean(data)
format short

figure(4)
% histogram with 30 bins
hist(data, 30)
%print -depsc 'spikes2.eps'


figure(5)
[f,x,u] = ksdensity(data);
plot(x,f,'linewidth',lw)
title('Density estimate for data')
hold on
[f,x] = ksdensity(data,'width',u/3);
plot(x,f,'r','linewidth',lw);
[f,x] = ksdensity(data,'width',u*3);
plot(x,f,'g','linewidth',lw);
legend('default width','1/3 default','3*default')
hold off
%print -depsc 'spikes3.eps'
 figure(6)
[f,xf] = ecdf(data);
stairs(xf,f)
title('Empirical Cumulative Distribution Function (CDF) for Spikes')
%print -depsc 'spikes4.eps'
%-------------------------------------------------
% Autocorrelations in 1-D data; Comparison naive/mothers
%
figure(7)
k=8;  %number of lags for autocorrelations
an = acorr(dn, k+1)  %need simple functions autocorr, autocov
am = acorr(dm, k+1)
nmac = [an ; am];
g = [0 1]; % labels for naive/mothers
nID = [1:2]; % Actual labels required on the x axis;

bar3(0:k,nmac') ; % Plot the 3D bar

% Put  labels
xlabel('naive/mothers', 'FontSize', 14);
ylabel('lags', 'FontSize', 14);
zlabel('autocorrelation', 'FontSize', 14);
title('Autocorrelations for Naive/Mother Mice', 'FontSize', 14);

% Set custom ticks and tick labels
set(gca, 'Xtick', nID);
set(gca, 'xTickLabel', g(nID));
colormap('white');
axis tight;





