%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%  What are the 95% confidence limits on the true variance of minimum
%  temperature during January, 1980 (the data from HW1)?

clear

%  Load Madison daily temperature
dat = load('~/matlab/Data/AOS575/HW1_prob2_data_noheader.txt');

%  Take the 4th column 
temp = dat(:,4);

%  Get the sample variance
varsamp = var(temp);

%  Get the sample size
N = length(temp);

%  Compute the confidence limits
varmin = (N-1)*varsamp/chi2inv(0.975, N-1);
varmax = (N-1)*varsamp/chi2inv(0.025, N-1);

%  So, we state that there is a 95% chance that the true variance is
%  between varmin and varmax
varlims = [varmin varmax]



%  Another example:  confidence limits for a power spectrum

clear

%  Define AR1 time series with autocorrelation rho = sqrt(0.5);
  ntim = 1000;
  rho = sqrt(0.5);
  x = randn(ntim,1); x = standardize(x);  %  Noise
  
  %  Generate time series
  y = repmat(NaN, [ntim, 1]);  %  Actual AR1 time series
  y(1) = x(1);
  for i = 2:ntim;
    y(i) = rho*y(i-1)+sqrt(1-rho^2)*x(i);  %  red noise w/ unit var
  end
  y = standardize(y);
  
  %  Now, add a true peak with variance 0.2:
  y = y + 0.4*sin(2*pi*[1:ntim]'/10);
  
  %  Generate power spectrum
  nfft = 100;
  [p,f] = spectrum(y, nfft, nfft/2, hamming(nfft), 1);
  
%  Now, construct theoretical AR1 spectrum
  pAR1 = 0.5./(1+rho^2-2*rho*cos(2*pi*[0:1/nfft:0.5]));
  

%  OK!  Let's go through hypothesis test

%  Significance level:
siglev = 0.90;

%  H0:  The spectral power at 0.1 is no different than red noise.
%  H1:  The spectrum has more power than red noise at f=0.1.

%  Statistic:  One-tailed Chi2 test with 'dof' degrees of freedom
dof = 2*(ntim/nfft);

%  Define critical region:  p < psig95 ==> accept H0
psig95 = pAR1*chi2inv(siglev, dof)/(dof-1);

%  Evaluate statistic:  or, plot curves
figure_tall(1); clf;  %  Set up "tall" figure window
global_axes(3, 1.7, 0.5, 0.5, 1.5);  %  Set up axis geometry

subplot2(1,1);
  pp = plot(log(f), log(p(:,1)), 'k', ...
	    log(f), log(pAR1), 'k', ...
	    log(f), log(psig95), '--k');
  set(pp, 'linewidth', 1);
  set(pp(1), 'linewidth', 2);
  axis([log(1/nfft) log(0.5) -2.5 2.5]);
  set(gca, 'XTick', log(2.^[-8:-1]), 'XTickLabel', 2.^[8:-1:1]);
  xl = xlabel('Period', 'fontsize', 10);
  yl = ylabel('log(P_y_y)', 'fontsize', 10);
  t1 = title2('{\bf\chi^2 Example}', 'fontsize', 12);
  grid on;
  
  %  For kicks, draw a vertical line at f=0.1
  hold on;
    ll = line(log(0.1)*[1 1], [-2.5 2.5]);
    set(ll, 'color', 'k', 'linewidth', 2, 'linestyle', '--');
  hold off;
  
%  Print it up