% script wrongMod.m
% estimates the relative bias that results from mis-specifying the 
% polynomial degree of the model, using 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 signed lognormal.
%
% To control which error distribution is used, set the parameter
% _noiseType_ to one of 'normal' 'logByU' 'logMixNorm'
% 'signLog' 'slash'.
%
% Defaults to "normal"
%
%
%
% P.B. Stark stark@stat.berkeley.edu
% 3 August 1997.

% parameters
ndat = 100;			% number of data
snr = 5;			% signal to noise ratio
bdim = 5;			% dimension of fitted model
loBdim = 3;			% small true model dimension
hiBdim = 7;			% large true model dimension
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 = 0;			% number of bootstrap replications
nrep = 100;			% number of replications in 1st generation
nmod = 20;			% number of models
mix = 0.02;			% fraction of lognormal to mix in

x = [1:ndat];
X = feval(ftype,ndat,bdim-1);	% build design matrix
XP = X(ceil(ndat/2),:);		% prediction matrix
Xlo = feval(ftype,ndat,loBdim-1); % design matrix for small dimension
XPlo = Xlo(ceil(ndat/2),:);
Xhi = feval(ftype,ndat,hiBdim-1); % design matrix for large dimension
XPhi = Xhi(ceil(ndat/2),:);

defaultNoise = 'normal';	  % default type of additive noise

aveBiasLo = [0, 0];		  % vector of biases for low dim truth
aveAbsBiasLo = [0, 0];		  % vector of abs bias, low dim truth
aveBias = [0,0];		  % vector of bias, true dim
aveAbsBias = [0, 0];		  % vector of abs bias, true dim
aveBiasHi = [0, 0];		  % vector of bias, high dim truth
aveAbsBiasHi = [0, 0];		  % vector of abs bias, high dim truth


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 = 1.0/snr*randn(ndat,nrep); 		% normal
%
elseif (strcmp(noiseType,'logByU')),
%
% Lognormal times uniform
	noise = 2.0/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 = 1.0/snr*randn(ndat,nrep);
	lgn = 1.0/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 = 1.0/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 = 1.0/snr*randn(ndat,nrep)./rand(ndat,nrep);  % "slash"
%
else
	error(['Noise type ' noiseType ' not implemented.'])
end

%
% loop over realizations of the model
for i = 1:nmod,
	disp(['model replication ' num2str(i)])
% generate random b
	b = 10*randn(hiBdim,1);		% normal (for now)

	ythi = Xhi*b;			% noise-free data for high dimension
	ytlo = Xlo*b(1:loBdim);		% noise-free data for low dimension
	yt = X* b(1:bdim);		% noise-free data, right dimension
	normyhi = norm(ythi);		% nominal signal level, high dim.
	normylo = norm(ytlo);		% nominal signal level, low dim.predLo = zeros(nrep,2);
	normy = norm(yt);
	predHi = zeros(nrep,2);
	pred = zeros(nrep, 2);
	predLo = zeros(nrep,2);

% loop over noise realizations for a fixed model
	for j = 1:nrep,
% lower dimension truth than fitted
		y = ytlo + normylo*noise(:,j);
		[b1, res1, se1, it1 ] = ...
			rwlsBoot(y,X,XP,nboot,wght1,param1,sc);
		predLo(j,1) = XP*b1;
		[b2, res2, se2, it2 ] = ...
			rwlsBoot(y,X,XP,nboot,wght2,param2,sc);
		predLo(j,2) = XP*b2;
% correct dimension truth
		y = yt + normy*noise(:,j);
		[b1, res1, se1, it1 ] = ...
			rwlsBoot(y,X,XP,nboot,wght1,param1,sc);
		pred(j,1) = XP*b1;
		[b2, res2, se2, it2 ] = ...
			rwlsBoot(y,X,XP,nboot,wght2,param2,sc);
		pred(j,2) = XP*b2;		
% larger dimension truth than fitted
		y = ythi + normyhi*noise(:,j);
		[b1, res1, se1, it1 ] = ...
			rwlsBoot(y,X,XP,nboot,wght1,param1,sc);
		predHi(j,1) = XP*b1;
		[b2, res2, se2, it2 ] = ...
			rwlsBoot(y,X,XP,nboot,wght2,param2,sc);
		predHi(j,2) = XP*b2;
	
	end;
	relBiasLo = (mean(predLo) - ...
		XPlo*b(1:loBdim))/(XPlo*b(1:loBdim));
	relBiasHi = (mean(predHi) - XPhi*b)/(XPhi*b);
	relBias = (mean(pred) - XP*b(1:bdim))/(XP*b(1:bdim));
	aveBiasLo = aveBiasLo + relBiasLo;
	aveAbsBiasLo = aveAbsBiasLo + abs(relBiasLo);
	aveBiasHi = aveBiasHi + relBiasHi;
	aveAbsBiasHi = aveAbsBiasHi + abs(relBiasHi);
	aveBias = aveBias + relBias;
	aveAbsBias = aveAbsBias + abs(relBias);
end;

aveBiasLo = aveBiasLo/nmod
aveAbsBiasLo = aveAbsBiasLo/nmod
aveBias = aveBias/nmod
aveAbsBias = aveAbsBias/nmod
aveBiasHi = aveBiasHi/nmod
aveAbsBiasHi = aveAbsBiasHi/nmod