# 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
     Jackknife uses the following transform of the coherence
     z=sqrt(2*dim-2)atanh(C). Asymptotically (and for Gaussian data) var(z)=1.

## CROSS-REFERENCE INFORMATION

This function calls:

This function is called by:

- coherencyc Multi-taper coherency,cross-spectrum and individual spectra - continuous process
- coherencyc_unequal_length_trials This routine computes the average multi-taper coherence for a given set of unequal length segments. It is
- cohmathelper Helper function called by coherency matrix computations.
- coherencycpb Multi-taper coherency,cross-spectrum and individual spectra - continuous and binned point process data
- coherencycpt Multi-taper coherency,cross-spectrum and individual spectra -continuous data and point process as times
- coherencypb Multi-taper coherency,cross-spectrum and individual spectra - binned point process
- coherencypt Multi-taper coherency - point process times

## 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 % 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;