Function to compute lower and upper confidence intervals on the coherency
given the tapered fourier transforms, errchk, trialave.
Usage: [confC,phistd,Cerr]=coherr(C,J1,J2,err,trialave,numsp1,numsp2)
Inputs:
C - coherence
J1,J2 - tapered fourier transforms
err - [errtype p] (errtype=1 - asymptotic estimates; errchk=2 - Jackknife estimates;
p - p value for error estimates)
trialave - 0: no averaging over trials/channels
1 : perform trial averaging
numsp1 - number of spikes for data1. supply only if finite size corrections are required
numsp2 - number of spikes for data2. supply only if finite size corrections are required
Outputs:
confC - confidence level for C - only for err(1)>=1
phistd - theoretical or jackknife standard deviation for phi for err(1)=1 and err(1)=2
respectively. returns zero if coherence is 1
Cerr - Jacknife error bars for C - only for err(1)=2
Jackknife uses the following transform of the coherence
z=sqrt(2*dim-2)atanh(C). Asymptotically (and for Gaussian data) var(z)=1.0001 function [confC,phistd,Cerr]=coherr(C,J1,J2,err,trialave,numsp1,numsp2) 0002 % Function to compute lower and upper confidence intervals on the coherency 0003 % given the tapered fourier transforms, errchk, trialave. 0004 % 0005 % Usage: [confC,phistd,Cerr]=coherr(C,J1,J2,err,trialave,numsp1,numsp2) 0006 % Inputs: 0007 % C - coherence 0008 % J1,J2 - tapered fourier transforms 0009 % err - [errtype p] (errtype=1 - asymptotic estimates; errchk=2 - Jackknife estimates; 0010 % p - p value for error estimates) 0011 % trialave - 0: no averaging over trials/channels 0012 % 1 : perform trial averaging 0013 % numsp1 - number of spikes for data1. supply only if finite size corrections are required 0014 % numsp2 - number of spikes for data2. supply only if finite size corrections are required 0015 % 0016 % Outputs: 0017 % confC - confidence level for C - only for err(1)>=1 0018 % phistd - theoretical or jackknife standard deviation for phi for err(1)=1 and err(1)=2 0019 % respectively. returns zero if coherence is 1 0020 % Cerr - Jacknife error bars for C - only for err(1)=2 0021 % Jackknife uses the following transform of the coherence 0022 % z=sqrt(2*dim-2)atanh(C). Asymptotically (and for Gaussian data) var(z)=1. 0023 0024 if nargin < 5; error('Need at least 5 input arguments'); end; 0025 if err(1)==0; error('Need err=[1 p] or [2 p] for error bar calculation'); end; 0026 if nargout==4 && err(1)==1; error('Cerr contains Jackknife errors: only computed when err(1) is 2'); end; 0027 [nf,K,Ch]=size(J1); 0028 errchk=err(1); 0029 p=err(2); 0030 pp=1-p/2; 0031 % 0032 % Find the number of degrees of freedom 0033 % 0034 if trialave; 0035 dim=K*Ch; 0036 dof=2*dim; 0037 dof1=dof; 0038 dof2=dof; 0039 Ch=1; 0040 if nargin>=6 && ~isempty(numsp1) 0041 totspikes1=sum(numsp1); 0042 dof1=fix(2*totspikes1*dof/(2*totspikes1+dof)); 0043 end 0044 if nargin==7 && ~isempty(numsp2); 0045 totspikes2=sum(numsp2); 0046 dof2=fix(2*totspikes2*dof/(2*totspikes2+dof)); 0047 end; 0048 dof=min(dof1,dof2); 0049 J1=reshape(J1,nf,dim); 0050 J2=reshape(J2,nf,dim); 0051 else 0052 dim=K; 0053 dof=2*dim; 0054 dof1=dof; 0055 dof2=dof; 0056 for ch=1:Ch; 0057 if nargin>=6 && ~isempty(numsp1); 0058 totspikes1=numsp1(ch); 0059 dof1=fix(2*totspikes1*dof/(2*totspikes1+dof)); 0060 end; 0061 if nargin==7 && ~isempty(numsp2); 0062 totspikes2=numsp2(ch); 0063 dof2=fix(2*totspikes2*dof/(2*totspikes2+dof)); 0064 end; 0065 dof(ch)=min(dof1,dof2); 0066 end; 0067 end; 0068 % 0069 % variance of the phase 0070 % 0071 % 0072 % Old code is the next few lines - new code is in the if statement below 0073 % beginning line 87 0074 % 0075 % if isempty(find((C-1).^2 < 10^-16)); 0076 % phierr = sqrt((2./dof(ones(nf,1),:)).*(1./(C.^2) - 1)); 0077 % else 0078 % phierr = zeros(nf,Ch); 0079 % end 0080 0081 % 0082 % theoretical, asymptotic confidence level 0083 % 0084 if dof <= 2 0085 confC = 1; 0086 else 0087 df = 1./((dof/2)-1); 0088 confC = sqrt(1 - p.^df); 0089 end; 0090 % 0091 % Phase standard deviation (theoretical and jackknife) and jackknife 0092 % confidence intervals for C 0093 % 0094 if errchk==1; 0095 totnum=nf*Ch; 0096 phistd=zeros(totnum,1); 0097 CC=reshape(C,[totnum,1]); 0098 indx=find(abs(CC-1)>=1.e-16); 0099 dof=repmat(dof,[nf,1]); 0100 dof=reshape(dof,[totnum 1]); 0101 phistd(indx)= sqrt((2./dof(indx).*(1./(C(indx).^2) - 1))); 0102 phistd=reshape(phistd,[nf Ch]); 0103 elseif errchk==2; 0104 tcrit=tinv(pp,dof-1); 0105 for k=1:dim; % dim is the number of 'independent' estimates 0106 indxk=setdiff(1:dim,k); 0107 J1k=J1(:,indxk,:); 0108 J2k=J2(:,indxk,:); 0109 eJ1k=squeeze(sum(J1k.*conj(J1k),2)); 0110 eJ2k=squeeze(sum(J2k.*conj(J2k),2)); 0111 eJ12k=squeeze(sum(conj(J1k).*J2k,2)); 0112 Cxyk=eJ12k./sqrt(eJ1k.*eJ2k); 0113 absCxyk=abs(Cxyk); 0114 atanhCxyk(k,:,:)=sqrt(2*dim-2)*atanh(absCxyk); % 1-drop estimate of z 0115 phasefactorxyk(k,:,:)=Cxyk./absCxyk; 0116 % indxk=setdiff(1:dim,k); 0117 % J1jk=J1(:,indxk,:); 0118 % J2jk=J2(:,indxk,:); 0119 % eJ1jk=squeeze(sum(J1jk.*conj(J1jk),2)); 0120 % eJ2jk=squeeze(sum(J2jk.*conj(J2jk),2)); 0121 % eJ12jk=squeeze(sum(conj(J1jk).*J2jk,2)); 0122 % atanhCxyjk(k,:,:)=sqrt(2*dim-2)*atanh(abs(eJ12jk)./sqrt(eJ1jk.*eJ2jk)); 0123 end; 0124 atanhC=sqrt(2*dim-2)*atanh(C); % z 0125 sigma12=sqrt(dim-1)*squeeze(std(atanhCxyk,1,1)); % Jackknife estimate std(z)=sqrt(dim-1)*std of 1-drop estimates 0126 % sigma12=sqrt(dim-1)*squeeze(std(atanhCxyjk,1,1)); 0127 if Ch==1; sigma12=sigma12'; end; 0128 Cu=atanhC+tcrit(ones(nf,1),:).*sigma12; 0129 Cl=atanhC-tcrit(ones(nf,1),:).*sigma12; 0130 %Cerr(1,:,:) = tanh(Cl/sqrt(2*dim-2)); 0131 Cerr(1,:,:) = max(tanh(Cl/sqrt(2*dim-2)),0); % This ensures that the lower confidence band remains positive 0132 Cerr(2,:,:) = tanh(Cu/sqrt(2*dim-2)); 0133 %phistd=(2*dim-2)*(1-abs(squeeze(mean(phasefactorxyk)))); 0134 phistd = sqrt( (2*dim-2)*(1-abs(squeeze(mean(phasefactorxyk)))) ); 0135 if trialave; phistd=phistd'; end; 0136 end 0137 % ncrit=norminv(pp); 0138 % phierr=zeros([2 size(phistd)]); 0139 % phierr(1,:,:)=phi-ncrit*phistd; phierr(2,:,:)=phi+ncrit*phistd;