% script magsTest
% loads the Magsat data and fits parts of it
% using robust regression, with bootstrap uncertainty estimates.
% Plots estimates +- 2SE, color-coded for coverage of the datum
%
%
% P.B. Stark  stark@stat.berkeley.edu
% 15 October 1997.

% parameters
ndat = 25;			% number of data in each window
mid = ceil(ndat/2);		% half-width of window
frst = 20;			% window to start with
bdim = 4;			% dimension of model (cubic)
ftype = 'bldCheb';		% subroutine to build design matrix
scalMeth = 'one-step';		% way to estimate the scale
wght1 = 'biwght';		% biweight is weight function 1
par1 = 8;			% parameter of weight
wght2 = 'l1wght';		% second weight function
par2 = 1;
wght3 = 'hampel';		% number 3
par3 = [2 4 8];
nints = 1;			% number of time intervals of length
				% ndat to compute estimates for
nboot = 100;			% number of bootstrap samples
baddat = 99999;			% distinctive value for missing MAGSAT
				% data
iplsym = 'r+';			% plotting symbols and colors for the
jplsym = 'b+';			% error bars
kplsym = 'c+';			%

% load calibration data

if (~exist('redo')),
   redo = 0;
end

if (redo == 0),
	redo = 1;
	load magsatQ.dat -ascii	% quiet MAGSAT data
	load magsatD.dat -ascii % disturbed MAGSAT data
