Сравните устойчивые методы регрессии

Этот пример сравнивает результаты среди методов регрессии, которые являются и не устойчивы к влиятельным выбросам.

Влиятельные выбросы являются экстремальным ответом или наблюдениями предиктора что оценки параметра влияния и выводы регрессионного анализа. Ответы, которые являются влиятельными выбросами обычно, происходят в экстремальных значениях области. Например, у вас может быть инструмент, который измеряет ответ плохо или беспорядочно на экстремальных уровнях температуры.

С достаточным количеством доказательства можно удалить влиятельные выбросы из данных. Если удаление не возможно, можно использовать методы регрессии, которые устойчивы к выбросам.

Симулируйте данные

Сгенерируйте случайную выборку 100 ответов из линейной модели

$y_{t} = 1 + 2 x_{t} + ε_{t}$

где:

$x$ вектор равномерно расположенных с интервалами значений от 0 до 2.
$ε_{t} \sim N (0, 0 . 5^{2}) .$

rng('default');
n = 100;
x = linspace(0,2,n)';
b0 = 1;
b1 = 2;
sigma = 0.5;
e = randn(n,1);
y = b0 + b1*x + sigma*e;

Создайте влиятельные выбросы путем расширения всего соответствия ответов $x < 0 . 25$ фактором 3.

y(x < 0.25) = y(x < 0.25)*3;

Отобразите данные на графике. Сохраните график для дальнейших графиков.

figure;
plot(x,y,'o');
h = gca;
xlim = h.XLim';
hl = legend('Data');
xlabel('x');
ylabel('y');
title('Regression Techniques Comparison')
hold on;

Оцените линейную модель

Оцените коэффициенты и ошибочное отклонение при помощи простой линейной регрессии. Постройте линию регрессии.

LSMdl = fitlm(x,y)

LSMdl = 
Linear regression model:
    y ~ 1 + x1

Estimated Coefficients:
                   Estimate      SE       tStat       pValue  
                   ________    _______    ______    __________

    (Intercept)     2.6814     0.28433    9.4304    2.0859e-15
    x1             0.78974     0.24562    3.2153     0.0017653


Number of observations: 100, Error degrees of freedom: 98
Root Mean Squared Error: 1.43
R-squared: 0.0954,  Adjusted R-Squared: 0.0862
F-statistic vs. constant model: 10.3, p-value = 0.00177

plot(xlim,[ones(2,1) xlim]*LSMdl.Coefficients.Estimate,'LineWidth',2);
hl.String{2} = 'Least Squares';

LSMdl подходящий LinearModel объект модели. Прерывание и наклон, кажется, соответственно выше и ниже, чем они должны быть. Линия регрессии может сделать плохие прогнозы для любого $x < 1$ и $x > 1.6$ .

Оцените байесовую модель линейной регрессии с рассеянным предшествующим распределением

Создайте Байесовую модель линейной регрессии с рассеянным соединением, предшествующим для ошибочного отклонения и коэффициентов регрессии. Задайте один предиктор для модели.

PriorDiffuseMdl = bayeslm(1);

PriorDiffuseMdl diffuseblm объект модели, который характеризует объединенное предшествующее распределение.

Оцените следующую из Байесовой модели линейной регрессии. Постройте линию регрессии.

PosteriorDiffuseMdl = estimate(PriorDiffuseMdl,x,y);

Method: Analytic posterior distributions
Number of observations: 100
Number of predictors:   2
 
           |  Mean     Std         CI95        Positive      Distribution     
------------------------------------------------------------------------------
 Intercept | 2.6814  0.2873  [ 2.117,  3.246]    1.000   t (2.68, 0.28^2, 98) 
 Beta      | 0.7897  0.2482  [ 0.302,  1.277]    0.999   t (0.79, 0.25^2, 98) 
 Sigma2    | 2.0943  0.3055  [ 1.580,  2.773]    1.000   IG(49.00, 0.0099)

plot(xlim,[ones(2,1) xlim]*PosteriorDiffuseMdl.Mu,'--','LineWidth',2);
hl.String{3} = 'Bayesian Diffuse';

PosteriorDiffuseMdl conjugateblm объект модели, который характеризует объединенное апостериорное распределение линейных параметров модели. Оценки Байесовой модели линейной регрессии с предшествующим рассеянным почти равны тем из простой модели линейной регрессии. Обе модели представляют наивный подход к влиятельным выбросам, то есть, методы обрабатывают выбросы как любое другое наблюдение.

Оцените модель регрессии с ошибками ARIMA

Создайте модель регрессии с ошибками ARIMA. Укажите, что ошибки следуют за t распределением с 3 степенями свободы, но никакими изолированными условиями. Эта спецификация является эффективно моделью регрессии с распределенными ошибками t.

regMdl = regARIMA(0,0,0);
regMdl.Distribution = struct('Name','t','DoF',3);

regMdl regARIMA объект модели. Это - шаблон для оценки.

Оцените модель регрессии с ошибками ARIMA. Постройте линию регрессии.

estRegMdl = estimate(regMdl,y,'X',x);

 
    Regression with ARMA(0,0) Error Model (t Distribution):
 
                  Value     StandardError    TStatistic      PValue  
                 _______    _____________    __________    __________

    Intercept     1.4613         0.154         9.4892       2.328e-21
    Beta(1)       1.6531       0.12939         12.776      2.2245e-37
    Variance     0.93116        0.1716         5.4263      5.7546e-08
    DoF                3             0            Inf               0

plot(xlim,[ones(2,1) xlim]*[estRegMdl.Intercept; estRegMdl.Beta],...
    'LineWidth',2);
hl.String{4} = 'regARIMA';

