


Function to compute lower and upper confidence intervals on the spectrum
Usage: Serr=specerr(S,J,err,trialave,numsp)
Outputs: Serr (Serr(1,...) - lower confidence level, Serr(2,...) upper confidence level)
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
Outputs:
Serr - error estimates. Only for err(1)>=1. If err=[1 p] or [2 p] Serr(...,1) and Serr(...,2)
contain the lower and upper error bars with the specified method.

0001 function Serr=specerr(S,J,err,trialave,numsp) 0002 % Function to compute lower and upper confidence intervals on the spectrum 0003 % Usage: Serr=specerr(S,J,err,trialave,numsp) 0004 % Outputs: Serr (Serr(1,...) - lower confidence level, Serr(2,...) upper confidence level) 0005 % 0006 % Inputs: 0007 % S - spectrum 0008 % J - 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 % numsp - number of spikes in each channel. specify only when finite 0014 % size correction required 0015 % 0016 % Outputs: 0017 % Serr - error estimates. Only for err(1)>=1. If err=[1 p] or [2 p] Serr(...,1) and Serr(...,2) 0018 % contain the lower and upper error bars with the specified method. 0019 if nargin < 4; error('Need at least 4 input arguments'); end; 0020 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; 0021 [nf,K,C]=size(J); 0022 errchk=err(1); 0023 p=err(2); 0024 pp=1-p/2; 0025 qq=1-pp; 0026 0027 if trialave 0028 dim=K*C; 0029 C=1; 0030 dof=2*dim; 0031 if nargin==5; dof = fix(1/(1/dof + 1/(2*sum(numsp)))); end 0032 J=reshape(J,nf,dim); 0033 else 0034 dim=K; 0035 dof=2*dim*ones(1,C); 0036 for ch=1:C; 0037 if nargin==5; dof(ch) = fix(1/(1/dof + 1/(2*numsp(ch)))); end 0038 end; 0039 end; 0040 Serr=zeros(2,nf,C); 0041 if errchk==1; 0042 Qp=chi2inv(pp,dof); 0043 Qq=chi2inv(qq,dof); 0044 Serr(1,:,:)=dof(ones(nf,1),:).*S./Qp(ones(nf,1),:); 0045 Serr(2,:,:)=dof(ones(nf,1),:).*S./Qq(ones(nf,1),:); 0046 elseif errchk==2; 0047 tcrit=tinv(pp,dof); 0048 for k=1:dim; 0049 indices=setdiff(1:dim,k); 0050 Jjk=J(:,indices,:); % 1-drop projection 0051 eJjk=squeeze(sum(Jjk.*conj(Jjk),2)); 0052 Sjk(k,:,:)=eJjk/(dim-1); % 1-drop spectrum 0053 end; 0054 sigma=sqrt(dim-1)*squeeze(std(log(Sjk),1,1)); if C==1; sigma=sigma'; end; 0055 conf=tcrit(ones(nf,1),:).*sigma; 0056 conf=squeeze(conf); 0057 Serr(1,:,:)=S.*exp(-conf); Serr(2,:,:)=S.*exp(conf); 0058 end; 0059 Serr=squeeze(Serr);