function null=discriminate()

%DISCRIMINATE
%This function generates a null expectation for the jack-knife
%reclassification success of a discriminant function analysis by
%randomizing user-provided data

%Copyright 2005 J. Wilson White
%Dept. of Ecology, Evolution, and Marine Biology
%University of California, Santa Barbara, CA 93106
%w_white@lifesci.ucsb.edu

%version: last modified March 8, 2007

%Input:  
%        Excel file DataImport.xls (must be completed and placed in same
%        directory as this .m file)

%Output: Vector 'null' containing (in order):
%        Mean of randomized distribution (null(1))
%        Standard deviation of randomized distribution (null(2))
%        Median of randomized distribution (null(3))
%        95% quantile of randomized distribution (null(4))
%        5% quantile of randomized distribution  (null(5))
%        Actual jackknife value (null(6))
%        One-tailed P-value (alpha = 0.05) for actual jackknife value
%        (p(7))

%This program uses linear DFA with uniform priors by default
%To change this setting, alter the input arguments to function 'classify'
%on line 82 (for more information see the Matlab help file for 'classify')

D = xlsread('DataImport.xls');  %Loads actual data to be used
i = 1; while isnan(D(i,2)), i = i+1;end %finds beginning of data
runs = D(i,2);                 %Number of bootstraps desired
D = D(i+4:end,:);               %Truncates empty cells out of data matrix


Groupid = D(:,1);               %vector Groupid holds group labels
Data = D(:,2:end);              %matrix Data holds raw data
vars = length(D(1,2:end));      %vars = #columns in Data
g = length(Groupid);           %g = #rows in Data


rdist=[];

%Begin randomization
for z=1:runs   
    
    %Assemble randomized dataset:
    newgroups = randperm(g);
    for b = 1:g %for each sample
        rgroups(b) = Groupid(newgroups(b)); %select random group ID
    end    %end for each sample
    

%the following loop performs the jackknife 
for i=1:length(Data(:,1)) %begin jackknifing
    data2=Data;
    group2=rgroups;
    group2(i)=NaN; %delete the group ID to be jackknifed out
    data2(i,:)=NaN; %delete the data to be jackknifed out
    sample=Data(i,:);
    class(i)=classify(sample,data2,group2,'linear'); %Matlab's DFA function
end %end jackknifing

q = 0;  %q is counter for DFA classification success
q = sum(class == rgroups); %tabulate correct classifications

dis=q/length(Data(:,1));    %calculate jackknife success rate

rdist(z)=dis;%store success rate in vector 'rdist'

end %end randomizations


%Now calculate observed jackknife success rate for the data:
 
for i=1:length(Data(:,1))
    data2=Data;
    group2=Groupid;
    group2(i)=NaN;
    data2(i, :) = NaN;
    sample=Data(i,:);
    class(i)=classify(sample,data2,group2,'linear');
end
  

q = 0; 
q = sum(class' == Groupid); %tabulate correct classifications
Uservalue=q/length(Groupid);    %calculate observed jackknife success rate

%Calculate Pvalue

rdist = sort(rdist); %sort rdist in ascending order

i = 1;
while rdist(i) < Uservalue && i < runs %determine where actual jackknife value lies in null distribution
    i = i + 1;
end

Pvalue = 1-(i)/runs; %one-tailed P-value based on position of actual jackknife value in null distribution
if Pvalue == 0       
    Pvalue = 1/runs;  %the smallest possible P-value provided is 1/runs
end

   

null=[mean(rdist) std(rdist) quantile(rdist,0.5) quantile(rdist,0.95) quantile(rdist,0.05) Uservalue Pvalue];
%Report results
fprintf('\n','Results:','\n')
Mean=null(1)
sd = null(2)
Median = null(3)
Upper95quantile = null(4)
Lower95quantile = null(5)
Observed = Uservalue
P = Pvalue

end

function y = quantile(x,p,dim);
%QUANTILE Quantiles of a sample.
%   Y = QUANTILE(X,P) returns quantiles of the values in X.  P is a scalar
%   or a vector of cumulative probability values.  When X is a vector, Y is
%   the same size as P, and Y(i) contains the P(i)-th quantile.  When X is
%   a matrix, the i-th row of Y contains the P(i)-th quantiles of each
%   column of X.  For N-D arrays, QUANTILE operates along the first
%   non-singleton dimension.
%
%   Y = QUANTILE(X,P,DIM) calculates quantiles along dimension DIM.  The
%   DIM'th dimension of Y has length LENGTH(P).
%
%   Quantiles are specified using cumulative probabilities, from 0 to 1.
%   For an N element vector X, QUANTILE computes quantiles as follows:
%      1) The sorted values in X are taken as the (0.5/N), (1.5/n),
%         ..., ((N-0.5)/N) quantiles.
%      2) Linear interpolation is used to compute quantiles for
%         probabilities between (0.5/N) and ((N-0.5)/N)
%      3) The minimum or maximum values in X are assigned to quantiles
%         for probabilities outside that range.
%
%   QUANTILE treats NaNs as missing values, and removes them.
%
%   Examples:
%      y = quantile(x,.50); % the median of x
%      y = quantile(x,[.025 .25 .50 .75 .975]); % a useful summary of x
%
%   See also PRCTILE, IQR, MEDIAN.

%   Copyright 1993-2004 The MathWorks, Inc. 
%   $Revision: 1.1.6.3 $  $Date: 2004/01/24 09:36:42 $

if ~isvector(p) || numel(p) == 0
    error('stats:quantile:BadProbs', ...
          'P must be a scalar or a non-empty vector.');
elseif any(p < 0 | p > 1)
    error('stats:quantile:BadProbs', ...
          'P must take values between 0.0 and 1.0');
end

if nargin < 3
    y = prctile(x,100.*p);
else
    y = prctile(x,100.*p,dim);
end

end