estRegMdl regARIMA объект модели, содержащий результаты оценки. Поскольку t распределение - больше рассеянное, линия регрессии приписывает больше изменчивости влиятельным выбросам, чем к другим наблюдениям. Поэтому линия регрессии, кажется, лучшая прогнозная модель, чем другие модели.

Реализуйте регрессию квантиля Используя мешок деревьев регрессии

Вырастите мешок 100 деревьев регрессии. Задайте 20 для минимального листового размера.

QRMdl = TreeBagger(100,x,y,'Method','regression','MinLeafSize',20);

QRMdl подходящий TreeBagger объект модели.

Предскажите средние ответы для всех наблюдаемых $x$ значения, то есть, регрессия квантиля реализации. Постройте прогнозы.

qrPred = quantilePredict(QRMdl,x);
plot(x,qrPred,'LineWidth',2);
hl.String{5} = 'Quantile';
hold off;

Линия регрессии, кажется, немного под влиянием выбросов в начале выборки, но затем быстро следует за regARIMA модельный ряд.

Можно настроить поведение линии путем определения различных значений для MinLeafSize когда вы обучаете мешок деревьев регрессии. Более низкий MinLeafSize значения имеют тенденцию следовать за данными в графике более тесно.

Реализуйте устойчивую байесовую линейную регрессию

Рассмотрите Байесовую модель линейной регрессии, содержащую один предиктор, t распределил отклонение воздействия профилируемым параметром степеней свободы $ν$ Пусть:

$λ_{j} \sim I G (ν / 2, 2 / ν)$ .
$ε_{j} | λ_{j} \sim N (0, λ_{j} σ^{2})$
$f (β, σ^{2}) \propto \frac{1}{σ^{2}}$

Эти предположения подразумевают:

$ε_{j} \sim t (0, σ^{2}, ν)$
$λ_{j} | ε_{j} \sim I G (\frac{ν + 1}{2}, \frac{2}{ν + ε_{j}^{2} / σ^{2}})$

$λ$ вектор скрытых масштабных коэффициентов, который приписывает низкую точность наблюдениям, которые далеки от линии регрессии. $ν$ гиперпараметр, управляющий влиянием $λ$ на наблюдениях.

Задайте сетку значений для $ν$ .

1 соответствует распределению Коши.
2.1 средние значения, что среднее значение четко определено.
4,1 средних значения, что отклонение четко определено.
100 средние значения, что распределение приблизительно нормально.

nu = [0.01 0.1 1 2.1 5 10 100];
numNu = numel(nu);

Для этой проблемы сэмплер Гиббса является подходящим, чтобы оценить коэффициенты, потому что можно симулировать параметры Байесовой модели линейной регрессии, обусловленной на $λ$ , и затем симулируйте $λ$ от его условного распределения.

Реализуйте этот сэмплер Гиббса.

Чертите параметры от апостериорного распределения $β, σ^{2} | y, x, λ$ . Выкачайте наблюдения $λ$ , создайте рассеянную предшествующую модель с двумя коэффициентами регрессии и чертите набор параметров от следующего. Первый коэффициент регрессии соответствует прерыванию, поэтому задайте тот bayeslm не включают тот.
Вычислите остаточные значения.
Чертите значения от условного выражения, следующего из $λ$ .

Для каждого значения $ν$ , запустите сэмплер Гиббса для 20 000 итераций и примените электротермотренировку 5 000. Предварительно выделите для следующих ничьих и инициализируйте $λ$ к вектору из единиц.

rng(1);
m = 20000;
burnin = 5000;
lambda = ones(n,m + 1,numNu); % Preallocation
estBeta = zeros(2,m + 1,numNu);
estSigma2 = zeros(1,m + 1,numNu);

% Create diffuse prior model.
PriorMdl = bayeslm(2,'Model','diffuse','Intercept',false);

for p = 1:numNu
    for j = 1:m
        % Scale observations.
        yDef = y./sqrt(lambda(:,j,p)); 
        xDef = [ones(n,1) x]./sqrt(lambda(:,j,p));
        
        % Simulate observations from conditional posterior of beta and
        % sigma2 given lambda and the data.
        [estBeta(:,j + 1,p),estSigma2(1,j + 1,p)] = simulate(PriorMdl,xDef,yDef);
        
        % Estimate residuals.
        ep = y - [ones(n,1) x]*estBeta(:,j + 1,p);
        
        % Specify shape and scale using conditional posterior of lambda
        % given beta, sigma2, and the data
        sp = (nu(p) + 1)/2;
        sc =  2./(nu(p) + ep.^2/estSigma2(1,j + 1,p));
        
        % Draw from conditional posterior of lambda given beta, sigma2, 
        % and the data
        lambda(:,j + 1,p) = 1./gamrnd(sp,sc);
    end
end

Оцените следующие средние значения для $β$ , $σ^{2}$ , и $λ$ .

postEstBeta = squeeze(mean(estBeta(:,(burnin+1):end,:),2));
postLambda = squeeze(mean(lambda(:,(burnin+1):end,:),2));

Для каждого $ν$ , отобразите на графике данные и линии регрессии.

figure;
plot(x,y,'o');
h = gca;
xlim = h.XLim';
ylim = h.YLim;
hl = legend('Data');
hold on;
for p = 1:numNu;
    plotY = [ones(2,1) xlim]*postEstBeta(:,p);
    plot(xlim,plotY,'LineWidth',2);
    hl.String{p+1} = sprintf('nu = %f',nu(p));
end
xlabel('x');
ylabel('y');
title('Robust Bayesian Linear Regression')

Низкие значения $ν$ будьте склонны приписывать высокую изменчивость влиятельным выбросам. Поэтому линия регрессии напоминает regARIMA строка. Как $ν$ увеличения, линия ведет себя больше как те из наивного подхода.

Документация