Description of TestSkewness

0001 function [h,p,stats] = TestSkewness(control,repeat,test,varargin)
0002 
0003 %TestSkewness - Test if firing field skewness changes between two conditions.
0004 %
0005 %  Current implementation only tests 1D environments.
0006 %
0007 %  USAGE
0008 %
0009 %    [h,p,cr,ca] = TestSkewness(control,repeat,test,<options>)
0010 %
0011 %    control        firing fields in control condition (MxN: M fields, N bins)
0012 %    repeat         optional firing fields in repeated control condition
0013 %    test           firing fields in test condition
0014 %    <options>      optional list of property-value pairs (see table below)
0015 %
0016 %    =========================================================================
0017 %     Properties    Values
0018 %    -------------------------------------------------------------------------
0019 %     'alpha'       significance level (default = 0.05)
0020 %     'iterations'  number of iterations for bootstrap (default = 150)
0021 %     'show'        plot results (default = 'off')
0022 %     'style'       'hist' for histograms (default) or 'curves' for smooth
0023 %                   curves
0024 %     'bins'        [m M n], lower and upper bounds, and number of bins,
0025 %                   respectively (default [min max 100])
0026 %    =========================================================================
0027 %
0028 %  OUTPUT
0029 %
0030 %    h                  1 if H0 (random remapping) can be rejected
0031 %    p                  p-value of bootstrap test
0032 %    stats.control.m       median skewness (control)
0033 %    stats.control.s       standard error of the median skewness (control)
0034 %    stats.repeat.m       median skewness (repeat)
0035 %    stats.repeat.s       standard error of the median skewness (repeat)
0036 %    stats.test.m          median skewness (test)
0037 %    stats.test.s          standard error of the median skewness (test)
0038 
0039 % Copyright (C) 2013 by Michaël Zugaro and Anne Cei
0040 %
0041 % This program is free software; you can redistribute it and/or modify
0042 % it under the terms of the GNU General Public License as published by
0043 % the Free Software Foundation; either version 3 of the License, or
0044 % (at your option) any later version.
0045 
0046 % Defaults
0047 show = 'off';
0048 nBootstrap = 150;
0049 alpha = 0.05;
0050 style = 'hist';
0051 bins = [];
0052 
0053 % Check number of parameters
0054 if nargin < 3 | mod(length(varargin),2) ~= 0,
0055   error('Incorrect number of parameters (type ''help <a href="matlab:help TestSkewness">TestSkewness</a>'' for details).');
0056 end
0057 
0058 % Check parameter sizes
0059 if ~isdmatrix(control) || ~isdmatrix(repeat) || ~isdmatrix(test),
0060     error('All firing fields should be MxN matrices (type ''help <a href="matlab:help TestSkewness">TestSkewness</a>'' for details).');
0061 end
0062 
0063 % Parse parameter list
0064 for i = 1:2:length(varargin),
0065     if ~ischar(varargin{i}),
0066         error(['Parameter ' num2str(i+2) ' is not a property (type ''help <a href="matlab:help TestSkewness">TestSkewness</a>'' for details).']);
0067     end
0068     switch(lower(varargin{i})),
0069         case 'alpha',
0070             alpha = varargin{i+1};
0071             if ~isdscalar(alpha,'>=0','<=1'),
0072                 error('Incorrect value for property ''alpha'' (type ''help <a href="matlab:help TestSkewness">TestSkewness</a>'' for details).');
0073             end
0074         case 'iterations',
0075             nBootstrap = varargin{i+1};
0076             if ~isiscalar(nBootstrap,'>0'),
0077                 error('Incorrect value for property ''iterations'' (type ''help <a href="matlab:help TestSkewness">TestSkewness</a>'' for details).');
0078             end
0079         case 'show',
0080             show = varargin{i+1};
0081             if ~isastring(show,'on','off'),
0082                 error('Incorrect value for property ''show'' (type ''help <a href="matlab:help TestSkewness">TestSkewness</a>'' for details).');
0083             end
0084         case 'style',
0085             style = varargin{i+1};
0086             if ~isastring(style,'hist','curves'),
0087                 error('Incorrect value for property ''style'' (type ''help <a href="matlab:help TestSkewness">TestSkewness</a>'' for details).');
0088             end
0089         case 'bins',
0090             bins = varargin{i+1};
0091             if ~isdvector(bins,'#3') | bins(1) > bins(2) | ~isiscalar(bins(3),'>0'),
0092                 error('Incorrect value for property ''bins'' (type ''help <a href="matlab:help TestSkewness">TestSkewness</a>'' for details).');
0093             end
0094         otherwise,
0095             error(['Unknown property ''' num2str(varargin{i}) ''' (type ''help <a href="matlab:help TestSkewness">TestSkewness</a>'' for details).']);
0096     end
0097 end
0098 
0099 % Smooth and bins
0100 smooth = 10;
0101 if ~isempty(bins),
0102     xMin = bins(1);
0103     xMax = bins(2);
0104     nBins = bins(3);
0105 else
0106     % Default values
0107     xMin = -1;
0108     xMax = 1;
0109     % (number of bins differ for smooth curves vs histograms)
0110     if strcmp(style,'hist'),
0111         nBins = 20;
0112     else
0113         nBins = 200;
0114     end
0115 end
0116 
0117 % Discard NaNs
0118 bad = any(isnan(control),2);
0119 control(bad,:) = [];
0120 bad = any(isnan(repeat),2);
0121 repeat(bad,:) = [];
0122 bad = any(isnan(test),2);
0123 test(bad,:) = [];
0124 
0125 % Matrix sizes
0126 [mc,nc] = size(control);
0127 [mr,nr] = size(repeat);
0128 [mt,nt] = size(test);
0129 
0130 % Normalize: transform firing rates into firing probabilities
0131 control = control ./ repmat(sum(control,2),1,nc);
0132 repeat = repeat ./ repmat(sum(repeat,2),1,nr);
0133 test = test ./ repmat(sum(test,2),1,nt);
0134 
0135 % Compute skewness
0136 skewnessControl = UnbiasedSkewness(control);
0137 skewnessRepeat = UnbiasedSkewness(repeat);
0138 skewnessTest = UnbiasedSkewness(test);
0139 
0140 % Test distribution differences and compute medians
0141 [hc,pc,ksc] = kstest2(skewnessControl,skewnessRepeat,alpha);
0142 [h,p,ks] = kstest2(skewnessControl,skewnessTest,alpha);
0143 stats.control.m = median(skewnessControl);
0144 stats.repeat.m = median(skewnessRepeat);
0145 stats.test.m = median(skewnessTest);
0146 stats.control.s = semedian(skewnessControl);
0147 stats.repeat.s = semedian(skewnessRepeat);
0148 stats.test.s = semedian(skewnessTest);
0149 
0150 % Plot results
0151 if strcmp(show,'on'),
0152 
0153     figure;hold on;
0154     % Test results
0155     sc = ['KS(' int2str(mc) ',' int2str(mr) ')=' sprintf('%.3f',ksc) ', p = ' num2str(pc)];
0156     st = ['KS(' int2str(mt) ',' int2str(mc) ')=' sprintf('%.3f',ks) ', p = ' num2str(p)];
0157 %      if hc,
0158 %          sc = ['p < ' num2str(alpha)];
0159 %      else
0160 %          sc = 'NS';
0161 %      end
0162 %      if h,
0163 %          st = ['p < ' num2str(alpha)];
0164 %      else
0165 %          st = 'NS';
0166 %      end
0167 
0168     % Histograms
0169     x = linspace(xMin,xMax,nBins);
0170     dx = x(2)-x(1);
0171     hControl = hist(skewnessControl,x)/mc;
0172     hRepeat = hist(skewnessRepeat,x)/mr;
0173     hTest = hist(skewnessTest,x)/mt;
0174     % Plot control and repeat
0175     subplot(2,1,1);hold on;
0176     if strcmp(style,'hist'),
0177         bar(x,hControl,0.5,'b');
0178         bar(x+dx/2,hRepeat,0.5,'k');
0179     else
0180         plot(x,Smooth(hControl,smooth),'b','linewidth',2);
0181         plot(x+dx/2,Smooth(hRepeat,smooth),'k','linewidth',2);
0182     end
0183     PlotHVLines(stats.control.m,'v','b');
0184     PlotHVLines(stats.repeat.m,'v','k');
0185     PlotHVLines(0,'v','k--');
0186     legend(['control ' mean2str(stats.control.m,stats.control.s,mc)],['repeat ' mean2str(stats.repeat.m,stats.repeat.s,mr)],'location','northwest');
0187     title(['Control vs Repeat (KS test: ' sc ')']);
0188     xlabel('Skewness');
0189     ylabel('Probability');
0190     % Plot control and test
0191     subplot(2,1,2);hold on;
0192     if strcmp(style,'hist'),
0193         bar(x,hControl,0.5,'b');
0194         bar(x+dx/2,hTest,0.5,'r');
0195     else
0196         plot(x,Smooth(hControl,smooth),'b','linewidth',2);
0197         plot(x+dx/2,Smooth(hTest,smooth),'r','linewidth',2);
0198     end
0199     PlotHVLines(stats.control.m,'v','b');
0200     PlotHVLines(stats.test.m,'v','r');
0201     PlotHVLines(0,'v','k--');
0202     legend(['control ' mean2str(stats.control.m,stats.control.s,mc)],['repeat ' mean2str(stats.test.m,stats.test.s,mt)],'location','northwest');
0203     title(['Control vs Test (KS test: ' st ')']);
0204     xlabel('Skewness');
0205     ylabel('Probability');
0206 
0207 end
0208 
0209 function s = UnbiasedSkewness(distributions)
0210 
0211 [m,n] = size(distributions);
0212 x = repmat(1:n,m,1);
0213 
0214 % Mean for each distribution
0215 mu = sum(x.*distributions,2);
0216 % Centered x
0217 x0 = x-repmat(mu,1,n);
0218 % Third central moment
0219 m3 = sum( x0.^3 .* distributions,2 );
0220 % Second central moment
0221 s2 = sum( x0.^2 .* distributions,2 );
0222 % Biased skewness
0223 s = m3./s2.^1.5;
0224 % Unbiased skewness
0225 s = sqrt(n*(n-1))/(n-2) * s;
TestSkewness

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
