uispecerr

 Function to compute lower and upper confidence intervals on the spectrum 
 Usage: [Serrp,Serrj]=uispecerr(S,J,err,trialave,numsp)
 Outputs: Serrp (Serrp(1,...) - lower confidence level, Serrp(2,...) upper
 confidence level), similarly for Serrj

 Inputs:
 S - spectrum
 J - 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
 numsp    - number of spikes in each channel. specify only when finite
            size correction required (and of course, only for point
            process data)

 Outputs:
 Serrp - population error estimates. Only for err(1)>=1. 
 Serrj - jackknife error estimates. Computed only for err(1)>=2, otherwise
 set to zero

0001 function [Serrp,Serrj]=uispecerr(S,J,err,trialave,numsp)
0002 % Function to compute lower and upper confidence intervals on the spectrum
0003 % Usage: [Serrp,Serrj]=uispecerr(S,J,err,trialave,numsp)
0004 % Outputs: Serrp (Serrp(1,...) - lower confidence level, Serrp(2,...) upper
0005 % confidence level), similarly for Serrj
0006 %
0007 % Inputs:
0008 % S - spectrum
0009 % J - tapered fourier transforms
0010 % err - [errtype p] (errtype=1 - asymptotic estimates; errchk=2 - Jackknife estimates;
0011 %                   p - p value for error estimates)
0012 % trialave - 0: no averaging over trials/channels
0013 %            1 : perform trial averaging
0014 % numsp    - number of spikes in each channel. specify only when finite
0015 %            size correction required (and of course, only for point
0016 %            process data)
0017 %
0018 % Outputs:
0019 % Serrp - population error estimates. Only for err(1)>=1.
0020 % Serrj - jackknife error estimates. Computed only for err(1)>=2, otherwise
0021 % set to zero
0022 if nargin < 4; error('Need at least 4 input arguments'); end;
0023 if err(1)==0; error('Need err=[1 p] or [2 p] for error bar calculation. Make sure you are not asking for the output of Serr'); end;
0024 [nf,K,C]=size(J);
0025 errchk=err(1);
0026 p=err(2);
0027 pp=1-p/2;
0028 qq=1-pp;
0029 
0030 if trialave
0031    dim=K*C;
0032    C=1;
0033    dof=2*dim;
0034    if nargin==5; dof = fix(1/(1/dof + 1/(2*sum(numsp)))); end
0035    J=reshape(J,nf,dim);
0036 else
0037    dim=K;
0038    dof=2*dim*ones(1,C);
0039    for ch=1:C;
0040      if nargin==5; dof(ch) = fix(1/(1/dof + 1/(2*numsp(ch)))); end 
0041    end;
0042 end;
0043 Serrp=zeros(2,nf,C);
0044 Serrj=zeros(2,nf,C);
0045 if errchk>1=;
0046    Qp=chi2inv(pp,dof);
0047    Qq=chi2inv(qq,dof);
0048    Serrp(1,:,:)=dof(ones(nf,1),:).*S./Qp(ones(nf,1),:);
0049    Serrp(2,:,:)=dof(ones(nf,1),:).*S./Qq(ones(nf,1),:);
0050 elseif errchk>=2;
0051    tcrit=tinv(pp,dim-1);
0052    for k=1:dim;
0053        indices=setdiff(1:dim,k);
0054        Jjk=J(:,indices,:); % 1-drop projection
0055        eJjk=squeeze(sum(Jjk.*conj(Jjk),2));
0056        Sjk(k,:,:)=eJjk/(dim-1); % 1-drop spectrum
0057    end;
0058    sigma=sqrt(dim-1)*squeeze(std(log(Sjk),1,1)); if C==1; sigma=sigma'; end; 
0059    conf=repmat(tcrit,nf,C).*sigma;
0060    conf=squeeze(conf); 
0061    Serrj(1,:,:)=S.*exp(-conf); Serrj(2,:,:)=S.*exp(conf);
0062 end;
0063 Serrp=squeeze(Serrp);
0064 Serrj=squeeze(Serrj);