# 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 (and of course, only for point
                process data)

     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:

- mtspectrumc Multi-taper spectrum - continuous process
- mtspectrumc_unequal_length_trials This routine computes the multi-taper spectrum for a given set of unequal length segments. It is
- mtspectrumsegc Multi-taper segmented spectrum for a univariate continuous process
- mtspectrumpb Multi-taper spectrum - binned point process
- mtspectrumsegpb Multi-taper segmented spectrum for a univariate binned point process
- mtspectrumpt Multi-taper spectrum - point process times
- mtspectrumsegpt Multi-taper segmented spectrum for a univariate binned point process

## 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 (and of course, only for point
    0015 %            process data)
    0016 %
    0017 % Outputs:
    0018 % Serr - error estimates. Only for err(1)>=1. If err=[1 p] or [2 p] Serr(...,1) and Serr(...,2)
    0019 % contain the lower and upper error bars with the specified method.
    0020 if nargin < 4; error('Need at least 4 input arguments'); end;
    0021 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;
    0022 [nf,K,C]=size(J);
    0023 errchk=err(1);
    0024 p=err(2);
    0025 pp=1-p/2;
    0026 qq=1-pp;
    0027
    0028 if trialave
    0029    dim=K*C;
    0030    C=1;
    0031    dof=2*dim;
    0032    if nargin==5; dof = fix(1/(1/dof + 1/(2*sum(numsp)))); end
    0033    J=reshape(J,nf,dim);
    0034 else
    0035    dim=K;
    0036    dof=2*dim*ones(1,C);
    0037    for ch=1:C;
    0038      if nargin==5; dof(ch) = fix(1/(1/dof + 1/(2*numsp(ch)))); end
    0039    end;
    0040 end;
    0041 Serr=zeros(2,nf,C);
    0042 if errchk==1;
    0043    Qp=chi2inv(pp,dof);
    0044    Qq=chi2inv(qq,dof);
    0045    Serr(1,:,:)=dof(ones(nf,1),:).*S./Qp(ones(nf,1),:);
    0046    Serr(2,:,:)=dof(ones(nf,1),:).*S./Qq(ones(nf,1),:);
    0047 elseif errchk==2;
    0048    tcrit=tinv(pp,dim-1);
    0049    for k=1:dim;
    0050        indices=setdiff(1:dim,k);
    0051        Jjk=J(:,indices,:); % 1-drop projection
    0052        eJjk=squeeze(sum(Jjk.*conj(Jjk),2));
    0053        Sjk(k,:,:)=eJjk/(dim-1); % 1-drop spectrum
    0054    end;
    0055    sigma=sqrt(dim-1)*squeeze(std(log(Sjk),1,1)); if C==1; sigma=sigma'; end;
    0056    conf=repmat(tcrit,nf,C).*sigma;
    0057    conf=squeeze(conf);
    0058    Serr(1,:,:)=S.*exp(-conf); Serr(2,:,:)=S.*exp(conf);
    0059 end;
    0060 Serr=squeeze(Serr);