Outlier rejection plot strategy

5 views (last 30 days)
HIRAKJYOTI BASUMATARY
HIRAKJYOTI BASUMATARY on 15 Oct 2017
I am new to matlab and i am struggling with a simple figure plot strategy. I have written a code for regression diagnostics and outlier rejection for Linear regression model using chebyshev inequality (i.e if the error is greater than mean+std deviation , then it can be regarded as outlier). Through the code, i want to show the outlier detected in each iteration by blue dot(or circle) and also show the fitting line in each iteration. How to go for this technique of matlab plot in a efficient way?
function trialregressiondiagnostics
xmin=-10;
xmax=10;
n=7;
x1=xmin+rand(n,1)*(xmax-xmin);
a=2; b=3;
for i =1:n
y(i,1)=a*x1(i,1)+b;
end
% Simulation of random noise around this line
noise=randn(7,1); % Generating numbers from a N(0,1)dist
S = std(noise); % Standard deviation of the noise
ynoisy=y+2*S*noise;% Making the noisy data
hold on;
plot(x1,ynoisy,'bx')
%%outlier
outlier=rand(3,1);
S1= std(outlier);
Youtlier = ynoisy(1:3,1)+10*S1*outlier;
hold on;
Ytotal = horzcat(ynoisy',Youtlier');
size(Ytotal')
%%line fit
x2=xmin+rand(3,1)*(xmax-xmin);
x=horzcat(x1',x2');
m=length(x);
X = [x' ones(m,1)];
Y = [Ytotal'];
P = (((X'*X)^-1)*X')*Y;
hold on;
plot(x,Ytotal,'rx');
%%line fit
b0=P(1); b1=P(2);
xfit= linspace(-10,10, 10);
yfit= b0*xfit+b1;
hold on;
plot(xfit,yfit,'k','LineWidth',2)
%%error
k=0;
while k<1000;
if length(Ytotal)<5
break
end
yerror=abs(Ytotal-yfit);
[yerror idx]=sort(yerror,'descend');
Ytotal=Ytotal(idx);
x=x(idx);
ymean = mean(yerror);
yerrorstd=std(yerror);
num_total = length(Ytotal);
num_ready = 0;
num_position = 1;
zz=Ytotal;
for num_position=1:length(zz)
if num_position>length(Ytotal)
break
end
if yerror(num_position) > ymean+yerrorstd
Ytotal(num_position)=[];
x(num_position)=[];
end
end
mm=length(x);
X1 = [x' ones(mm,1)];
Y1 = [Ytotal]';
PP = (((X1'*X1)^-1)*X1')*Y1;
B0=PP(1); B1=PP(2);
xfit= linspace(-10,10,length(Ytotal));
yfit= B0*xfit+B1;
k=k+1;
figure
plot(x,Ytotal,'rx');
hold on;
plot(xfit,yfit,'r','LineWidth',2);
pause(1);
end
end

Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!