Home > chronux_1_0 > statistical_tests > two_group_test_spectrum.m

two_group_test_spectrum

PURPOSE ^

function [dz,vdz,Adz]=two_group_test_spectrum(J1,J2,p)

SYNOPSIS ^

function [dz,vdz,Adz]=two_group_test_spectrum(J1,J2,p)

DESCRIPTION ^

 function [dz,vdz,Adz]=two_group_test_spectrum(J1,J2,p)
 Test the null hypothesis (H0) that data sets J1, J2 in 
 two conditions c1,c2 have equal population spectrum

 Usage:
 [dz,vdz,Adz]=two_sample_test_spectrum(J1,J2,p)

 Inputs:
 J1   tapered fourier transform in condition 1
 J2   tapered fourier transform in condition 2
 p      p value for test (default: 0.05)


 Dimensions: J1: frequencies x number of samples in condition 1
             J2: frequencies x number of samples in condition 2
              number of samples = number of trials x number of tapers
 Outputs:
 dz    test statistic (will be distributed as N(0,1) under H0
 vdz   Arvesen estimate of the variance of dz
 Adz   1/0 for accept/reject null hypothesis of equal population
       coherences based dz ~ N(0,1)
 
 
 Note: all outputs are functions of frequency

 References: Arvesen, Jackkknifing U-statistics, Annals of Mathematical
 Statisitics, vol 40, no. 6, pg 2076-2100 (1969)

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [dz,vdz,Adz]=two_group_test_spectrum(J1,J2,p)
0002 % function [dz,vdz,Adz]=two_group_test_spectrum(J1,J2,p)
0003 % Test the null hypothesis (H0) that data sets J1, J2 in
0004 % two conditions c1,c2 have equal population spectrum
0005 %
0006 % Usage:
0007 % [dz,vdz,Adz]=two_sample_test_spectrum(J1,J2,p)
0008 %
0009 % Inputs:
0010 % J1   tapered fourier transform in condition 1
0011 % J2   tapered fourier transform in condition 2
0012 % p      p value for test (default: 0.05)
0013 %
0014 %
0015 % Dimensions: J1: frequencies x number of samples in condition 1
0016 %             J2: frequencies x number of samples in condition 2
0017 %              number of samples = number of trials x number of tapers
0018 % Outputs:
0019 % dz    test statistic (will be distributed as N(0,1) under H0
0020 % vdz   Arvesen estimate of the variance of dz
0021 % Adz   1/0 for accept/reject null hypothesis of equal population
0022 %       coherences based dz ~ N(0,1)
0023 %
0024 %
0025 % Note: all outputs are functions of frequency
0026 %
0027 % References: Arvesen, Jackkknifing U-statistics, Annals of Mathematical
0028 % Statisitics, vol 40, no. 6, pg 2076-2100 (1969)
0029 
0030 if nargin < 2; error('Need four sets of Fourier transforms'); end;
0031 %
0032 % Test for matching dimensionalities
0033 %
0034 m1=size(J1,2); % number of samples, condition 1
0035 m2=size(J2,2); % number of samples, condition 2
0036 dof1=m1; % degrees of freedom, condition 1
0037 dof2=m2; % degrees of freedom, condition 2
0038 
0039 if nargin < 3; p=0.05; end; % set the default p value
0040 
0041 %
0042 % Compute the individual condition spectra, coherences
0043 %
0044 S1=conj(J1).*J1; % spectrum, condition 1
0045 S2=conj(J2).*J2; % spectrum, condition 2
0046 
0047 Sm1=squeeze(mean(S1,2)); % mean spectrum, condition 1
0048 Sm2=squeeze(mean(S2,2)); % mean spectrum, condition 2
0049 %
0050 % Compute the statistic dz, and the probability of observing the value dz
0051 % given an N(0,1) distribution i.e. under the null hypothesis
0052 %
0053 bias1=psi(dof1)-log(dof1); bias2=psi(dof2)-log(dof2); % bias from Thomson & Chave
0054 var1=psi(1,dof1); var2=psi(1,dof2); % variance from Thomson & Chave
0055 z1=log(Sm1)-bias1; % Bias-corrected Fisher z, condition 1
0056 z2=log(Sm2)-bias2; % Bias-corrected Fisher z, condition 2
0057 dz=(z1-z2)/sqrt(var1+var2); % z statistic
0058 pdz=normpdf(dz,0,1); % probability of observing value dz
0059 %
0060 % The remaining portion of the program computes Jackknife estimates of the mean (mdz) and variance (vdz) of dz
0061 %
0062 samples1=[1:m1];
0063 samples2=[1:m2];
0064 %
0065 % Leave one out of one sample
0066 %
0067 bias11=psi(dof1-1)-log(dof1-1); var11=psi(1,dof1-1);
0068 for i=1:m1;
0069     ikeep=setdiff(samples1,i); % all samples except i
0070     Sm1=squeeze(mean(S1(:,ikeep),2)); % 1 drop mean spectrum, data 1, condition 1
0071     z1i(:,i)=log(Sm1)-bias11; % 1 drop, bias-corrected Fisher z, condition 1
0072     dz1i(:,i)=(z1i(:,i)-z2)/sqrt(var11+var2); % 1 drop, z statistic, condition 1
0073     ps1(:,i)=m1*dz-(m1-1)*dz1i(:,i);
0074 end; 
0075 ps1m=mean(ps1,2);
0076 bias21=psi(dof2-1)-log(dof2-1); var21=psi(1,dof2-1);
0077 for j=1:m2;
0078     jkeep=setdiff(samples2,j); % all samples except j
0079     Sm2=squeeze(mean(S2(:,jkeep),2)); % 1 drop mean spectrum, data 2, condition 2
0080     z2j(:,j)=log(Sm2)-bias21; % 1 drop, bias-corrected Fisher z, condition 2
0081     dz2j(:,j)=(z1-z2j(:,j))/sqrt(var1+var21); % 1 drop, z statistic, condition 2
0082     ps2(:,j)=m2*dz-(m2-1)*dz2j(:,j);
0083 end;
0084 %
0085 % Leave one out, both samples
0086 % and pseudo values
0087 % for i=1:m1;
0088 %     for j=1:m2;
0089 %         dzij(:,i,j)=(z1i(:,i)-z2j(:,j))/sqrt(var11+var21);
0090 %         dzpseudoval(:,i,j)=m1*m2*dz-(m1-1)*m2*dz1i(:,i)-m1*(m2-1)*dz2j(:,j)+(m1-1)*(m2-1)*dzij(:,i,j);
0091 %     end;
0092 % end;
0093 %
0094 % Jackknife mean and variance
0095 %
0096 % dzah=sum(sum(dzpseudoval,3),2)/(m1*m2);
0097 ps2m=mean(ps2,2);
0098 % dzar=(sum(ps1,2)+sum(ps2,2))/(m1+m2);
0099 vdz=sum((ps1-ps1m(:,ones(1,m1))).*(ps1-ps1m(:,ones(1,m1))),2)/(m1*(m1-1))+sum((ps2-ps2m(:,ones(1,m2))).*(ps2-ps2m(:,ones(1,m2))),2)/(m2*(m2-1));
0100 % vdzah=sum(sum((dzpseudoval-dzah(:,ones(1,m1),ones(1,m2))).*(dzpseudoval-dzah(:,ones(1,m1),ones(1,m2))),3),2)/(m1*m2);
0101 %
0102 % Test whether H0 is accepted at the specified p value
0103 %
0104 Adz=zeros(size(dz));
0105 x=norminv([p/2 1-p/2],0,1);
0106 indx=find(dz>=x(1) & dz<=x(2)); 
0107 Adz(indx)=1;
0108 
0109 % Adzar=zeros(size(dzar));
0110 % indx=find(dzar>=x(1) & dzar<=x(2));
0111 % Adzar(indx)=1;
0112 %
0113 % Adzah=zeros(size(dzah));
0114 % indx=find(dzah>=x(1) & dzah<=x(2));
0115 % Adzah(indx)=1;

Generated on Fri 09-Jun-2006 23:38:05 by m2html © 2003