close all
clear all
%reset random number generator - will generate same random sequences every
%time
rng default

%make x vector
x=[-1:.002:1];
n=length(x);

%set up and make "signal"
m=4;
b=-1;
noise_amp=1;
y=m*x+b;

%make noise
noise=noise_amp*randn(1,n);
%fit noise - hopefully will have zero slope and 0 intercept! this is a test
%it is not used anywhere else in script.
[pnoise,Snoise] = polyfit(x,noise,1);
%The coefficients in p are in descending powers, and the length of p is n+1
%structure S to be used as an input to polyval to obtain model

%add noise to signal
noisy=y+noise;

%fit signal plus noise, get model from fit
[pnoisy,Snoisy] = polyfit(x,noisy,1);
[modlnoisy deltanoisy]=polyval(pnoisy,x,Snoisy);

%set limit on standard deviation for selecting/rejecting data
SigScale=2.;
%calc standard deviation
stdnoisy=std(noisy-modlnoisy);

figure(1)
FullPlot(x,y,noisy,modlnoisy,SigScale*stdnoisy,pnoisy)

%add some really bad points - generate a random set of x points and then
%random data for those points
%percent of bad points
PctBad=.01;
NBad=fix(PctBad*n);
%generate random indices for bad points
IndxBad=fix(n*rand(NBad,1));
%generate bad data for bad points (have to scale it up from regular range
%of randn output
BadScale=10;
BadPts=BadScale*randn(1,NBad); %all data has to be row vector
%now make new vector with the bad points added
BadNoisy=noisy; %have to fill it first - so can change indexed points
%indices that have bad points - note where we need and don't need IndxBad
BadNoisy(IndxBad)=noisy(IndxBad)+BadPts; %watch the indexing here

%now fit and get model
[pBadNoisy,SBadNoisy] = polyfit(x,BadNoisy,1);
[modlBadNoisy deltaBadNoisy]=polyval(pBadNoisy,x,SBadNoisy);
stdBadNoisy=std(BadNoisy-modlBadNoisy);

figure(2)
FullPlot(x,y,BadNoisy,modlBadNoisy,SigScale*stdBadNoisy,pBadNoisy)

%find points outside 2 sd and remove them - do this by SELECTING AND
%KEEPING points within 2 sd.

%here you need to put the test into the find function
[row col val]=find(XXXX);
%here you need to calculate the linear indicies lf the selected elements
XXXX = sub2ind(XXXX);
%here you need to make a new set of input vectors, with only the selected
%elements, for the fit (polyfit), the calculation of the model (polyval)
%from the fit, and the part of the original y vector with just the
%selected elelments for plotting
cleanedBadNoisy=BadNoisy(XXXX);
%need to make new x vector (and y vector for plot - but only use once 
%below) with just the cleaned points
xclean=x(XXXX);
yclean=y(XXXX);

%now fit and get model
[pcleanedBadNoisy,ScleanedBadNoisy] = polyfit(xclean,cleanedBadNoisy,1);
[modlcleanedBadNoisy deltacleanedBadNoisy]=polyval(pcleanedBadNoisy,xclean,ScleanedBadNoisy);
stdmodlcleanedBadNoisy=std(cleanedBadNoisy-modlcleanedBadNoisy);

figure(3)
FullPlot(xclean,yclean,cleanedBadNoisy,modlcleanedBadNoisy,SigScale*stdmodlcleanedBadNoisy,pcleanedBadNoisy)

%

function FullPlot(x,truth,data,model,delta,p)
plot(x,truth,'b.')
hold on
grid
xlabel('x')
ylabel('data - blue, noisy data - red, model green')
plot(x,data,'r.')
plot(x,model,'g.')
plot(x,model+delta,'m--',x,model-delta,'m--')
titletext=strcat("Lin Fit with 95% Prediction Interval, p= ",num2str(p)," delta ",num2str(delta));
title(titletext)
legend('Data','Linear Fit','95% Prediction Interval')
end