


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

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 0022 if nargin < 5; error('Need at least 5 input arguments'); end; 0023 if err(1)==0; error('Need err=[1 p] or [2 p] for error bar calculation'); end; 0024 if nargout==4 && err(1)==1; error('Cerr contains Jackknife errors: only computed when err(1) is 2'); end; 0025 [nf,K,Ch]=size(J1); 0026 errchk=err(1); 0027 p=err(2); 0028 pp=1-p/2; 0029 % 0030 % Find the number of degrees of freedom 0031 % 0032 if trialave; 0033 dim=K*Ch; 0034 dof=2*dim; 0035 dof1=dof; 0036 dof2=dof; 0037 Ch=1; 0038 if nargin>=6 && ~isempty(numsp1) 0039 totspikes1=sum(numsp1); 0040 dof1=fix(2*totspikes1*dof/(2*totspikes1+dof)); 0041 end 0042 if nargin==7 && ~isempty(numsp2); 0043 totspikes2=sum(numsp2); 0044 dof2=fix(2*totspikes2*dof/(2*totspikes2+dof)); 0045 end; 0046 dof=min(dof1,dof2); 0047 J1=reshape(J1,nf,dim); 0048 J2=reshape(J2,nf,dim); 0049 else 0050 dim=K; 0051 dof=2*dim; 0052 dof1=dof; 0053 dof2=dof; 0054 for ch=1:Ch; 0055 if nargin>=6 && ~isempty(numsp1); 0056 totspikes1=numsp1(ch); 0057 dof1=fix(2*totspikes1*dof/(2*totspikes1+dof)); 0058 end; 0059 if nargin==7 && ~isempty(numsp2); 0060 totspikes2=numsp2(ch); 0061 dof2=fix(2*totspikes2*dof/(2*totspikes2+dof)); 0062 end; 0063 dof(ch)=min(dof1,dof2); 0064 end; 0065 end; 0066 % 0067 % variance of the phase 0068 % 0069 % 0070 % Old code is the next few lines - new code is in the if statement below 0071 % beginning line 87 0072 % 0073 % if isempty(find((C-1).^2 < 10^-16)); 0074 % phierr = sqrt((2./dof(ones(nf,1),:)).*(1./(C.^2) - 1)); 0075 % else 0076 % phierr = zeros(nf,Ch); 0077 % end 0078 0079 % 0080 % theoretical, asymptotic confidence level 0081 % 0082 if dof <= 2 0083 confC = 1; 0084 else 0085 df = 1./((dof/2)-1); 0086 confC = sqrt(1 - p.^df); 0087 end; 0088 % 0089 % Phase standard deviation (theoretical and jackknife) and jackknife 0090 % confidence intervals for C 0091 % 0092 if errchk==1; 0093 totnum=nf*Ch; 0094 phistd=zeros(totnum,1); 0095 CC=reshape(C,[totnum,1]); 0096 indx=find(abs(CC-1)>=1.e-16); 0097 dof=repmat(dof,[nf,1]); 0098 dof=reshape(dof,[totnum 1]); 0099 phistd(indx)= sqrt((2./dof(indx).*(1./(C(indx).^2) - 1))); 0100 phistd=reshape(phistd,[nf Ch]); 0101 elseif errchk==2; 0102 tcrit=tinv(pp,dof-1); 0103 for k=1:dim; 0104 indxk=setdiff(1:dim,k); 0105 J1k=J1(:,indxk,:); 0106 J2k=J2(:,indxk,:); 0107 eJ1k=squeeze(sum(J1k.*conj(J1k),2)); 0108 eJ2k=squeeze(sum(J2k.*conj(J2k),2)); 0109 eJ12k=squeeze(sum(conj(J1k).*J2k,2)); 0110 Cxyk=eJ12k./sqrt(eJ1k.*eJ2k); 0111 absCxyk=abs(Cxyk); 0112 atanhCxyk(k,:,:)=sqrt(2*dim-2)*atanh(absCxyk); 0113 phasefactorxyk(k,:,:)=Cxyk./absCxyk; 0114 % indxk=setdiff(1:dim,k); 0115 % J1jk=J1(:,indxk,:); 0116 % J2jk=J2(:,indxk,:); 0117 % eJ1jk=squeeze(sum(J1jk.*conj(J1jk),2)); 0118 % eJ2jk=squeeze(sum(J2jk.*conj(J2jk),2)); 0119 % eJ12jk=squeeze(sum(conj(J1jk).*J2jk,2)); 0120 % atanhCxyjk(k,:,:)=sqrt(2*dim-2)*atanh(abs(eJ12jk)./sqrt(eJ1jk.*eJ2jk)); 0121 end; 0122 atanhC=sqrt(2*dim-2)*atanh(C); 0123 sigma12=sqrt(dim-1)*squeeze(std(atanhCxyk,1,1)); 0124 % sigma12=sqrt(dim-1)*squeeze(std(atanhCxyjk,1,1)); 0125 if Ch==1; sigma12=sigma12'; end; 0126 Cu=atanhC+tcrit(ones(nf,1),:).*sigma12; 0127 Cl=atanhC-tcrit(ones(nf,1),:).*sigma12; 0128 Cerr(1,:,:) = tanh(Cl/sqrt(2*dim-2)); 0129 Cerr(2,:,:) = tanh(Cu/sqrt(2*dim-2)); 0130 phistd=(2*dim-2)*(1-abs(squeeze(mean(phasefactorxyk)))); 0131 if trialave; phistd=phistd'; end; 0132 end 0133 % ncrit=norminv(pp); 0134 % phierr=zeros([2 size(phistd)]); 0135 % phierr(1,:,:)=phi-ncrit*phistd; phierr(2,:,:)=phi+ncrit*phistd;