Home > chronux_1_0 > helper > specerr.m

specerr

PURPOSE ^

Function to compute lower and upper confidence intervals on the spectrum

SYNOPSIS ^

function Serr=specerr(S,J,err,trialave,numsp)

DESCRIPTION ^

 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.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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);

Generated on Fri 09-Jun-2006 23:38:05 by m2html © 2003