function [yes, pct_in] = Diagnostic3(resid)
%DIAGNOSTIC3: GIDE digital image deblurring
% function [yes, pct_in] = Diagnostic3(resid)
% 
% This function determines whether "Diagnostic 3" is satisfied,
% testing whether the input vector has elements that are samples 
% from a normal distribution with mean 0 and standard deviation 1.
%
% It also graphs the cumulative periodogram of the vector resid.
%
% Input:
%
%  resid - vector of length m
%
% Diagnostic 3: The cumulative periodogram of resid should
% be close to a straight line, with at most 5% of its points 
% lying outside its 95% confidence interval. 
%
% Output: 
%
%   yes    - 1 if diagnostic is satisfied, 0 otherwise
%
%   pct_in - percentage of points in the cumulative periodogram
%            that lie inside the 95% confidence interval for
%            a normally distributed residual
%
% Based on 
% B.W. Rust and D.P. O'Leary, Residual periodograms for choosing 
% regularization parameters for ill-posed problems, 
% Inverse Problmes, 24 (2008) 034005 
%
% Brianna Cash and Dianne P. O'Leary 05/2013, 
% following a Fortran program by Bert Rust
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

n = length(resid);

%  Compute the periodogram, the absolute value of the first 
%  half of the entries in the FFT of resid.

fresid = fft(resid);
fresid = fresid .* conj(fresid)*2/n; %conj finds the the complex conjugate
mhalf = ceil((n-1)/2);
fresid = fresid(1:mhalf);

%  Compute the cumulative periodogram and its length.

cumper = cumsum(fresid);     %cumulative sum
cumper = cumper/cumper(mhalf); 

% Compute test = number of entries outside the confidence interval.
%  1.36/sqrt(mhalf-1) gives 95% point for mhalf>31; 
%  1.63/sqrt(mhalf-1) gives 99%
%  See WA Fuller, Introduction to Statistical Time Series, Wiley 1996, p.287.

ss = 1.36/sqrt(mhalf-1);
upperbnd = (2*[1:mhalf]/(2*mhalf) + ss)';
lowerbnd = (2*[1:mhalf]/(2*mhalf) - ss)';

test = sum(cumper > upperbnd) + sum(cumper < lowerbnd);

% Convert to a percentage inside confidence interval.

pct_in = (1-(test/mhalf))*100;

yes = (pct_in >= 95);

% Graph the cumulative periodogram.

plot([0:mhalf-1]/(2*mhalf),cumper,'r')
hold on
plot([0:mhalf-1]/(2*mhalf),upperbnd,'-')
plot([0:mhalf-1]/(2*mhalf),lowerbnd,'-')
axis([0 .5 0 1])
text(.03,.93,'95% confidence interval')
hold off
drawnow