% resize the data vectors
	[m n] = size(magsatq);
	magsatQ = reshape(magsatq',m*n,1);
	magsatD = reshape(magsatd',m*n,1);
	QN = magsatQ(1:1024);	% north component measurement
	QE = magsatQ(1025:2048);% east component
	QV = magsatQ(2049:3072);% vertical component
	Q = [QN QE QV];		% all quiet data
	DN = magsatD(1:1024);	% north component measurement
	DE = magsatD(1025:2048);% east component
	DV = magsatD(2049:3072);% vertical component
	D = [DN DE DV];		% all disturbed data
	t = 1:1024;
end


x = [1:ndat];
X = feval(ftype,ndat,bdim-1);	% build design matrix
XP = X(mid,:);			% prediction vector for middle sample

ncoveri = 0;
ncoverj = 0;
ncoverk = 0;
titlab = ['N' 'E' 'V'];
for i=1:nints,
	for j = 1:3,
		start = (frst+i-1)*ndat + 1;
		y = Q(start : start+ndat-1, j);
% censor the irrelevant rows of the design matrix caused by missing data
% the flag for missing data is baddat
		XU = X((y ~= baddat),:);
		yu = y(y ~= baddat);
		xu = x(y ~= baddat);
		figure
		plot(xu+start-1,yu,'go')
		title(['Quiet MAGSAT data: 17 Feb 1980, ' titlab(j) ...
			' component'])
		xlabel('Sample after 23,751s on MJD 44286')
		hold
		[bbi, resbi, sebi ] = rwlsBoot(yu,XU,XP,nboot,wght1,par1, ...
			scalMeth);		% irwls w/ weight 1
		[bbj, resbj, sebj ] = rwlsBoot(yu,XU,XP,nboot,wght2,par2, ...
			scalMeth);		% irwls weight 2
		[bbk, resbk, sebk ] = rwlsBoot(yu,XU,XP,nboot,wght3,par3, ...
			scalMeth);		% irwls weight 3		
		[bols, resols, wssols] = wls(yu, XU, ones(size(yu))); 
			% ordinary least-squares
%
		plot(xu+start-1, XU*bols,'w')
		plot(xu+start-1, XU*bbi,'r--')
		plot(xu+start-1, XU*bbj,'b-.')
		plot(xu+start-1, XU*bbk,'c:')

		ncoveri = ncoveri + (abs(y(mid) - XP*bbi) <= 2*sebi);
		ncoverj = ncoverj + (abs(y(mid) - XP*bbj) <= 2*sebj);
		ncoverk = ncoverk + (abs(y(mid) - XP*bbk) <= 2*sebk);
		plot(x(mid)+start-1, XP*bbi + 2*sebi,iplsym)
		plot(x(mid)+start-1, XP*bbi - 2*sebi,iplsym)
		plot(x(mid)+start-1, XP*bbj + 2*sebj,jplsym)
		plot(x(mid)+start-1, XP*bbj - 2*sebj,jplsym)
		plot(x(mid)+start-1, XP*bbk + 2*sebk,kplsym)
		plot(x(mid)+start-1, XP*bbk - 2*sebk,kplsym)	
% plot residuals
		figure
		hold	
		title(['Quiet MAGSAT residuals: 17 Feb 1980, ' titlab(j) ...
			' component'])
		xlabel('Sample after 23,751s on MJD 44286')
		plot(xu+start-1, resols,'w')
		plot(xu+start-1, resbi,'r--')
		plot(xu+start-1, resbj,'b-.')
		plot(xu+start-1, resbk,'c:')

		plot(x(mid)+start-1, 2*sebi,iplsym)
		plot(x(mid)+start-1, -2*sebi,iplsym)
		plot(x(mid)+start-1, 2*sebj,jplsym)
		plot(x(mid)+start-1, -2*sebj,jplsym)
		plot(x(mid)+start-1, 2*sebk,kplsym)
		plot(x(mid)+start-1, -2*sebk,kplsym)	

	end;
	for j = 1:3,
		start = (frst+i-1)*ndat + 1;
		y = D(start : start+ndat-1, j);
% censor the irrelevant rows of the design matrix caused by missing data
% the flag for missing data is baddat
		XU = X((y ~= baddat),:);
		yu = y(y ~= baddat);
		xu = x(y ~= baddat);
		figure
		plot(xu+start-1,yu,'go')
		title(['Disturbed MAGSAT Data: 16 Feb 1980, ' titlab(j) ...
			' component'])
		xlabel('sample after 9,939s on MJD 44285')
		hold
		[bbi, resbi, sebi ] = rwlsBoot(yu,XU,XP,nboot,wght1,par1, ...
			scalMeth);		% irwls w/ weight 1
		[bbj, resbj, sebj ] = rwlsBoot(yu,XU,XP,nboot,wght2,par2, ...
			scalMeth);		% irwls weight 2
		[bbk, resbk, sebk ] = rwlsBoot(yu,XU,XP,nboot,wght3,par3, ...
			scalMeth);		% irwls weight 3	
		[bols, resols, wssols] = wls(yu, XU, ones(size(yu))); 
			% ordinary least-squares
%
		plot(xu+start-1, XU*bols,'w')
		plot(xu+start-1, XU*bbi,'r--')
		plot(xu+start-1, XU*bbj,'b-.')
		plot(xu+start-1, XU*bbk,'c:')

		ncoveri = ncoveri + (abs(y(mid) - XP*bbi) <= 2*sebi);
		ncoverj = ncoverj + (abs(y(mid) - XP*bbj) <= 2*sebj);
		ncoverk = ncoverk + (abs(y(mid) - XP*bbk) <= 2*sebk);
		plot(x(mid)+start-1, XP*bbi + 2*sebi,iplsym)
		plot(x(mid)+start-1, XP*bbi - 2*sebi,iplsym)
		plot(x(mid)+start-1, XP*bbj + 2*sebj,jplsym)
		plot(x(mid)+start-1, XP*bbj - 2*sebj,jplsym)
		plot(x(mid)+start-1, XP*bbk + 2*sebk,kplsym)
		plot(x(mid)+start-1, XP*bbk - 2*sebk,kplsym)

		figure
		hold	
		title(['Disturbed MAGSAT residuals: 17 Feb 1980, ' titlab(j) ...
			' component'])
		xlabel('Sample after 23,751s on MJD 44286')
		plot(xu+start-1, resols,'w')
		plot(xu+start-1, resbi,'r--')
		plot(xu+start-1, resbj,'b-.')
		plot(xu+start-1, resbk,'c:')

		plot(x(mid)+start-1, 2*sebi,iplsym)
		plot(x(mid)+start-1, -2*sebi,iplsym)
		plot(x(mid)+start-1, 2*sebj,jplsym)
		plot(x(mid)+start-1, -2*sebj,jplsym)
		plot(x(mid)+start-1, 2*sebk,kplsym)
		plot(x(mid)+start-1, -2*sebk,kplsym)		
	end;
end;

coveri = ncoveri/(6*nints)
coverj = ncoverj/(6*nints)
coverk = ncoverk/(6*nints)
redo = 1;