Home > chronux_2_00 > spectral_analysis > helper > coherr.m

coherr

PURPOSE ^

Function to compute lower and upper confidence intervals on the coherency

SYNOPSIS ^

function [confC,phistd,Cerr]=coherr(C,J1,J2,err,trialave,numsp1,numsp2)

DESCRIPTION ^

 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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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;

Generated on Fri 15-Aug-2008 11:35:42 by m2html © 2003