-
Notifications
You must be signed in to change notification settings - Fork 4
/
xyrobustregress.m
executable file
·37 lines (34 loc) · 1.38 KB
/
xyrobustregress.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
function [b,bstat] = xyrobustregress(x,y,alpha)
%ROBUSTREGRESS combination of regress(y(:),[ones(size(x,1),1) x(:)]) robustfit(x(:),y(:)) to generate confidence interval not sensitive to outliers
% Syntax: b = robustregress(x,y [,alpha])
% [b,bstat] = robustregress(...)
% x,y: nx1 vector
% alpha: tp compute 100*(1-alpha)% confidence level (default = 0.05)
% b: [intercept slope]
% bstat: [ interceptmin interceptmax
% slopemin slopemax]
% MS 2.1 - 02/03/12 - Olivier Vitrac - rev. 07/03/2016
% revision history
% 07/03/2016 fix n=2 (robustfit was not working), manage NaN values
if nargin<2, error('2 arguments are required'), end
if nargin<3, alpha = 0.05; end
x = x(:); y = y(:);
ok = ~isnan(y); x = x(ok); y = y(ok);
if numel(x) ~=numel(y), error('x and y must have the same size'); end
n = length(x);
[b1,b1stat] = regress(y,[ones(n,1) x]);
if n>2
[b2,b2stat] = robustfit(x,y);
tval = tinv((1-alpha(1)/2),b2stat.dfe);
b2stat = [b2-tval*b2stat.se, b2+tval*b2stat.se];
performance = mean([(b1stat(:,2)-b2)./abs(b1) (b2stat(:,2)-b2)./abs(b2)]);
else
performance = [0 Inf];
end
if performance(2)<performance(1)/1.05 % confidence interval reduced by at least 5% when robustfit is used instead regress
b = b2;
if nargout>1, bstat = b2stat; end
else
b = b1;
if nargout>1, bstat = b1stat; end
end