% script testSe.m
% sets up dummy regression problems with various kinds of
% additive noise: normal, uniform(-1,1)*lognormal, signed lognormal with
% random sign, Tukey's "slash" (normal/uniform(0,1)), and a mixture of
% a normal and a lognormal with mixture fraction _mix_.
%
% To control which error distribution is used, set the parameter
% _noiseType_ to one of 'normal' 'logByU' 'logMixNorm'
% 'signLog' 'slash'.
%
% Defaults to "normal"
%
% Estimates the bias and the bias of the bootstrap SE estimate
%
%
% P.B. Stark stark@stat.berkeley.edu
% 2 August 1997

% parameters
ndat = 100;			% number of data
snr = 5;			% signal to noise ratio
bdim = 5;			% dimension of model
sc = 'one-step';		% procedure for estimating scale of
				% the residuals			
ftype = 'bldCheb';		% subroutine to build design matrix
wght1 = 'biwght';		% subroutine for robust weights
param1 = 8;			% parameter for the weight
wght2 = 'hampel';
param2 = [2 4 8];
wght3 = 'l1wght';
param3 = 1;
nboot = 100;			% number of bootstrap replications
nrep = 100;			% number of replications in 1st generation
mix = 0.02;			% the fraction of lognormal to mix with the
				% normal in the mixture error model

x = [1:ndat];
X = feval(ftype,ndat,bdim-1);	% build design matrix
XP = X(ceil(ndat/2),:);		% prediction matrix
pred = zeros(nrep,2);		% matrix of predictions
seHat = zeros(nrep,2);		% matrix of resampling-estimated se's

% generate random b
b = 10*randn(bdim,1);		% normal (for now)

yt = X*b;
normy = norm(yt);		% nominal signal level
defaultNoise = 'normal';

if(~exist('noiseType')),
	disp(['This script requires a value for noiseType.'])
	disp(['Setting noiseType to default value: ' defaultnoise])
	noiseType = defaultNoise;
end

if(strcmp(noiseType,'normal')),
% Normal
	noise = normy/snr*randn(ndat,nrep); 		% normal
%
elseif (strcmp(noiseType,'logByU')),
%
% Lognormal times uniform
	noise = 2*normy/snr*(rand(ndat,nrep) -0.5).* ...% U*lognormal
	exp(randn(ndat,nrep));
%
elseif (strcmp(noiseType,'logMixNorm')),
%
% Signed lognormal mixed with normal, mixture fraction of
% _mix_ of the lognormal
	umat = rand(ndat,nrep);				% uniform mix prob
	noise = normy/snr*randn(ndat,nrep);
	lgn = normy/snr*exp(randn(ndat,nrep))*2 ...	% lognormal w/ rand
		.*((rand(ndat,nrep) >= 0)-0.5);		% sign
	noise(umat < mix) = lgn(umat < mix);		% make the mix
%
elseif(strcmp(noiseType,'signLog')),
% Lognormal with random sign
	noise = normy/snr*exp(randn(ndat,nrep))*2 ...	% lognormal w/ rand
		.*((rand(ndat,nrep) >= 0)-0.5);		% sign
%
elseif(strcmp(noiseType,'slash')),
%
% Tukey's "slash:" normal divided by uniform
	noise = normy/snr*randn(ndat,nrep)./rand(ndat,nrep);  % "slash"
%
else
	error(['Noise type ' noiseType ' not implemented.'])
end

% 
% loop over replications

for i = 1:nrep,
	disp(['replication ' num2str(i)])
	y = yt + noise(:,i);
	[b1, res1, seHat(i,1), it1 ] = ...
		rwlsBoot(y,X,XP,nboot,wght1,param1,sc);
	pred(i,1) = XP*b1;
	[b2, res2, seHat(i,2), it2 ] = ...
		rwlsBoot(y,X,XP,nboot,wght2,param2,sc);
	pred(i,2) = XP*b2;
end;

ePred = mean(pred);
truth = XP*b;
bias = ePred - truth
relBias = bias/truth
seTrue = std(pred)
meanSeHat = mean(seHat)
seBias = meanSeHat-seTrue
seRelBias = seBias ./ seTrue